(私も)Rで基盤地図情報を扱いたい

はい,どーも僕です.

戸田恵梨香の朝ドラが楽しみです.

はじめに

先日,R言語のライブラリーfgdrgithubにて公開されました. fgdrでは国土地理院が公開している基盤地図情報XMLファイルを 読み込むための関数を提供しています.

uribo.hatenablog.com

基盤地図情報の変換ライブラリーがあったら便利だなと以前から思っていたため, ライブラリーの公開はありがたいです. さっそくfgdr使ってみたところ,いくつか修正したい点ありました. ここではそれらを指摘するとともに, 修正点に対応したスクリプト案を示します.

read_fgd

修正したい点

DEM以外のデータを読み込みsfで返す read_fgd関数に対して,次に示す2点を 修正したと感じました. このとき検証に使用したデータはニ次メッシュ563804です.

  1. データサイズが10000を超える場合にリストで返す
  2. データが欠測しているときに読み込みエラーとなる

1点目はファイルに含まれている要素数が10,000より多い場合と 少ない場合とで異なる構造で返されるというものです. 作者のブログでも言及されています. クラスに依りますが基盤地図情報では10,000要素を超える 要素数のファイルは少なくないため,修正の必要があると思います.

2点目はデータが仕様と異なる場合に生じるものです. たとえばFG-GML-563804-RdCompt-20180101-0001.xmlでは データのIDにあたる<RdCompt>ノードが, データの属性となる<vis>ノードを持っていない場合があります. その場合<vis>ノードをテーブルの属性にしようとすると ID数とデータサイズが異なりエラーとなります. 当該状況でどのような挙動に実装するのかは作者に依りますが, 私としては欠測する場合でもNAを値としてエラーを避けて欲しいと 考えています.

修正案

本来はプルリクエストを出すべきとは思いますが, もとの実装のどこを直すべきかがわかりませんでした. そのため自分で書いたものを示します.

もともとの実装で10,000要素を境に処理を分けている理由はわかっていませんが, ここで示すスクリプトでは50,000要素程度でもsfのテーブルとして 返せることを確認しています. データが欠測する場合の対応としてxml_find_allではなく, xml_find_firstを利用しています.

こんな感じで使えます.sfで返ってきます.

rdcompt <- readfgd("FG-GML-563804-ALL-20180101/FG-GML-563804-RdCompt-20180101-0001.xml", 4612)
select(rdcompt, 1) %>% plot()

f:id:aaaazzzz036:20181221214552p:plain

> rdcompt
Simple feature collection with 55169 features and 8 fields
geometry type:  LINESTRING
dimension:      XY
bbox:           xmin: 138.5 ymin: 37.33333 xmax: 138.625 ymax: 37.41667
epsg (SRID):    4612
proj4string:    +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
First 10 features:
              gmlid                fid   lfSpanFr    devDate orgGILvl  vis   type admOffice                       geometry
1  K24_4986013443_1  49860-13443-s-571 2014-08-11 2017-12-27     2500 表示 分離帯      不明 LINESTRING (138.5083 37.347...
2  K24_4986013443_2  49860-13443-s-572 2014-08-11 2017-12-27     2500 表示 分離帯      不明 LINESTRING (138.5015 37.348...
3  K24_4986013443_3  49860-13443-s-573 2014-08-11 2017-12-27     2500 表示 分離帯      不明 LINESTRING (138.5015 37.348...
4  K24_4986013443_4  49860-13443-s-574 2014-08-11 2017-12-27     2500 表示 分離帯      不明 LINESTRING (138.5015 37.348...
5  K24_4986013443_5  49860-13443-s-953 2016-05-11 2017-12-27     2500 表示 分離帯      不明 LINESTRING (138.5 37.34831,...
6  K24_4986013443_6  49860-13443-s-954 2016-05-11 2017-12-27     2500 表示 分離帯      不明 LINESTRING (138.5 37.34835,...
7  K24_4986013443_7  49860-13443-s-955 2016-05-11 2017-12-27     2500 表示 分離帯      不明 LINESTRING (138.5 37.34831,...
8  K24_4986013443_8  49860-13443-s-956 2016-05-11 2017-12-27     2500 表示 分離帯      不明 LINESTRING (138.5 37.34835,...
9  K24_4986013446_1 49860-13446-s-1124 2014-08-11 2017-12-27     1000 表示   側溝      不明 LINESTRING (138.5083 37.353...
10 K24_4986013446_2 49860-13446-s-1125 2014-08-11 2017-12-27     1000 表示   側溝      不明 LINESTRING (138.5081 37.353...

不安というか,TODOというか穴あきポリゴンがある場合には対応していません.

おまけ

read_fgd_dem

read_fgd_demには不具合がありましたが プルリクエスト(人生初の)を出したところ マージして頂けたため,私が気づいている不具合は解消されています.

蛇足ですが,自分用に作成していたDEMのラスター化のスクリプトを示します. read_fgd_demとの違いはCRSを自分で設定できることです. 現在はもうないと思いますが作成当時のGISではJGD2011に 対応していないものもあったので,あえてJGD2000を設定するという場合に 役だってました. まっ,こんな実装もありますよ程度です.

OGR

上で示したreadfgdは正直遅いです.手元のノートパソコンでFG-GML-563804-BldA-20180101-0001.xmlを 読み込んだところ約40秒必要でした. ギリギリ許容範囲ではあるけどな,という感じです. そのため,ogr2ogrコマンドでxmlをシェープに変換してしまってから, sf::st_readで読み込んだ 方が良いのかと思います.