【R】数値標高データの結合

こちらでは、国土地理院の基盤地図情報ダウンロードサービスでダウンロードできる「数値標高モデル」を表示してみました。ただ、データが1km四方の１つのデータに限られていました。
そこで今回は、もっと広い範囲を表示すべくデータの結合を行います。

今回も、こちらのサイトを参考にさせていただきました。

データは、同じく陣馬山周辺のデータで下記の通りです。

> files
 [1] "./FG-GML-5339-31-60-DEM5A-20161001.xml" "./FG-GML-5339-31-61-DEM5A-20161001.xml"
 [3] "./FG-GML-5339-31-62-DEM5A-20161001.xml" "./FG-GML-5339-31-63-DEM5A-20161001.xml"
 [5] "./FG-GML-5339-31-64-DEM5A-20161001.xml" "./FG-GML-5339-31-65-DEM5A-20161001.xml"
 [7] "./FG-GML-5339-31-66-DEM5A-20161001.xml" "./FG-GML-5339-31-67-DEM5A-20161001.xml"
 [9] "./FG-GML-5339-31-68-DEM5A-20161001.xml" "./FG-GML-5339-31-70-DEM5A-20161001.xml"
[11] "./FG-GML-5339-31-71-DEM5A-20161001.xml" "./FG-GML-5339-31-72-DEM5A-20161001.xml"
[13] "./FG-GML-5339-31-73-DEM5A-20161001.xml" "./FG-GML-5339-31-74-DEM5A-20161001.xml"
[15] "./FG-GML-5339-31-75-DEM5A-20161001.xml" "./FG-GML-5339-31-76-DEM5A-20161001.xml"
[17] "./FG-GML-5339-31-77-DEM5A-20161001.xml" "./FG-GML-5339-31-78-DEM5A-20161001.xml"
[19] "./FG-GML-5339-31-80-DEM5A-20161001.xml" "./FG-GML-5339-31-81-DEM5A-20161001.xml"
[21] "./FG-GML-5339-31-82-DEM5A-20161001.xml" "./FG-GML-5339-31-83-DEM5A-20161001.xml"
[23] "./FG-GML-5339-31-84-DEM5A-20161001.xml" "./FG-GML-5339-31-85-DEM5A-20161001.xml"
[25] "./FG-GML-5339-31-86-DEM5A-20161001.xml" "./FG-GML-5339-31-87-DEM5A-20161001.xml"
[27] "./FG-GML-5339-31-88-DEM5A-20161001.xml"

読み込むデータのリストを作ります。

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

# 結合したいxmlファイルをベクトルにします
files <- 
  list.files(pattern = ".+5339-31-[6-8][0-8].+.xml$",
                        full.names = TRUE)

データからラスターデータを作り出す処理を関数にします。

dem_xml2df <- function(file) {
  xml2::read_xml(file) %>% 
    xml2::xml_find_all("/*/*/*/gml:rangeSet/gml:DataBlock/gml:tupleList") %>% 
    xml2::xml_contents() %>% 
    as.character() %>% 
    readr::read_delim(delim = ",", col_names = c("type", "value"))
}
set_coords <- function(raster, meshcode){
  m <- 
    jpmesh::export_mesh(meshcode)
  
  bb <-
    m %>%
    sf::st_bbox() %>%
    as.numeric()
  
  raster::extent(raster) <- 
    raster::extent(bb[1], bb[3], bb[2], bb[4])
  raster::crs(raster) <- 
    sp::proj4string(as_Spatial(m))
  raster
}

結合データをつくります。

raster_dem_list <- 
  purrr::map2(
    files,
    stringr::str_remove_all(basename(files), "FG-GML-") %>% 
      stringr::str_remove_all("-DEM5A-20161001.xml") %>% 
      stringr::str_remove_all("-"),
    ~ dem_xml2df(.x) %>% 
      dplyr::pull(value) %>% 
      matrix(nrow = 225, ncol = 150) %>% 
      t() %>% 
      raster() %>% 
      set_coords(.y)
  )

raster_dem_merged <- 
  raster_dem_list %>% 
  purrr::reduce(merge)

プロットしてみます。

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

ちゃんとデータが結合されて、陣馬山だけではなく、その周辺の山塊も表示できています。OpenStreetMapでも見てみます。

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

もう少し、美しく表示してみます。

df_alt <-
  as.data.frame(raster_dem_merged, xy = TRUE) %>%
  tibble::as_tibble() %>%
  dplyr::rename("Elevation" = layer)

ggplot() +
  geom_raster(data = df_alt ,
              aes(x, y, fill = Elevation),
              hjust = 0, vjust = 0) +
  scale_fill_gradientn(colours = cptcity::cpt(pal = cptcity::cpt_names[[4]], n = 50)) +
  geom_contour(data = df_alt, aes(x, y, z = Elevation), col = "white", size = 0.2) +
  coord_quickmap() +
  theme_void(base_family = "IPAexGothic") +
  theme(legend.key.height = unit(2, "lines")) +
  labs(title = "陣馬山周辺",
       caption = "データ：国土地理院 基盤地図情報数値標高モデル")