Rで国土地理院基盤地図情報データを扱う

この記事は、FOSS4Gアドベントカレンダー2018の16日目の記事です。タイトルにある通り、Rから国土地理院基盤地図情報ダウンロードサービスでダウンロードしたファイルを扱う方法について書きます。

基盤地図情報が提供するデータはXMLファイルとなっており、GISアプリケーション等で扱う際には変換処理が必要になります。データ変換ツールとして国土地理院のオフィシャルなものやオープンソースのPythonライブラリ、標高DEMデータに対応したものなど様々あり、車輪の再発明でもあるのですが、今回はR用のパッケージとしてまとめました。その名もfgdrパッケージです。fgdは基盤地図情報の英語表記 (Fundermental Geographic Data) からとりました。この記事ではfgdrパッケージの基本的な使い方を紹介します。

なおfgdrパッケージでは、基盤地図情報のデータをダウンロードする関数は用意していません。読み込むデータは各自でダウンロードしてください。また、読み込んだデータの扱いについても基盤地図情報ダウンロードサービスの利用規約に従ってください。（この記事中での国土地理院コンテンツの出典はページ末尾に記載しています。）

fgdrパッケージ

github.com

fdgrパッケージでは、基盤地図情報データダウロードサービスが提供している

基本項目
数値標高モデル (5m, 10mメッシュ)

について、データフレーム、Rの地理空間データを扱うためのクラスであるsfまたはraster (starsについては準備中)として返却する関数が備わっています。

数値標高モデルの変換方法は以前に書いたコードを整理したものとなります。5mだけ出なく10mメッシュに対応しました。

uribo.hatenablog.com

インストールはGitHubから行います。

remotes::install_github("uribo/fgdr")

library(sf)
library(raster)
library(fgdr)

基本項目

基盤地図情報の13項目のうち、「測量の基準点」、「海岸線」、「行政区画の境界線及び代表点」、「建築物の外周線」など10項目が基本項目として扱われています。基盤地図情報ダウンロードサービスでは、これらの項目をぞれぞれXMLファイルとして提供します。

基本項目のデータはread_fgd()のfile引数にファイルのパスを指定して読み込みます。この関数は対象のファイルの種類に応じて、自動的に地物の種類を判別して返却します。例えば、行政区画代表点 (AdmPt) では下記のようにポイントデータです。

sf_admpt <-
  read_fgd(file = "FG-GML-523346-AdmPt-20180701-0001.xml")

sf_admpt

## Simple feature collection with 4 features and 2 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 133.7835 ymin: 35.00697 xmax: 133.8736 ymax: 35.06013
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##              gml_id adm_name                  geometry
## 1 K6_4816212612_1-2   真庭市  POINT (133.7835 35.0416)
## 2 K6_4818012606_1-2   美咲町 POINT (133.8394 35.02376)
## 3 K6_4818312621_1-2   津山市 POINT (133.8496 35.06013)
## 4 K6_4819212600_1-2   津山市 POINT (133.8736 35.00697)

水域 (WL) はライン、行政区画 (AdmArea) はポリゴンデータになります。

sf_admarea <- 
  read_fgd("FG-GML-523346-AdmArea-20180701-0001.xml")

sf_admarea

## Simple feature collection with 4 features and 2 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 133.75 ymin: 35 xmax: 133.875 ymax: 35.08333
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##          gml_id adm_name                       geometry
## 1 K4_523346_1-2   美咲町 POLYGON ((133.7869 35, 133....
## 2 K4_523346_2-2   真庭市 POLYGON ((133.7879 35.00639...
## 3 K4_523346_3-2   津山市 POLYGON ((133.8742 35.00319...
## 4 K4_523346_4-2   津山市 POLYGON ((133.8746 35.03678...

読み込んだデータをggplot2で表示してみましょう。（といってもこの情報だけでは物足りませんが）

library(ggplot2)
ggplot() +
  geom_sf(data = sf_admarea) +
  geom_sf_label(data = sf_admpt, aes(label = adm_name))

f:id:u_ribo:20181216082822p:plain

ファイル内のデータ件数が多い場合や、いくつかの項目 (建築物の外周線、水涯線など)を読み込む際はちょっと時間がかかります。

# 件数が多い場合はリストになります...
sf_blda <- 
  read_fgd("FG-GML-523346-BldA-20180701-0001.xml")

mapviewパッケージを使って確認してみましょう。実際の建物の位置を参照したいため、ベースタイルにはEsriのWorldImageryを使います。

m <- 
  mapview::mapview(sf_blda[[1]], map.types = "Esri.WorldImagery")

xy <- 
  sf_blda[[1]] %>% 
  st_transform(crs = 6673) %>% 
  st_centroid() %>% 
  st_transform(crs = 4326) %>% 
  st_coordinates()

m@map <- 
  m@map %>% 
  leaflet::setView(lng = xy[1], lat = xy[2], zoom = 16)

m

f:id:u_ribo:20181216082555p:plain

数値標高モデル

標高のメッシュデータである数値標高モデル(DEM) のデータはデータフレームまたはrasterとして読み込みます。対象のファイルが保存されているパスおよび数値標高データの種類を指定したread_fgd_dem()実行します。

データがない地点や値が不明な場合、-9999に置き換わっている点に注意してください。

5mメッシュ

read_fgd_dem("FG-GML-5135-63-00-DEM5A-20161001.xml", 
             resolution = 5)

## # A tibble: 33,750 x 2
##    type       value
##    <chr>      <dbl>
##  1 データなし -9999
##  2 データなし -9999
##  3 データなし -9999
##  4 データなし -9999
##  5 データなし -9999
##  6 データなし -9999
##  7 データなし -9999
##  8 データなし -9999
##  9 データなし -9999
## 10 データなし -9999
## # ... with 33,740 more rows

デフォルトでは返り値のオブジェクトがデータフレームですが、これは引数return_class =によりrasterを選ぶことも可能です。rasterとして読み込み、可視化する例を示します。

r <- 
  read_fgd_dem("FG-GML-5135-63-00-DEM5A-20161001.xml", 
             resolution = 5,
             return_class = "raster")

r <- 
  r %>% 
  raster::setValues(r %>% 
  raster::getValues() %>% 
  dplyr::if_else(. == -9999, -1, .))

plot(r, col = cptcity::cpt("arendal_arctic", n = 10))
title(main = "数値標高モデル (5mメッシュ): 51356300",
      sub = "「基盤地図情報 数値標高モデル 5mメッシュ」（国土地理院）\n（https://fgd.gsi.go.jp/download/menu.php）\nをもとに瓜生真也(@uribo)が作成",
      cex.main = 1,
      cex.sub = 0.85,
      adj = 1)

f:id:u_ribo:20181216082806p:plain

10mメッシュ

同様に10mメッシュのデータを読み込むにはresolution =を10に変更してください。

read_fgd_dem("FG-GML-5440-10-dem10b-20161001.xml", 
             resolution = 10,
             return_class = "raster")

## class       : RasterLayer 
## dimensions  : 750, 1125, 843750  (nrow, ncol, ncell)
## resolution  : 0.0001111111, 0.00011112  (x, y)
## extent      : 140, 140.125, 36.08333, 36.16667  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : layer 
## values      : -9999, 316  (min, max)