【R】jpmesh

1. はじめに

国交省が提供する基盤地図情報など、標準地域メッシュを扱いたい時があります。僕は、基盤地図情報ダウンロードサービスの地図から目星をつけてデータをダウンロードしてくることが多かったのですが、せっかくダウンロードしても本当に欲しい位置のデータではないこともしばしばありました。このjpmeshパッケージは、得たい地図の座標からメッシュを得ることができてとっても便利です。逆に、メッシュから座標の取得もできるという素晴らしい機能もあります。使ってみましょう。作者が回折してくれているページ「標準地域メッシュを扱うRパッケージを更新しました: jpmesh v.1.0.0」を参考にさせていただきました。

2. インストール

CRANに登録されているので簡単にインストールできます。

install.packages("jpmesh")

3. 使ってみる

では、使ってみます。赤城山のメッシュをとってみます。

library(jpmesh)
library(magrittr)

coords_to_mesh(longitude = 139.168197, latitude = 36.539656) #赤城山

> coords_to_mesh(longitude = 139.168197, latitude = 36.539656)
[1] "54396143"

おお！簡単！素晴らしい！

逆に、メッシュから座標も求められます。

mesh_to_coords("54396143")

> mesh_to_coords("54396143")
# A tibble: 1 x 4
  lng_center lat_center lng_error lat_error
       <dbl>      <dbl>     <dbl>     <dbl>
1       139.       36.5   0.00625   0.00417

すごい！

さらに、隣接したメッシュも求められます。広い範囲の地図を使うときには有効ですね。

coords_to_mesh(longitude = 139.168197, latitude = 36.539656) %>%
  fine_separate()

> coords_to_mesh(longitude = 139.168197, latitude = 36.539656) %>%
+   fine_separate()
[1] "543961431" "543961432" "543961433" "543961434"

coords_to_mesh(longitude = 139.168197, latitude = 36.539656) %>%
  find_neighbor_mesh()

> coords_to_mesh(longitude = 139.168197, latitude = 36.539656) %>%
+   find_neighbor_mesh()
[1] "54396132" "54396133" "54396134" "54396142" "54396143" "54396144" "54396152"
[8] "54396153" "54396154"

では、実際に赤城山を描いてみましょう！

上記の場所から、ちょっとだけ位置をずらして赤城山総合観光案内所周辺にしてみます。具体的には、次の通りです。

> coords_to_mesh(longitude = 139.168044, latitude = 36.543367) 
[1] "54396153"

とりあえず、描いてみます。

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

# 5mメッシュを格納したファイル
tgt_file <- "FG-GML-5439-61-53-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)

xy_lc


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

head(df_dem, 3)


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メッシュ): 5439-61-53")

次に周辺の地形も含めます。範囲は、下記の通りです。

> coords_to_mesh(longitude = 139.168044, latitude = 36.543367) %>%
+   find_neighbor_mesh()
[1] "54396142" "54396143" "54396144" "54396152" "54396153" "54396154" "54396162"
[8] "54396163" "54396164"

これらのデータをマージして図形を描画します。

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

# 結合したいxmlファイルをベクトルにします
files <- 
  list.files(pattern = ".+5439-61-[4-6][2-4].+.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メッシュ): 543961",
      sub = "データ：国土地理院 基盤地図情報数値標高モデル")

さらに、rayshaderで。

# 3D model
library(rayshader)
elmat = raster_to_matrix(raster_dem)

#We use another one of rayshader's built-in textures:
elmat %>%
  sphere_shade(texture = "imhof1") %>%
  plot_map()

elmat %>%
  sphere_shade(texture = "imhof1") %>%
  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(clear=TRUE)

4. さいごに

メッシュ情報を簡単に取得できる素晴らしいパッケージです。ありがとうございます！