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)Rで富士山周辺の3D地形図を作成
R, RStudio, Quarto and STEM(Science, Technology, Engineering and Mathematics)
Rで富士山周辺の3D地形図を作成します。
rayshader パッケージを利用して、富士山周辺の実際の標高データ（DEM）から陰影付きの立体地形図を生成します。
1 基本の流れ
- 範囲指定：
sf::st_bbox()で対象範囲を作成し、st_as_sfc() %>% st_sf()で正式なsfオブジェクトに変換（sfcのままだとget_elev_raster()内部でエラーになるため必須）
- 標高データ取得：
elevatr::get_elev_raster()で富士山周辺のDEM（数値標高モデル）を取得
- 行列変換：
rayshader::raster_to_matrix()でrayshader用の行列に変換
- 陰影処理：
sphere_shade()/ray_shade()/ambient_shade()で陰影をつける
- 3Dレンダリング：
plot_3d()でrgl(OpenGL)ウィンドウに立体表示
- HTML埋め込み：
- チャンク先頭で
options(rgl.useNULL = TRUE)とrgl::setupKnitr(autoprint = TRUE)を設定し、レンダリングチャンクに#| webgl: trueを付けることで、インタラクティブな3DシーンをそのままHTMLページに自動埋め込み
- チャンク先頭で
3～5（陰影処理〜レンダリング）は%>%パイプで一気につなげられるので、実質的には「①範囲・sf整形 → ②DEM取得 → ③〜⑤陰影＆描画 → ⑥HTML埋め込み」という4ブロック構成になります。
2 初期設定
3 富士山周辺の3D地形
下記のRコードにより、標高データの取得からレンダリングまでを行い、チャンクオプション webgl: true によって結果がそのままページ内にインタラクティブなWebGLとして埋め込まれます（マウスドラッグで回転、スクロールでズーム可能）。
# -----------------------------------------------------
# 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--- sfオブジェクトの確認 ---
Simple feature collection with 1 feature and 0 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 138.5 ymin: 35.1 xmax: 139.2 ymax: 35.6
Geodetic CRS: WGS 84
geometry
1 POLYGON ((138.5 35.1, 139.2...# -----------------------------------------------------
# 2. 標高データ(DEM)を取得
# z: ズームレベル(解像度)。大きいほど精細になりますが
# データ量・処理時間が指数的に増えます
# 全国規模なら z = 5〜7 程度が現実的
# -----------------------------------------------------
dem <- get_elev_raster(bbox, z = 9, clip = "bbox")
elmat <- raster_to_matrix(dem)
cat('--- 数値標高モデルの確認 ---\n')
dem--- 数値標高モデルの確認 ---
class : RasterLayer
dimensions : 416, 575, 239200 (nrow, ncol, ncell)
resolution : 0.001252417, 0.001252417 (x, y)
extent : 138.4901, 139.2102, 35.08951, 35.61052 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
source : memory
names : file24e81bca3bac
values : -942, 3707 (min, max)# -----------------------------------------------------
# 3. rayshader用の行列へ変換
# -----------------------------------------------------
cat('--- 行列変換の結果を確認(一部) ---\n')
str(elmat)--- 行列変換の結果を確認(一部) ---
num [1:575, 1:416] 262 261 259 258 258 256 255 254 255 257 ...# -----------------------------------------------------
# 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))算出したzscale: 113.7# -----------------------------------------------------
# 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()を明示的に呼ぶ必要はない）