rayshaderでDEM+GPSデータ可視化 #1
GPSデータを3次元で可視化したい
以下のTwitterをみて、やりたくなった。
‼️New #rayshader release (v0.32.5)! This version brings 3rd person views to your maps, generated directly from a spatial path: turn your GPS bike rides/runs/hikes into cool 3D animations, directly in #RStats! You can even steady the camera or add rotation for a cinematic effect. pic.twitter.com/rYE2dG6n75
— Tyler Morgan-Wall @tylermorganwall@fosstodon.org (@tylermorganwall) December 13, 2022
Population within 20 km of Japan's Tōkaidō Shinkansen, the world's first high-speed rail line. I swear every time I share something, I get new ideas I wish I'd implemented.
— terence fosstodon (@researchremora) December 18, 2022
But I need to fix the spelling error in my Germany map first.#rayshader adventures, an #rstats tale pic.twitter.com/WWJUthr8SF
Rのrayshaderというパッケージで、3次元で可視化する、という物だ。
よくクロスカントリースキーやランニングをするため、その軌跡を3Dで視覚化できればな、と思っていたため、試してみることにした。Rにはすこしなじみがあったため、ちょっと触って、Twitterで自慢する、くらいの気持ちだった。
が、割とスタックした。その記録。
最初
以下の記事を見て、rayshaderanimateというパッケージ入れてコードパクればいけるやん、と思っていた。
ただこのパッケージは、元となるDTMをどこかのサイトからダウンロードしてくる。自分は日本の地理院のデータ、ゆくゆくはUAV-LiDARで撮影したデータを可視化したかった。
そのため、ビデオを作るという目的で上のパッケージ使用はあきらめ(後程GPSデータ変換では使う)、本家rayshaderを使うことにした。インストールはgithubにある最新版から行った方が良い。
構築環境
R 4.2.2
Rstudio最新版(2023年1月16日時点)
rayshader 0.34.2
rayrender 0.29.0
terra1.6-47
DTMの可視化
地理院の基盤地図情報をQGISに読み込み、描画したいエリアのみ切り取ってGeoTiffに保存した。
(これ以外の方法には、Rのみでやる方法(↓)もあるっぽい。こちらの方法に従うと、これから説明する問題は発生しない可能性ありだが未確認。)
DTM自体はすぐに読みこめ、あとはrayshaderHPのチュートリアルに従うとすぐに3次元にできた。
library(raster)
library(rayshader)
raster_path <- "path.tif" #
r <- raster(raster_path)
elmat <- rayshader::raster_to_matrix(r)
#あとはチュートリアルのとおり
![](https://assets.st-note.com/img/1673773270734-Hr7Dd7WqWP.png?width=1200)
GPSデータの準備
手元にあったデータを使うことにした。なお、ここで上記使わないことにした、といったrayshaderanimateを使う。
library(rayshaderanimate)
gpx_path <- "path.gpx"
#緯度経度高度時間+追加情報のリストを作る。自分の場合、追加情報には心拍数データが入った。
gpx_table <- rayshaderanimate::get_table_from_gpx(gpx_path)
#緯度経度高度に分けてリスト作成。もっと良い方法がある気が?
{
ski_track_lat = list()
ski_track_long = list()
ski_track_alt = list()
nrow <- NROW(gpx_table)
for(i in 1:nrow) {
ski_track_lat[[i]] = gpx_table$lat[i]
ski_track_long[[i]] = gpx_table$lon[i]
ski_track_alt[[i]] = as.numeric(gpx_table$ele[i])
}
}
ぶち当たった問題
チュートリアルには、3次元で点およびラインを可視化する方法がある。ただ、いきなり「montereybay」というデータを使い始めるのだが、これの素性が分かりにくい。
要は、上記「elmat」、つまり
elmat <- rayshader::raster_to_matrix(r)
でrayshaderで扱うデータのスタイル?に変換したものだ。
これを元に、パスを以下のコードで描画しようと思った。
elmat %>%
sphere_shade(texture = "desert") %>%
#add_water(detect_water(elmat), color = "desert") %>%
add_shadow(ray_shade(elmat, zscale = 3), 0.5) %>%
add_shadow(ambient_shade(elmat), 0) %>%
plot_3d(elmat, zscale = 10, fov = 0, theta = 135, zoom = 0.75, phi = 45, windowsize = c(1000, 800))
render_path(extent = attr(elmat,"extent"),
lat = unlist(ski_track_lat), long = unlist(ski_track_long),
altitude = unlist(ski_track_alt), zscale=10,color="white", antialias=TRUE)
すると無限に
Error in get_extent(extent) :
# class of extent (`NULL`) not one of supported types (`Extent`, `bbox`, `numeric`, `SpatExtent`)
というエラーが出てわけわからんくなった。
解決策
結論から言うと、montereybayには"extent"という属性が付されており、これは要はラスタの範囲を示すものだ(と思われる)。
自分が読み込んだelmatには属性がなく、そこでエラーが発生していた。
attr(elmat, "extent") <- terra::ext(r, cells=NULL)
以上のコードを実行、その結果きちんと描画出来た。
![](https://assets.st-note.com/img/1673775180793-9FocUTIbNV.png?width=1200)
コードの意味を理解せず実行したので、エラーが出るのも当然だっ
た。
以上!
この記事が気に入ったらサポートをしてみませんか?