【R】太白山

2020年8月4日

仙台市太白区に太白山という山があります。（僕の記憶では）結構きれいな独立峰です。標高は低い（保育園の遠足で行くらしい）けども、遠くから見るとなかなか感じのいい山です。

国土地理院の基盤地図情報ダウンロードサービスでこの太白山の「数値標高モデル」をデータとしてダウンロードして表示してみます。

こちらのサイトを参考にRにて表示してみます。

太白山のデータのメッシュコードは「574026」（５ｍメッシュ）です。データは「FG-GML-5740-26-84-DEM5A-20161001.xml」でした。これを見つけるのに時間がかかった・・・。

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

tgt_file <- "./FG-GML-5740-26-DEM5A/FG-GML-5740-26-84-DEM5A-20161001.xml"

x <- read_xml(tgt_file)

* lower-left cornerを取得
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])

地表面の高さ情報取得。

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()

raster_dem

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))

プロットしてみます。

plot(raster_dem)
plot(rasterToContour(raster_dem), add = TRUE)
title(main = "太白山",
      sub = "データ：国土地理院 基盤地図情報数値標高モデル　(5mメッシュ): 5740-26-84")

左下に太白山が！きれいな同心円の等高線です（笑）。学生時代には、ここの地図も頑張って作ったなあ。。。。

せっかくなので、rayshaderを使って、3次元表示してみます。

library(rayshader)
elmat = raster_to_matrix(raster_dem)

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()