RからWFSを利用する

R で地理データを扱う時に、WFS から直接取り込むことができたら便利に違いない。 幾つか方法があるのかも知れないが、google でまず GDAL を使う方法が ヒットしたので試してみる。

gdalUtils と rgdal というパッケージを使うのだが、実態は、GDAL のコマンドを 呼び出すラッパーのようなので、先に直接コマンドを実行してみる。

さて、まず判り難かったのが、GDAL で WFS ドライバをビルドするには configure で curl オプションを有効にしなければならないということ。 OSX の macports、および FreeBSD の pkg で素直にインストールすると curl は有効になっていないので WFS が使えない。 macports では次のようにする。

$ sudo port install gdal +curl

FreeBSD では ports でリビルドしなければならないので少々面倒臭い。

curl オプション付きで GDAL をインストールすると、ogrinfo で WFS のメタデータを取得できるようになる。以下で利用可能なレイヤの一覧が出力される。

$ ogrinfo "WFS:http://hogehoge/wfs?service=wfs&version=1.0.0&request=GetCapabilities&srsname=EPSG:4326"

ogr2ogr で WFS から取得した GIS データをローカルファイルに保存できる。 TypeName に欲しいレイヤ名を、 request は GetFeature にする。 なお、元データにマルチバイト文字があるとエンコードエラーが出るのと、 線とポリゴンが含まれているとジオメトリ変換エラーが出る。 どちらも、デフォルトの出力フォーマットが shp であることが原因のようだ。 前者は -lco オプションでエンコードを指定して回避する。 後者は、shp がそもそも複数のジオメトリ型をひとつのファイルに 共存させられないから回避できない(処理するジオメトリ型を明示的に指定して 複数回実行するしかない)。出力フォーマットはオプションで 変更可能だが、まずはエラーを無視するように -skipfailures を指定してみる。

$ ogr2ogr test.shp "WFS:http://hogehoge/wfs?service=wfs&version=1.0.0&request=GetFeature&typename=mylayer&srsname=EPSG:4326" -lco "encoding=utf-8" -skipfailures

出力フォーマットは -f オプション で指定する。GML だとエンコーディング問題も 複数型共存問題も生じない。

$ ogr2ogr test.gml "WFS:http://hogehoge/wfs?service=wfs&version=1.0.0&request=GetFeature&typename=PlotNet:plots&srsname=EPSG:4326" -f gml