sfオブジェクトを描画する方法まとめ: ggplot2::geom_sf()もあるよ

Rの作図パッケージとして人気なggplot2の時期バージョンの2.3.0が間も無くリリースされるそうです。ggplot2は前回の更新が2016年12月末のバージョン2.2.1なので、久しぶりのバージョンアップとなります。

バージョンアップに伴う変更点はこのページを見て欲しいですが、私としてはsfオブジェクトを描画するための新たな関数geom_sf()に触れずには要られません。これまでずっと開発版でのみ利用可能でしたが、やっとCRAN版でも利用できるようになります。（人に紹介するときはGitHubの説明からやったり面倒でもあった。）

... と思っていたらid: yutannihilationさんに先を越されてしまいました。

ggplot2 2.3.0（RC版）を使ってみた - Technically, technophobic.

notchained.hatenablog.com

同じ話題をしても仕方がないので、ここでは、ggplot2に実装されるgeom_sf()の紹介を兼ねて、sfパッケージにより作成されるsfオブジェクトを描画する方法を整理しようと思います。また、sfオブジェクトの可視化に祭して、私の知っているいくつかのtipsを紹介します。(私の持ち合わせる知識の範囲です。他にあればぜひ教えてください。)

公式のvignettesにもsfオブジェクトを描画する方法について整理されています。

https://cran.r-project.org/web/packages/sf/vignettes/sf5.html

sfオブジェクトを描画する際はまず、静的あるいは動的に描画するという目的によって使い分けます。静的な描画はRの標準グラフィックスに利用されるplot()、冒頭のggplot2などです。一方、動的な描画は、興味のある箇所を拡大したり、表示・非表示をユーザが切り替え可能な描画の方法です。具体的にはleaflet、plotlyがこの方法での選択肢となります。

私は、静的な出力の場合はgeom_sf()、動的にズームしたり動かしたい場合はmapview()というような使いわけをしています（少し前はmapviewではなくてleafletでしたが鞍替えしました）。

用意
plot
ggplot2
tmap
leaflet
mapview
plotly

用意

作図用のデータとして、国土数値情報の提供する行政区域のデータと市区町村役場データを用います。このデータはjpndistrictによりRから取得可能です。

library(sf)
library(jpndistrict)

## Loading required package: jpmesh
## This package provide map data is based on the Digital Map 25000(Map Image) published by
## Geospatial Information Authorityof Japan (Approval No.603FY2017 information usage
## <http://www.gsi.go.jp>)

# 岡山県の行政区域データ
# この記事の中の岡山県のポリゴンデータを使った図は、国土地理院長の承認を得て、
# 同院発行の数値地図（国土基本情報）電子国土基本図（地図情報）を使用したものです (承認番号　平28情使、第603号)
pref_33 <- 
  jpn_pref(33)
head(pref_33)

## Simple feature collection with 6 features and 4 fields
## geometry type:  GEOMETRY
## dimension:      XY
## bbox:           xmin: 133.6025 ymin: 34.4173 xmax: 134.169 ymax: 35.30721
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## # A tibble: 6 x 5
##   pref_code prefecture city_code city                             geometry
##   <chr>     <chr>      <chr>     <chr>                      <GEOMETRY [°]>
## 1 33        岡山県     33101     岡山市 北区 POLYGON ((133.9098 34.948, 1…
## 2 33        岡山県     33102     岡山市 中区 MULTIPOLYGON (((133.9747 34.…
## 3 33        岡山県     33103     岡山市 東区 MULTIPOLYGON (((133.9958 34.…
## 4 33        岡山県     33104     岡山市 南区 MULTIPOLYGON (((133.9109 34.…
## 5 33        岡山県     33202     倉敷市      MULTIPOLYGON (((133.6414 34.…
## 6 33        岡山県     33203     津山市      POLYGON ((134.0871 35.30721,…

# 岡山県内の市役所等のポイントデータ
df_office33 <- 
  jpn_admins(33)

## options:        ENCODING=cp932 
## Reading layer `P34-14_33' from data source `/private/var/folders/0x/mb63hycs4k30_7httqyxh2rr0000gn/T/RtmpVmQLFw/P34-14_33_GML/P34-14_33.shp' using driver `ESRI Shapefile'
## Simple feature collection with 105 features and 4 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 133.3315 ymin: 34.35497 xmax: 134.3596 ymax: 35.28274
## epsg (SRID):    NA
## proj4string:    +proj=longlat +ellps=GRS80 +no_defs

head(df_office33)

## Simple feature collection with 6 features and 4 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 133.8795 ymin: 34.65511 xmax: 134.039 ymax: 34.90654
## epsg (SRID):    NA
## proj4string:    +proj=longlat +ellps=GRS80 +no_defs
##   jis_code type           name                    address
## 1    33100    1     岡山市役所        岡山市北区大供1-1-1
## 2    33101    2       御津支所     岡山市北区御津金川1020
## 3    33101    2       建部支所    岡山市北区建部町福渡489
## 4    33103    2       瀬戸支所     岡山市東区瀬戸町瀬戸45
## 5    33101    2     鶴田連絡所 岡山市北区建部町角石谷2063
## 6    33101    1 岡山市北区役所        岡山市北区大供1-1-1
##                    geometry
## 1 POINT (133.9196 34.65511)
## 2  POINT (133.9354 34.7992)
## 3 POINT (133.9038 34.86932)
## 4  POINT (134.039 34.73388)
## 5 POINT (133.8795 34.90654)
## 6 POINT (133.9196 34.65511)

plot

sfパッケージには、作図関数のplot()の総称関数が用意されています。そのため、sfが生成する、sfg (simple feature geometry)、sfc (simple feature geometry list-column)、sf (simple feature)のいずれかのクラスに属するsfオブジェクトであれば、plot()を適用することで、難なく地物を描画することが可能です。

par(mfrow = c(1, 2))
# sfg
st_point(0:1) %>% 
  plot()
# sfc
st_sfc(st_point(0:1), 
       st_point(2:1)) %>% 
  plot(col = "red")

f:id:u_ribo:20180528070948p:plain

sfcを含んだsfオブジェクト（データフレームの変数 + geometryのようなもの)を描画すると全ての変数を対象に描画が行われます。

# 既定では全ての変数について描画を行う
plot(pref_33)

f:id:u_ribo:20180528070609p:plain

これは一旦の確認には良いですが、最終的な出力では一変数だけを出力することが多いかと思います。その際はst_geometry()をかまして地物情報だけを抽出し、sfcとして描画するか、特定列を参照した状態で描画を行うと良いでしょう。

# 地物情報だけを描画
st_geometry(pref_33) %>% 
  plot()

# 特定列だけの描画
# 凡例の表示のため、余白を調節します
par(oma = c(3, 3, 3, 6.5), mar = c(2, 1, 1, 2))
pref_33["city"] %>% 
  plot(family = "IPAPGothic")

f:id:u_ribo:20180528071544p:plain

ggplot2

次に今回の目玉となるggplot2でのgeom_sf()を使う方法です。基本的なことは冒頭にあげたyutannihilationさんの記事に書かれています。

library(ggplot2)

ggplot(pref_33) +
  geom_sf()

f:id:u_ribo:20180528071604p:plain

ggplot2に詳しい方ならご存知かと存じますが、ggplot2 では描画したいデータをggplot()の第一引数dataに指定します（ここでは岡山県の行政区域データ）。それを+で連結させ、どのように描画するかという際にgeom_*()を利用します。今回実装されるgeom_sf()はsfオブジェクトを描画する関数で、これまで散布図ならgeom_point()、棒グラフならgeom_bar()のように使い分けてきたように、sfオブジェクトであればgeom_sf()という塩梅となります。

ggplot(data = )のgeom_sf(data = )に渡す方法もあります。これは、描画したいオブジェクトが複数（例えばポリゴンとポイント）ある際に、それぞれのgeom_sf()で対象のデータを指定する時に、どのデータをどのにレイヤーに渡すかを区別させるために利用されます。

ggplot() +
  geom_sf(data = pref_33, aes(fill = city_code)) +
  geom_sf(data = df_office33, color = "white")

f:id:u_ribo:20180528071619p:plain

次のコードでも同じ図が得られますが、geom_sf()にどのデータが渡っているかがわかりにくいので、上記の方法が好みです。

ggplot(pref_33) +
  geom_sf(aes(fill = city_code)) +
  geom_sf(data = df_office33, color = "white")

より地図っぽく描画するには nozmaさんの記事にあるようにcoords_sf()やtheme_void()を使うと良いでしょう。

ggplot() +
  geom_sf(data = pref_33, aes(fill = city_code)) +
  # 緯度経度の線、ラベルを描画しない
  coord_sf(datum = NA) +
  # 背景色を白にする
  theme_void()

f:id:u_ribo:20180528071723p:plain

tmap

tmapは主題図を描くために適したパッケージですが、こちらもsfオブジェクトに対応しています。描画のための関数はqtm()およびtm_shape()です。

library(tmap)

qtm(pref_33)

f:id:u_ribo:20180528071956p:plain

tm_shape()ではmodeと呼ばれる機能があり、それを切り替えることで静的・動的な描画の切り替えが可能です。

# 結果は上の図と同じなので省略
tmap_mode("plot")

tm_shape(pref_33, projection = "longlat") + 
    tm_polygons()

インタラクティブに動かせる図を描画するにはモードをviewに切り替えます。動作の様子をgif画像で表示します。

tmap_mode("view")

tm_shape(pref_33) + 
  tm_fill("city_code", palette = sf.colors(20))

leaflet

leafletはJavaScript製の地図描画ライブラリですが、Rに実装されたleafletパッケージはsfオブジェクトを描画し、ユーザが任意の箇所を拡大したり表示箇所を移動することが可能です。tmapやmapviewの基礎として利用されているライブラリでもあります。

leafletでsfを表示するには、対象とするsfのgeometryの種類を把握した上で適した関数に渡す必要があります。すなわちポイントの場合は、addCircles()、ラインはaddPolylines()、ポリゴンはaddPolygons()といった具合です。

library(leaflet)

leaflet() %>% 
  addTiles() %>% 
  # ポイントの表示なのでaddCircles()
  # 日本のへそ公園の位置を表示
  addCircles(data = st_point(c(135, 35.0)) %>% 
               st_sfc() %>% 
               st_sf(crs = 4326))

f:id:u_ribo:20180528072339p:plain

leaflet() %>% 
  addTiles() %>% 
  # ポリゴンの表示
  addPolygons(data = pref_33)

f:id:u_ribo:20180528072414p:plain

mapview

ggplot2のようにオブジェクトを重ねていくことができますが、種類に応じて関数を分ける必要がある点や、ベースタイルの表示や表示の切り替えに別の関数を利用しなければならないという点でちょっと面倒です。そこで、もっとサクッと表示したいという時にmapviewが便利です。

ggplot2のように複数のオブジェクトを+演算子で重ねていくことができ、地物の種類も関係なくmapview()で済むので楽チンです（spオブジェクトにも対応します）。

library(mapview)

mapview(pref_33)

地図描画のためのオプションも豊富で一押しです。

mapview(pref_33, 
        zcol = "city_code",
        legend = FALSE) +
  mapview(st_point(c(133.9196, 34.65511)),
          color = "red", 
          col.regions = "red")

plotly

最後に、leafletライブラリとは別のアプローチでsfを動的に可視化する方法としてplotlyの例です。こちらはまだ開発版の機能となりますが、plotlyの開発ブログ記事に書書かれているように、plotlyの特徴を活かした3Dでの表示や他の種類のグラフと組み合わせが可能なようです。

こちらは開発版の機能ですので、将来実装が変更される可能性があります。

remotes::install_github("ropensci/plotly")
library(plotly)

nc <- 
  sf::st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)
plot_geo(nc)