国土地理院の数値標高モデルデータをラスタとしてRで扱う
国土地理院が提供している「基盤地図情報ダウンロードサービス」の中に「数値標高モデル」データがあります。このデータは、標高のメッシュデータです。標高点格子（メッシュサイズ）が5m、10mのものがそれぞれあります。利用登録をすればデータをダウンロードできます。
大変有益なデータなのですが、データを閲覧するには専用のビューアアプリが必要になります。しかしこれはWindows専用です。そのため、これまで私は、下記に挙げるような第三者が開発するデータ変換用のツール、QGISのプラグイン等を利用してデータを変換していました。
一方でこれらのツールはそれぞれ独立しており、Rと組み合わせて使うことはできません。しかし「数値標高モデル」自体はXMLファイルですので、標高値を紐解けばRで利用可能な形式へ変換することが可能です。今日はこの方法を紹介します。すでにPythonを使って同様のことをしている方がいらっしゃり、それのR版ということになります。
以下、数値標高データ本体とデータをRで扱うためのコードを説明をします。また簡単な可視化もしてみます。
「数値標高モデル」データの抽出
提供されるデータにはいくつかのバージョンがありますがそれに応じて、XMLスキーマ定義およびファイルの仕様書が用意されています。その中から重要なものを説明します。なおここで扱うファイルの仕様は、最新の基盤地図情報ダウンロードデータファイル仕様書4.1
に準じます。
なにはともあれ、まずはダウンロードしたXMLファイルを読み込みましょう。ファイルの読み込みにはxml2パッケージを利用します。
国土地理院のページからダウンロードした圧縮ファイル(一つのファイルでカバーする領域は10kmメッシュ)を解凍すると、10mメッシュのファイルでは1つのファイル、5mメッシュでは領域に含まれる複数のファイルが展開されます。分割されていない分、10mメッシュを扱う場合にはデータサイズが大きくなってしまうので注意が必要です。
library(purrr) library(xml2) # 5mメッシュを格納したファイル tgt_file <- "FG-GML-5440-34-00-DEM5A-20161001.xml"
XMLファイルはxml2::read_xml()でファイルを読み込みます。なお10mメッシュを対象とする際は引数optionsで"HUGE"を指定してください。
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
## [1] 36.25 140.50
# 座標からメッシュコードを特定 (ダウンロードしたファイルのメッシュコードと一致することを確認) mesh <- jpmesh::coords_to_mesh(xy_lc[2], xy_lc[1], mesh_size = "1km") mesh
## [1] "54403400"
メッシュに含まれるグリッドセルの数は以下のように格納されています。0が起点となっているので、5mメッシュの場合には225*150の標高値が格納されていることを示します。この値は10mメッシュでは1124,749となります(1125*750)。
x %>% xml_find_first("/*/*/*/gml:gridDomain/gml:Grid/gml:limits/gml:GridEnvelope")
## {xml_node}
## <GridEnvelope>
## [1] <gml:low>0 0</gml:low>
## [2] <gml:high>224 149</gml:high>
標高値の配列順序と開始点は次のコードで確認できます。ここで示した例は基盤地図情報データでの典型的なもので、先頭セルを北西端におき、配列順序がx軸の正方向へ進み、東端に達するとy軸の負方向に進み、南東端に至ることを示しています。くだけた表現でいうと、グリッドの西から東、北から南へとジグザグに進む方式です。
# +x-y が標準の配列順序 x %>% xml_find_all("/*/*/*/gml:coverageFunction/gml:GridFunction/gml:sequenceRule")
## {xml_nodeset (1)}
## [1] <gml:sequenceRule order="+x-y">Linear</gml:sequenceRule>
x %>% xml_find_all("/*/*/*/gml:coverageFunction/gml:GridFunction/gml:startPoint")
## {xml_nodeset (1)}
## [1] <gml:startPoint>0 0</gml:startPoint>
肝心の標高値は次のように抽出します。ここではデータフレームとして後の操作を扱いやすくしています。type列にDEM構成点の種別、value列に標高値を格納するようにしました。
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)
## # A tibble: 3 x 2 ## type value ## <chr> <dbl> ## 1 地表面 24.3 ## 2 地表面 23.9 ## 3 地表面 23.6
ちなみにXMLファイルの構造を把握するには、CRANには未登録のパッケージですが
xmltools を使うと便利でした。
ラスタへの変換
次は読み込んだ値を元にしてラスタデータへ変換するという処理です。また、地理空間データとして扱うための処理も行なっていきます。ラスタへの変換は、標高値を記録可能なサイズの行列オブジェクトを作り、raster::raster()を実行するだけです。これによりrasterに対応した可視化や地形解析のための関数が適用できるようになります。
library(sf) library(raster)
raster_dem <- df_dem$value %>% matrix(nrow = 225, ncol = 150) %>% t() %>% raster() raster_dem
## class : RasterLayer ## dimensions : 150, 225, 33750 (nrow, ncol, ncell) ## resolution : 0.00444444444444, 0.00666666666667 (x, y) ## extent : 0, 1, 0, 1 (xmin, xmax, ymin, ymax) ## coord. ref. : NA ## data source : in memory ## names : layer ## values : 12.81, 31.98 (min, max)
一方でこの状態では参照座標系が与えられていないので、メッシュデータの座標を元に定義します。
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メッシュ): 54403400", sub = "国土地理院 基盤地図情報数値標高モデルのデータを元に作成")
mapview::mapview(raster_dem, alpha.regions = 0.2)
OpenStreetMapのタイルに重ねてみます。
library(tmap) slope <- terrain(raster_dem, opt = "slope", unit = "radians", neighbors = 8) aspect <- terrain(raster_dem, opt = "aspect", unit = "radians", neighbors = 8) shade <- hillShade(slope, aspect) tm_shape(shade) + tm_raster(style = "cont", palette = "Greys") + tm_legend(show = FALSE)
これらの処理を関数化しておきました。
複数のファイルを結合する
さて、5mメッシュでは対象領域のメッシュが複数のファイルに分割されており、1ファイルの領域は1kmメッシュの範囲となっています。そこで複数のラスタを結合し、範囲を拡大することを試みます。XMLファイルの読み込みからラスタへの変換部分は先ほどの処理を関数化して使いまわします。
# 結合したいxmlファイルをベクトルにします files <- list.files(pattern = ".+5440-20-[6-8][7-9].+.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)
今度はggplot2を使って可視化を行います。標高のグラデーションを表現するために
http://soliton.vm.bytemark.co.uk/pub/cpt-city/index.html
のカラーパレットをRで利用可能にするcptcityパッケージを使いました。
library(ggplot2) 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 = "国土地理院 基盤地図情報数値標高モデルのデータを元に作成")
筑波山周辺の標高データを利用しました。綺麗な二つの頂が表現できていますね！
Enjoy!
次回はこのデータを使って、rayshaderパッケージによる下の図のような立体表現に挑戦します。
