【R】数値標高データの表示

国土地理院の基盤地図情報ダウンロードサービスで「数値標高モデル」をデータとしてダウンロードできます。僕は、トレランや山登りを楽しんでいるので標高付きの地図データというのが大好きです。これをRを使って表現したいと思っていました。

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

陣馬山の地図を表示させたいと思います。メッシュコードは、533931です。データは、FG-GML-5339-31-83-DEM5A-20161001.xmlで、５ｍのメッシュのデータです。
まず、データを読み込みます。

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

地表面のデータを取り出します。

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

プロットしてみます。

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

地図の表示ができました。陣馬山っぽいですね。

OpenStreetMapに重ねてみます。

mapview::mapview(raster_dem, alpha.regions = 0.2)