【R】数値標高データをRayshaderで表示

2020年7月12日

タイトルの通り、国土地理院の基盤地図情報ダウンロードサービスでダウンロードできる「数値標高モデル」をrayshaderで表示してみました。またも、陣馬山のデータを使用します。

前半の数値標高モデルからラスターデータを作成する箇所は、こちらと同じです。一応、コードを書いてきます。

library(tidyverse)
library(xml2)
library(jpmesh)
library(sf)
library(raster)

tgt_file <- "FG-GML-5339-31-83-DEM5A-20161001.xml"

x <- read_xml(tgt_file)

xy_lc <- 
  x %>% 
  xml_find_all("/*/*/*/gml:boundedBy/gml:Envelope/gml:lowerCorner") %>% 
  xml_contents() %>% 
  as.character() %>% 
  strsplit(" ") %>% 
  unlist() %>% 
  map_dbl(as.numeric)

mesh <- 
  jpmesh::coords_to_mesh(xy_lc[2], xy_lc[1])
mesh

x %>% 
  xml_find_first("/*/*/*/gml:gridDomain/gml:Grid/gml:limits/gml:GridEnvelope")


df_dem <- 
  x %>% 
  xml_find_first("/*/*/*/gml:rangeSet/gml:DataBlock/gml:tupleList") %>% 
  xml_contents() %>% 
  as.character() %>% 
  readr::read_delim(delim = ",", col_names = c("type", "value"))

raster_dem <- 
  df_dem$value %>% 
  matrix(nrow = 225, ncol = 150) %>% 
  t() %>% 
  raster()

mesh <- 
  jpmesh::export_mesh(mesh)

bb <-
  mesh %>%
  st_bbox() %>%
  as.numeric()

extent(raster_dem) <- extent(bb[1], bb[3], bb[2], bb[4])

crs(raster_dem) <- sp::proj4string(as_Spatial(mesh))

そして、rayshaderにて表示してみます。

library(rayshader)
elmat = raster_to_matrix(raster_dem)

まずは、平面的に。

elmat %>%
  sphere_shade(texture = "desert") %>%
  plot_map()

次は立体的に表示してみます。

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))
Sys.sleep(0.2)
render_snapshot()