ggplot2 で地図を描くとき sf 以外のオブジェクトも簡単に扱えるようになったので試してください

R GIS ggplot2 sf

少し前、 ggplot2 の開発版に Claus Wilke が半年くらい開発していた pull request がマージされ、sf 以外のオブジェクトを簡単に扱えるようになりました。

しばらくリリースはないんですが、内部的にはけっこう大きな変更だったので、 ggplot2 を地図を描くのに使っている方はぜひ試してみてフィードバックをください（GitHub に直接でも、私に Twitter ででも、r-wakalang ででも）。

CRAN 版（3.3.2）

まずは今の CRAN 版だとなにが不足なのか確認してみましょう。データは何でもいいんですが、適当に地理院地図のベクトルタイルを使います。詳しくは以下のブログ記事を参照。

library(ggplot2)
library(patchwork)

sec <- function(x) {
  1 / cos(x)
}

lonlat2xy <- function(lat_deg, lon_deg, zoom) {
  n <- 2^zoom
  
  x <- (n * (lon_deg + 180)) %/% 360
  lat_rad <- lat_deg * pi / 180
  y <- (n * (1 - log(tan(lat_rad) + sec(lat_rad)) / pi)) %/% 2
  
  list(x = x, y = y)
}

zoom_level <- 5
xy <- lonlat2xy(35, 139, zoom_level)

url <- glue::glue("https://cyberjapandata.gsi.go.jp/xyz/experimental_bvmap/{zoom_level}/{xy$x}/{xy$y}.pbf")
v <- protolite::read_mvt_sf(url)

このデータをプロットして、北緯36度、東経140度に赤い点を打ったり、その周辺に特定の領域にズームしたりしてみましょう。

p <- ggplot() +
  geom_sf(data = v$coastline) +
  theme_minimal()

p1 <- p +
  annotate("point", x = 140, y = 36, colour = "red")


p2 <- p1 +
  coord_sf(xlim = c(139, 141), ylim = c(35, 37))

p1 + p2 + plot_annotation(caption = "地理院地図ベクトルタイルを加工して作成")

これはうまくいきました。このデータの座標系が緯度経度で表されるものだからです。

sf::st_is_longlat(v$coastline)
#> [1] TRUE

一方、UTM 座標系では座標はメートルで表されるので、うまくいきません。EPSG 6691 を指定して試してみましょう。

p3 <- p1 + coord_sf(crs = 6691)
p4 <- p1 + coord_sf(crs = 6691, xlim = c(139, 141), ylim = c(35, 37))

p3 + p4 + plot_annotation(caption = "地理院地図ベクトルタイルを加工して作成")

海の上がズームされてしまいました。 もちろん、この緯度経度を一度 sf のオブジェクトにして、sf::st_transform() すればメートルでの位置を知ることはできますが、 ちょっと点とかテキストを置きたいだけなのにそこまでするのはめんどくさいですよね。

開発版

この煩雑さを解消すべく、 geom_sf() 以外のレイヤーはすべて緯度経度で表された座標だと解釈されるようになりました。 （正確には、 coord_sf() に default_crs という引数に指定された CRS に変換されるようになっていて、このデフォルトが 4326 （WSG 84）なので緯度経度になります。もちろん任意の CRS を指定することが可能です）

なので、同じコードを開発版で実行すると、以下のようになります。

p3 <- p1 + coord_sf(crs = 6691)
p4 <- p1 + coord_sf(crs = 6691, xlim = c(139, 141), ylim = c(35, 37))

p3 + p4 + plot_annotation(caption = "地理院地図ベクトルタイルを加工して作成")

思ったとおりに点を打ったりズームできてそうです。

（重要） geom_map()、coord_map()、coord_quickmap()は非推奨に

これをもって、 geom_map()/coord_map() にできて geom_sf()/coord_sf() にできないことはなくなったはずなので、これらはやがて非推奨（というか superseded ）になります。 もし、geom_sf()/coord_sf() ではできないようなパターンをご存知の方はご一報ください。