Go to Qiita Advent Calendar 2024 Top
LoginSignup
18
18

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

@uri(Shinya Uryu)in株式会社ホクソエム

【まとめ用】ggplot2を使って地図を書く

Posted at

Rを使って地図を表示することがしばしばある。その度にやり方を忘れていて、調べ直すことが多い。Google先生に尋ねると、当然のごとく結果が変わっている。時間が経つと、より評価され・より良い内容の記事が検索のトップに並ぶからだ。なので毎回方法が変わってしまう。自分なりの手法を確立しておきたい。

Rを用いた地理情報データとデータマッピングの流れ

まずは先行して書かれている既存のブログ記事などを参照に、Rを使って地図を描画する手順を整理してみる。大まかな流れとして、

  1. 必要なファイル（Shapefile）を読み込む
  2. 描画用のパッケージ関数を使って地図を描く
  3. 目的に応じてデータを地図に追加する（マッピング）

となっている。順を追って手を動かしてみよう。なお個人的に、地図をプロットする際にはこうしたい、という願望があるので、

  1. Shapefileは使いたくない。GeoJSONかTopoJSONファイルを使用したい

    • Shapefileの仕様が謎。複数のファイルが必要だったりして複雑
    • JSONベースのファイルは軽量かつ読みやすい

  2. プロットには**{ggplot2}**を使用したい

  3. 複雑ではなく、限りなくシンプルなコードにしたい

これらをできるだけ達成できるよう頑張る。

地図描画のためのパッケージあれこれ

適宜インストール。

library(magrittr)
library(sp)
library(ggplot2)
library(rgdal)
library(broom)
library(ggmap)
library(viridis)

1 必要なファイルを用意し、Rで扱えるようにする

幾つかの方法があるが、ここでは**{rgdal}**のreadOGR関数を利用してファイルを読み込む。readOGRでは、ベクター形式のファイルを読み込むことができる。readOGRの引数をみると、

args(readOGR) %>% as.list() %>%  names()

##  [1] "dsn"                       "layer"                    
##  [3] "verbose"                   "p4s"                      
##  [5] "stringsAsFactors"          "drop_unsupported_fields"  
##  [7] "input_field_name_encoding" "pointDropZ"               
##  [9] "dropNULLGeometries"        "useC"                     
## [11] "disambiguateFIDs"          "addCommentsToPolygons"    
## [13] "encoding"                  "use_iconv"                
## [15] "swapAxisOrder"             "require_geomType"         
## [17] "integer64"                 ""

が用意されている。このうち重要になってくるのは、描画するソースを定義するdsnとlayerである。最低限、これだけを指定すれば読み込んでくれる。今回は**{rgdal}**に内臓されているデータを使うことにする。次のコマンドを実行して出力されるディレクトリにデモ用のShapefileがある。

path <- vector()
system.file("vectors", package = "rgdal") %>% {
  path <<- .
  print(.)
}

## [1] "/Library/Frameworks/R.framework/Versions/3.2/Resources/library/rgdal/vectors"

次のようなディレクトリ構造になっている（一部のみ出力）。

list.files(path) %>% grep(pattern = "^trin_inca_pl03", x = ., value = TRUE)

## [1] "trin_inca_pl03.dbf" "trin_inca_pl03.shp" "trin_inca_pl03.shx"

trin_inca_pl03.shpに付属して、幾つかのファイルが用意されている。Shapefileを出力させる場合、これらのファイルも必要になる。では、このtrin_inca_pl03ファイルをRに取り込んでみる。

map <- data.frame()
readOGR(dsn = path, layer = "trin_inca_pl03") %>% {
  map <<- .
  print(class(.))
  str(., max.level = 2)
}

## OGR data source with driver: ESRI Shapefile 
## Source: "/Library/Frameworks/R.framework/Versions/3.2/Resources/library/rgdal/vectors", layer: "trin_inca_pl03"
## with 3 features
## It has 12 fields
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  2 obs. of  12 variables:
##   ..@ polygons   :List of 2
##   ..@ plotOrder  : int [1:2] 1 2
##   ..@ bbox       : num [1:2, 1:2] -87.1 5.5 -79 10.1
##   .. ..- attr(*, "dimnames")=List of 2
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot

readOGRで読み込んだオブジェクトは**{sp}**パッケージで定義されるSpatialPointsDataFrameというものになる。S4メソッドで定義されているので、各オブジェクトへのアクセスは@スロットを利用する。strで示したように、SpatialPointsDataFrameは５つのスロットを持つ。@dataは属性値を持つデータフレームであり、@coordsは位置座標の行列オブジェクトである。

dplyr::tbl_df(map@data)

## Source: local data frame [2 x 12]
## 
##           ENGL_NAME      SCI_NAME      AUTHORITY       FAMILY PRESENCE
##              (fctr)        (fctr)         (fctr)       (fctr)    (int)
## 1 Wandering Tattler Tringa incana (Gmelin, 1789) Scolopacidae        1
## 2 Wandering Tattler Tringa incana (Gmelin, 1789) Scolopacidae        1
## Variables not shown: ORIGIN (int), COMPILER (fctr), SCALE (fctr), TAX_COM
##   (fctr), DIST_COM (fctr), REFERENCES (fctr), REVIEWERS (fctr)

地図を描画する場合、元になるファイルが必要になるが、多くの場合Shapefileを使うことになるだろう。しかし、readOGRではShapefile以外のJSON(GeoJSON, TopoJSON)形式のファイルでも読み込むことが可能である。いずれの方法もShapefileの時と同じである。

# topojson
map <- data.frame()
# 国土交通省 国土地理院 ２０万分１地勢図図を元にしたファイルを用意する
file.path <- "https://gist.github.com/uribo/b09d642351c03dd975aa/raw/6cd2f7e70922b24f89b442c96c001dc3ab794317/japan_one_twenty_map.topojson"
readOGR(dsn              = file.path, 
        layer            = "OGRGeoJSON", 
        stringsAsFactors = FALSE) %>% {
          map <<- .
          print(class(.))
          str(., max.level = 2)
}

## OGR data source with driver: GeoJSON 
## Source: "https://gist.github.com/uribo/b09d642351c03dd975aa/raw/6cd2f7e70922b24f89b442c96c001dc3ab794317/japan_one_twenty_map.topojson", layer: "OGRGeoJSON"
## with 130 features
## It has 2 fields
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  130 obs. of  2 variables:
##   ..@ polygons   :List of 130
##   .. .. [list output truncated]
##   ..@ plotOrder  : int [1:130] 75 78 59 41 79 93 8 11 61 31 ...
##   ..@ bbox       : num [1:2, 1:2] 123 24 149 46
##   .. ..- attr(*, "dimnames")=List of 2
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot

GeoJSONあるいはTopoJSONの場合、readOGRのlayer引数に渡す値がOGRGeoJSONとなっているが、これは

ogrListLayers(dsn = file.path)[1]

## [1] "OGRGeoJSON"

で出力されるものである。出力したいlayers名が不明な場合、こちらで確認すると良いだろう。

対象のファイルがGeoJSONあるいはTopoJSON形式である場合、**{rgdal}の代わりに{geojsonio}**を利用することも可能である。

geojsonio::topojson_read(file.path) %>% {
  print(class(.))
  str(., max.level = 2)
}

2 地図をプロットする

readOGRにより、SpatialPolygonsDataFrameオブジェクトを作成したら、ひとまず地図がプロットできるようになる。plot(map)のように、plot関数にオブジェクト名を指定するだけで地図が描画される。ggplotによるプロットは次のようにすれば実行されるが、これはポリゴンとして描画しているため解像度が荒い。

## # ggplot2::ggplot による描画: よくない例。荒い。data.frame, SpatialPolygonsDataFrameどちらでもおk
ggplot() +  geom_polygon(data = map, 
                          aes(x = long, y = lat, group = id))

unnamed-chunk-11-1.png

{ggplot2}を使って描画させる場合は、SpatialPolygonsDataFrameをdata.frameに変換する必要がある。その際、{broom} のtidyを利用すると便利である（これまではggplot2::fortifyが使われてきたが、fortify は将来廃止予定とのことであり、 代わりに{broom}パッケージを使うことをHadleyは推奨している）。

broom::tidy(map) %>% {
  map <<- .
  print(class(.))
  dplyr::glimpse(.)
}

## [1] "data.frame"
## Observations: 650
## Variables: 7
## $ long  (dbl) 148, 149, 149, 148, 148, 148, 149, 149, 148, 148, 147, 1...
## $ lat   (dbl) 46.00000, 46.00000, 45.33333, 45.33333, 46.00000, 45.333...
## $ order (int) 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4,...
## $ hole  (lgl) FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, ...
## $ piece (fctr) 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
## $ group (fctr) 1.1, 1.1, 1.1, 1.1, 1.1, 2.1, 2.1, 2.1, 2.1, 2.1, 3.1, ...
## $ id    (chr) "1", "1", "1", "1", "1", "2", "2", "2", "2", "2", "3", "...

tidy処理を施したオブジェクトは共通の列名を持つデータフレームとなる。このオブジェクトを用いてggplotの上に地図をプロットしていく。

ggplot() + 
  geom_map(data = map, 
           map = map,
           aes(x = long, y = lat, map_id = id), 
           fill = "white", color = "black")

ggplotのdata引数及びgeom_mapのmapに与えるデータフレームはtidy処理を施したオブジェクトである。審美的属性に与える要素も共通である。これが基本型となるが、これだと地図が塗りつぶしされていたり、プロットが正方形になっていたりと、地図として格好が悪いので、次のように幾つかの属性を付加しておくと良い。

ggplot() +
  geom_map(data = map, 
           map = map,
           aes(x = long, y = lat, map_id = id),
           color = "black", fill = "white", # 境界線を黒く、領域を白で塗りつぶす
           size = 0.5) +　
  coord_map(projection = "mercator") -> plot_map # 投影法を指定（メルカトル図法）
# ggthemes::theme_map() も地図描画用のggplotテーマなので、指定しておくと良いかもしれない。
# 背景が白くなり、軸などが表示されなくなる

plot_map

gg_map_basic-1.png

plot_mapとして、再帰的に地図を描画できるようにした。ここから必要なデータを重ねたりしていく。

3 マッピング

いよいよ本番。自分のデータを地図上に追加していく。まず、データを用意するところから始め、マッピングを進めていく。地図上にデータを重ねていくマッピングの手法として、ラベルの出力、ポイントの表示、値や水準による塗り分けがあると考えられるので、それぞれについて日本の地図を例を元にして説明する。

データの用意

マッピングを行う際に必要になるのは、どの座標（緯度と経度）にどのような値（文字や種類）を出力するか、塗り分けの基準となる値や水準といったものなので、これらを含んだデータフレームを用意しておく。

df_jp <- data.frame(address          = c("Tokyo, Japan", 
                                         "Okayama, Japan", 
                                         "Sapporo, Japan",
                                         "Okinawa, Japan"),
                    id               = c(56, 90, 26, 125),
                    population       = c(8945695, 709584, 1913545, 315954),
                    # https://ja.wikipedia.org/wiki/都道府県庁所在地と政令指定都市の人口順位 から人口を入力
                    stringsAsFactors = FALSE) %>% 
  mutate_geocode(address) %>% # addressを元に、緯度経度を求める
  dplyr::mutate(address = gsub(pattern = "[[:punct:]] Japan$", replacement = "", address))

dplyr::tbl_df(df_jp)

## Source: local data frame [4 x 5]
## 
##   address    id population      lon      lat
##     (chr) (dbl)      (dbl)    (dbl)    (dbl)
## 1   Tokyo    56    8945695 139.6917 35.68950
## 2 Okayama    90     709584 133.9350 34.66167
## 3 Sapporo    26    1913545 141.3469 43.06417
## 4 Okinawa   125     315954 127.8014 26.33583

ラベル及びポイントの表示

緯度・経度、出力する文字列はそれぞれ、df_jp$lat, df_jp$lon, df_jp(addressである。また、ポイントの大きさを人口 df_jp)population`を元にして表示してみる。

plot_map + 
  # ポイントの出力
  geom_point(data = df_jp, 
             aes(x = lon, y = lat),
             fill = "skyblue",
             shape = 21, alpha = 0.8, size = 2 * log(df_jp$population) / 10) + 
  geom_text(data = df_jp, 
            aes(label = address, x = lon + 0.2, y = lat + -0.4), 
            size = 3, 
            color = "tomato",
            angle = 20)

unnamed-chunk-15-1.png

塗り分け地図（コロプレス図）

塗り分けは少し複雑になる。ベースとしている地図での塗り分けの基準と、データフレームでの基準を一致させていないと、塗り分けされる箇所が変わってしまう。

plot_map + 
  geom_map(data = df_jp, 
           map = map,
           aes(fill = population, map_id = id, stat = "identity")) +
  scale_fill_viridis()

unnamed-chunk-16-1.png

こうした塗り分け地図を描画する専門のパッケージもある。

https://github.com/trulia/choroplethr

おまけ: 拡大・縮小できる地図

実際に動作しているのはこちらを参考にしてほしい。

RPubs - Choropleth Map of Japan with svgPanZoom, viridis and lineworkmaps packges.

library(svgPanZoom)
library(SVGAnnotation)

svgPanZoom((plot_map), 
           controlIconsEnabled = TRUE)

次回は応用編とさまざまな地図描画手法の紹介です

参考

18
18
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
18
18

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?