(私も)Rで基盤地図情報を扱いたい
はい,どーも僕です.
戸田恵梨香の朝ドラが楽しみです.
はじめに
先日,R言語のライブラリーfgdrが githubにて公開されました. fgdrでは国土地理院が公開している基盤地図情報のXMLファイルを 読み込むための関数を提供しています.
基盤地図情報の変換ライブラリーがあったら便利だなと以前から思っていたため, ライブラリーの公開はありがたいです. さっそくfgdr使ってみたところ,いくつか修正したい点ありました. ここではそれらを指摘するとともに, 修正点に対応したスクリプト案を示します.
read_fgd
修正したい点
DEM以外のデータを読み込みsfで返す
read_fgd
関数に対して,次に示す2点を
修正したと感じました.
このとき検証に使用したデータはニ次メッシュ563804です.
- データサイズが10000を超える場合にリストで返す
- データが欠測しているときに読み込みエラーとなる
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()
> 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
で読み込んだ
方が良いのかと思います.