メモ:fgdrパッケージで基盤地図情報のDEMデータをstarsとして読み込みrayshaderする

fgdrパッケージがCRANに登録されたというツイートを見かけて、そういえば試してないなと思って。

詳しい使い方はこっちの記事が詳しいです。

データの入手

まず基盤地図情報のデータをここからダウンロードします(データは利用規約に従って使いましょう)

今回は、5mメッシュの標高データを使います。

データの読み込み

ダウンロードしたデータはZIPファイルになっていて、展開すると複数のXMLファイルが入っています。 read_fgd_dem()は1度に1つのファイルしか読めないので、map()とかlapply()とかで全ファイルに対してread_fgd_dem()を繰り返します。 ちなみに、read_fgd_dem()は読み込み形式としてdata.frame、raster、starsをサポートしていますが、今回はstarsを使ってみたいのでreturn_class = "stars"を指定します。

library(fgdr)
library(purrr)

# ZIPファイルを展開
tmp <- tempfile()
dir.create(tmp)
unzip("~/Downloads/FG-GML-5339-45-DEM5A.zip", exdir = tmp)

# 展開したファイルをリストアップ
xml_files <- list.files(tmp, full.names = TRUE)

# 今回は5mメッシュなのでresolutionは5
r <- map(xml_files, ~ read_fgd_dem(., resolution = 5, return_class = "stars"))

データの結合

すべてのデータをstars形式で読み込むことができたので、今度はこれを1つのラスタに合成します。 ここの操作がわからず(Geocomputation with Rにもrasterの操作しか書いてなくて...)半日くらい悩んでたんですが、stars::st_mosaic()が正解みたいです。 ラスタデータの操作まったく理解できてないので、変なところがあればご指摘ください。。

library(stars)
x <- do.call(st_mosaic, r)

rayshader

rayshaderを使うにはmatrix形式のデータが必要です。rasterの場合はこんな感じで変換するみたいですが:

#And convert it to a matrix:
elmat = matrix(raster::extract(localtif, raster::extent(localtif), buffer = 1000),
               nrow = ncol(localtif), ncol = nrow(localtif))

starsの場合は、内部的にmatrixでデータを持っているのでそのまま使えるっぽいです。 たぶんメッシュが均等でない場合はもうちょっといろいろ変換が必要そう? ここもよくわからず雰囲気でやってるので変なところがあればご指摘ください。

m <- x$layer

あとはrayshaderに突っ込むだけです。

ambmat <- ambient_shade(m)

m %>%
  sphere_shade(texture = "imhof1") %>%
  # 今回あまり高低差がないとこだったので水が見えてなさそう
  add_water(detect_water(m), color = "imhof1") %>%
  add_shadow(ray_shade(m, zscale = 3, maxsearch = 300), 0.5) %>%
  add_shadow(ambmat, 0.5) %>%
  plot_3d(m, zscale = 50, fov = 0, theta = -45, phi = 45, windowsize = c(1000, 800), zoom = 0.75,
          water = TRUE, waterdepth = 0, wateralpha = 0.5, watercolor = "lightblue",
          waterlinecolor = "white", waterlinealpha = 0.5)
render_snapshot(clear = TRUE)

f:id:yutannihilation:20191020130134p:plain
出典:基盤地図情報DEMデータを加工して作成

以上、fgdrパッケージで読み込んだデータを雰囲気でrayshaderする例でした。