第13章 地図のグラフ

library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.2     ✓ purrr   0.3.4
## ✓ tibble  3.0.4     ✓ dplyr   1.0.2
## ✓ tidyr   1.1.2     ✓ stringr 1.4.0
## ✓ readr   1.4.0     ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(maps)
## 
## Attaching package: 'maps'
## The following object is masked from 'package:purrr':
## 
##     map
library(jpndistrict)
## This package provide map data is based on the Digital Map 25000 (Map
## Image) published by Geospatial Information Authority of Japan (Approval
## No.603FY2017 information usage <https://www.gsi.go.jp>)
library(countrycode)
  • mapsは世界地図情報のパッケージです。
    • 読み込むとコンフリクトを起こしてpurrr::map()が使えなくなるので気をつけて下さいー
  • jpndistrictは都道府県地図情報のパッケージです。
  • countrycodeは国名を統一するパッケージです。

13.1 地図

13.1.1 世界地図

まず、世界地図のデータを呼び出します。

world <- map_data("world")
  • mapパッケージには世界、アメリカ、フランス、イタリア、ニュージーランドのデータが入っています。

こんな感じで緯度、軽度、グループ、地域などの情報が入っています。

head(world)

境界線が緯度経度でgeom_polygon()で描くことができます。

world %>% 
  ggplot(aes(x = long, y = lat, group = group)) +
  geom_polygon(fill = "lightgray", colour = "black", size = 0.1)

  • 国境を描くために、groupで指定しています。

13.1.2 日本地図

もちろん、filter()でデータを限定すれば、一部の国の地図を描くこともできます。

world %>% 
  filter(region == "Japan") %>% 
  ggplot(aes(x = long, y = lat, group = group)) +
  geom_polygon(fill = "lightgray", colour = "black", size = 0.1)

13.1.3 都道府県地図

残念ながら、mapパッケージには国レベルの情報しかありません。

world %>% 
  arrange(region) %>% 
  pull(region) %>% 
  unique() %>% 
  head()
## [1] "Afghanistan"    "Albania"        "Algeria"        "American Samoa"
## [5] "Andorra"        "Angola"

都道府県の地図を描きたい場合はjpndistrictというパッケージを使用します。 パッケージ付属のjpnprefsに都道府県一覧が存在します。

head(jpnprefs)

都道府県データを引っ張るにはjpn_pref()に該当する都道府県の数字を入力します。

kyoto <- jpn_pref(26)
## Linking to GEOS 3.6.2, GDAL 2.2.3, PROJ 4.9.3
head(kyoto)

今回は市区町村の境界線がsimple features (sf)と呼ばれる規格なので、geom_sf()を使います。

kyoto %>% 
  ggplot() + 
  geom_sf()

13.2 地域の違い

地域ごとの違いを見るために、色を塗り分けたいともいます。 tidyrの中にpopulationというWHOの人口データがあります。

head(population)

2013年時点の人口データを取り出します。

pop2013 <- population %>% 
  filter(year == 2013) %>% 
  rename(region = country)
  • populationでは国名がcountryですが、worldではregionなので、rename()で後者に合わせます。

先程のデータとworldをマージして、aes()の中のfillで塗り分けする変数を指定します。 それぞれのデータセットで国名が異なっているので、coutnrycodeというパッケージでISOコードに合わせます。

pop2013 <- pop2013 %>% 
  mutate(iso3c = countrycode(sourcevar = region, origin = "country.name", destination = "iso3c"))

world <- world %>% 
  mutate(iso3c = countrycode(sourcevar = region, origin = "country.name", destination = "iso3c"))
## Warning: Problem with `mutate()` input `iso3c`.
## ℹ Some values were not matched unambiguously: Ascension Island, Azores, Barbuda, Bonaire, Canary Islands, Chagos Archipelago, Grenadines, Heard Island, Kosovo, Madeira Islands, Micronesia, Saba, Saint Martin, Siachen Glacier, Sint Eustatius, Virgin Islands
## 
## ℹ Input `iso3c` is `countrycode(sourcevar = region, origin = "country.name", destination = "iso3c")`.
## Warning in countrycode(sourcevar = region, origin = "country.name", destination = "iso3c"): Some values were not matched unambiguously: Ascension Island, Azores, Barbuda, Bonaire, Canary Islands, Chagos Archipelago, Grenadines, Heard Island, Kosovo, Madeira Islands, Micronesia, Saba, Saint Martin, Siachen Glacier, Sint Eustatius, Virgin Islands
world %>% 
  filter(is.na(iso3c) == TRUE) %>% 
  select(region) %>% 
  distinct()
left_join(world, pop2013, by = "iso3c") %>% 
  ggplot() + 
  geom_polygon(aes(x = long, y = lat, group = group, fill = log(population)),
               colour = "black", size = 0.1) + 
  scale_fill_continuous(name = "Population (log)", low = "lightblue", high = "darkblue")

geom_sf()でも同じです。

kyoto %>% 
  ggplot() + 
  geom_sf(aes(fill = city), show.legend = FALSE)

13.3 空間の分布

地図上に点を打つ場合はgeom_point()を使うだけです。 mapsパッケージの中にworld.citiesというデータセットがあるので、各国の都市を表示します。

world %>% 
  ggplot() + 
  geom_polygon(aes(x = long, y = lat, group = group),
               fill = "lightgray", colour = "black", size = 0.1) + 
  geom_point(aes(x = long, y = lat, alpha = pop/10000),
             size = 0.1, data = world.cities)

  • ただし、地図データとは別のデータセットのはずなので、別途指定する必要があります。

geom_sf()の場合も同様にできます。

kyoto %>% 
  ggplot() + 
  geom_sf() + 
  geom_point(aes(x = capital_longitude, y = capital_latitude),
             colour = "red", data = jpnprefs %>% filter(prefecture == "京都府"))

また、geom_density_2d()で2次元の分布を描くこともできます。 試しに、都市の分布を示します。

world %>% 
  ggplot() + 
  geom_polygon(aes(x = long, y = lat, group = group),
               fill = "lightgray", colour = "black", size = 0.1) + 
  geom_density_2d(aes(x = long, y = lat),
                  data = world.cities)

  • 一般的には、geom_contour()で等高線を描くことができます。

色をつけることもできます。

world %>% 
  ggplot() + 
  geom_polygon(aes(x = long, y = lat, group = group),
               fill = "lightgray", colour = "black", size = 0.1) +
  stat_density2d(aes(x = long, y = lat, fill = stat(level)), 
                 geom = "polygon", alpha = 0.5) 

13.4 データの繋がり

13.4.1 二点間の繋がりを見る

データポイントを打った後は、それぞれを線で繋ぎたいのが人情だと思います。 geom_segment()を使うと、始点と終点のx軸とy軸の値を決めると線を引いてくれます。

まず、準備としてOpenFlightsから世界の空港と航路のデータを取得し、結合させます。

airports <- read_csv("data/airports.dat", col_names = FALSE) %>% 
  select(id = X1, lat = X7, long = X8)
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   X1 = col_double(),
##   X2 = col_character(),
##   X3 = col_character(),
##   X4 = col_character(),
##   X5 = col_character(),
##   X6 = col_character(),
##   X7 = col_double(),
##   X8 = col_double(),
##   X9 = col_double(),
##   X10 = col_double(),
##   X11 = col_character(),
##   X12 = col_character(),
##   X13 = col_character(),
##   X14 = col_character()
## )
## Warning: 353 parsing failures.
##  row col expected actual                file
## 6982 X10 a double    \N 'data/airports.dat'
## 6983 X10 a double    \N 'data/airports.dat'
## 6984 X10 a double    \N 'data/airports.dat'
## 6985 X10 a double    \N 'data/airports.dat'
## 6986 X10 a double    \N 'data/airports.dat'
## .... ... ........ ...... ...................
## See problems(...) for more details.
route <- read_csv("data/routes.dat", col_names = FALSE) %>% 
  select(source_id = X4, dest_id = X6) %>% 
  mutate_all(as.numeric)
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   X1 = col_character(),
##   X2 = col_character(),
##   X3 = col_character(),
##   X4 = col_character(),
##   X5 = col_character(),
##   X6 = col_character(),
##   X7 = col_character(),
##   X8 = col_double(),
##   X9 = col_character()
## )
## Warning: Problem with `mutate()` input `source_id`.
## ℹ NAs introduced by coercion
## ℹ Input `source_id` is `.Primitive("as.double")(source_id)`.
## Warning in mask$eval_all_mutate(dots[[i]]): NAs introduced by coercion
## Warning: Problem with `mutate()` input `dest_id`.
## ℹ NAs introduced by coercion
## ℹ Input `dest_id` is `.Primitive("as.double")(dest_id)`.
## Warning in mask$eval_all_mutate(dots[[i]]): NAs introduced by coercion
df_flight <- route %>% 
  left_join(airports %>% 
              rename(source_id = id, source_lat = lat, source_long = long),
            by = "source_id") %>% 
  left_join(airports %>% 
              rename(dest_id = id,dest_lat = lat, dest_long = long),
            by = "dest_id") %>% 
  filter(!(source_long == dest_long & source_lat == dest_lat))
  • source_longsource_latが出発地の緯度経度、dest_longdest_latが到着地の緯度経度を示しています。

それを世界地図の上に重ねてプロットします。

world %>% 
  ggplot() +
  geom_polygon(aes(x = long, y = lat, group = group), 
               fill = "lightgray", colour = "black", size = 0.1) + 
  geom_segment(aes(x = source_long, y = source_lat, xend = dest_long, yend = dest_lat),
               size = 0.1, alpha = 0.1, 
               data = df_flight)

  • 太平洋を横断する路線も西回りになってしまうのは課題です……

geom_curve()で曲線になります。

world %>% 
  ggplot() +
  geom_polygon(aes(x = long, y = lat, group = group), 
               fill = "lightgray", colour = "black", size = 0.1) + 
  geom_curve(aes(x = source_long, y = source_lat, xend = dest_long, yend = dest_lat),
               size = 0.1, alpha = 0.1, 
               data = df_flight)

13.4.2 経路を見る

二点間を繋ぐのではなく、一連の経路を描きたい場合は、この方法では非効率的です。 例として、2019年の台風の経路をプロットするために、デジタル台風のサイト()からデータを収集します。

library(rvest)
## Loading required package: xml2
## 
## Attaching package: 'rvest'
## The following object is masked from 'package:purrr':
## 
##     pluck
## The following object is masked from 'package:readr':
## 
##     guess_encoding
df_typhoon <- tibble()
for (i in formatC(1:28, width = 2, flag = "0")) {
  url <- str_glue("http://agora.ex.nii.ac.jp/digital-typhoon/summary/wnp/l/2019{i}.html.ja")
  df_typhoon <- bind_rows(df_typhoon,
                          url %>% 
                            read_html() %>% 
                            html_node("table.TRACKINFO") %>% 
                            html_table(fill = TRUE) %>% 
                            select(long = "経度", lat = "緯度", hpa = "中心気圧 (hPa)") %>%
                            mutate(id = i))
}

このように時系列に緯度経度が並んでいます。

head(df_typhoon)

経路の場合はgeom_path()を使います。

world %>% 
  ggplot() +
  geom_polygon(aes(x = long, y = lat, group = group), 
               fill = "lightgray", colour = "black", size = 0.1) + 
  geom_path(aes(x = long, y = lat, colour = hpa, group = id),
            data = df_typhoon) + 
  scale_colour_continuous(low = "lightblue", high = "darkblue")

  • 複数の経路がある場合はgroupで指定するのを忘れずに。