ggplot2::geom_sf()を使ってみる - Technically, technophobic.

ggplot2 2.3.0全体の変更点については以下に書きましたが、その中でも目玉の一つgeom_sf()について軽く紹介します。

追記： uriboがもっと網羅的な記事を書いたみたいなので、こっちを読んどけばばっちりです！

sfパッケージとは

sfパッケージは、GISデータをRで扱うためのパッケージです。 dplyrやtidyr、leafletなどのモダンなパッケージに対応していて、tidyなエコシステムの中で使うことができます。

しかし、「モダンなパッケージ」と言いましたが、CRAN版のggplot2ではまだsfをプロットすることはできません。

ggplot2のsf対応

実は、geom_sf()がggplot2に追加されたのはもう1年半ほども前です。そんなわけで、GitHub版のggplot2を使っている人にはもうおなじみの機能なんですが、この1年半ggplot2のリリースがなかったために（なんと前回のリリースは2016年12月...）まだCRAN版には入っていないのです。

`geom_sf()`

そのsfを扱うためにggplot2に爆誕した機能というのが、geom_sf()です。

これは、他のgeom_*()と比べてかなり特殊です。渡されるsfオブジェクトの種類によって、点でも線でも面でも、なんでも描いてくれます。また、それでいてgeom_point()やgeom_line()やgeom_polygon()を扱うときと同じように、fillやcolour、sizeといったaesに自在に変数をマッピングすることができます。

library(ggplot2)

# MULTIPOLYGONのデータ
nc <- sf::st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)

# POINTのデータ
data(meuse, package = "sp")
sp::coordinates(meuse) <- ~x+y
sp::proj4string(meuse) <- sp::CRS("+init=epsg:28992")

p1 <- ggplot(nc) +
  geom_sf(aes(fill = AREA))

p2 <- ggplot(sf::st_as_sf(meuse)) +
  geom_sf(aes(colour = cadmium))

gridExtra::grid.arrange(p1, p2, ncol = 1)

凡例

ただし、sfのプロットが点になるのか線になるのか面になるのかは、描いて見ないとわかりません。その結果を凡例を描く処理に伝える方法がないらしく、不正確なものになってしまいます。

具体的にはこんな感じです。

ggplot(sf::st_as_sf(meuse)) +
  geom_sf(aes(size = cadmium))

これはggplot2のアーキテクチャの問題らしく、解決する見込みはないそうです。

ワークアラウンドとしてshow.legendという引数が用意されているのでこれを使いましょう。デフォルトだとポリゴン用の凡例が表示されていますが、これを"line"や"point"を指定することで凡例の種類を変えられます。

ggplot(sf::st_as_sf(meuse)) +
  geom_sf(aes(size = cadmium), show.legend = "point")

ラベル

地図を描くと、地名をテキストで示したくなることもあるでしょう。しかし、それはオプションが複雑になりすぎるという理由でgeom_sf()には実装されませんでした。

なので、こんな感じで自分でX座標とY座標を取り出したデータをつくってgeom_text()あたりに渡す必要があります。

library(ggplot2)

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

# 投影座標系に
nc_3857 <- sf::st_transform(nc, 3857)
# 重心を取り出す
centroids <- sf::st_point_on_surface(nc_3857$geometry)
# 座標を取り出す（X, Yという名前の列になる）
xy <- sf::st_coordinates(centroids)
# 元のデータと合わせる
nc_3857 <- cbind(nc_3857, xy)

ggplot(nc_3857) +
  geom_sf() +
  geom_text(aes(X, Y, label = NAME))

（例が汚くてすいません...。たぶんggrepelを使えばいい感じになるはず）

CRS

座標系は自分で合わせなくても、geom_sf()、というかcoord_sf()がすべてのレイヤーの座標系を一致させてくれます。ヘルプの例に載ってるのはこんな感じ。ぴったり一致してるのがわかるでしょうか。

nc <- sf::st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)
nc_3857 <- sf::st_transform(nc, "+init=epsg:3857")

ggplot() +
  geom_sf(data = nc, colour = "black", fill = NA, size = 3) +
  geom_sf(data = nc_3857, colour = "white", fill = NA) +
  theme_minimal()

これもヘルプに出ている例。

library(maps)

world1 <- sf::st_as_sf(map('world', plot = FALSE, fill = TRUE))
world2 <- sf::st_transform(
  world1,
  "+proj=laea +y_0=0 +lon_0=155 +lat_0=-90 +ellps=WGS84 +no_defs"
)
ggplot() + geom_sf(data = world2)

感想

という感じで、地図を描くのがだいぶお手軽になりそうです。と、GIS素人的には思ったんですが、どうなんでしょう...？あんまりその辺ドメイン知識がないので、詳しい方のツッコミお待ちしています。