Rで富士山周辺の3D地形図を作成

library(rayshader)
library(elevatr)
library(sf)
library(terra)
library(magrittr)
library(rgl)

# 重要: quarto render をCUI/サーバー等のヘッドレス環境で実行する場合、
# 物理ディスプレイがないとrglが初期化に失敗することがあるため
# 必ずNULLデバイスを使うよう指定しておく
options(rgl.useNULL = TRUE)

# knitrにrglシーンを自動的にHTML(WebGL)へ変換させるためのセットアップ
# これにより、各チャンクの最後に残っているrglシーンが
# チャンクオプション webgl: true 付きで自動的に埋め込まれる
rgl::setupKnitr(autoprint = TRUE)

# -----------------------------------------------------
# 1. 対象範囲の境界を取得
# -----------------------------------------------------

# 注意: st_as_sfc()のみですとsfc(ジオメトリのみ)になり、
#       get_elev_raster()内部のrbind処理でエラーになるため
#       st_sf()で正式なsfオブジェクトに変換する
bbox <- st_bbox(c(
  xmin = 138.5, xmax = 139.2,
  ymin = 35.1,  ymax = 35.6
), crs = 4326) %>%
  st_as_sfc() %>%
  st_sf()

cat('--- sfオブジェクトの確認 ---\n')
bbox

# -----------------------------------------------------
# 2. 標高データ(DEM)を取得
#    z: ズームレベル(解像度)。大きいほど精細になりますが
#       データ量・処理時間が指数的に増えます
#       全国規模なら z = 5〜7 程度が現実的
# -----------------------------------------------------

dem <- get_elev_raster(bbox, z = 9, clip = "bbox")
elmat <- raster_to_matrix(dem)

cat('--- 数値標高モデルの確認 ---\n')
dem

# -----------------------------------------------------
# 3. rayshader用の行列へ変換
# -----------------------------------------------------

cat('--- 行列変換の結果を確認(一部) ---\n')
str(elmat)

# -----------------------------------------------------
# 3.5 zscaleを実際のラスタ解像度から動的に計算
#     固定値(例: 30)を使うと、elevatrのzoom levelによって
#     実際の解像度(緯度依存)と乖離し、地形が不自然に
#     誇張/平坦化されるため、必ずデータから算出する
#       zscale = 標高(m)に対する、格子間の実際の水平距離(m)
# -----------------------------------------------------

dem_t <- rast(dem)
res_deg <- res(dem_t)[1]
lat_center <- mean(c(ymin(dem_t), ymax(dem_t)))
zscale <- res_deg * 111320 * cos(lat_center * pi / 180)
cat(sprintf("算出したzscale: %.1f\n", zscale))

# -----------------------------------------------------
# 4. 陰影をつけて3Dレンダリング
#    plot_3d() はrgl(OpenGL)のインタラクティブウィンドウを開く
# -----------------------------------------------------

elmat %>%
  sphere_shade(texture = "imhof1") %>%
  add_shadow(ray_shade(elmat, zscale = zscale), 0.5) %>%
  add_shadow(ambient_shade(elmat), 0) %>%
  plot_3d(
    elmat,
    zscale = zscale,
    fov = 0,
    theta = 135,
    zoom = 0.75,
    phi = 45,
    windowsize = c(1000, 800)
  )

# plot_3d()実行後、チャンクの終わりでrglシーンが残っていれば
# setupKnitr()により自動的にwebglとして埋め込まれる
# （rglwidget()を明示的に呼ぶ必要はない）