RでGISをやるときにはsfパッケージ、という世の中になるらしい。

RでGISといえばspパッケージです。でした。いま、時代が動こうとしています。

...という記事を書けるほどのGISの知識が私にはないので、ほんとはもっと勉強してから書くべきなんですけど、とりあえず勉強のためにも調べたことをまとめとこう、と思って書きます。怪しいところがあればツッコミをください。

spパッケージとは

spは、data.frameを拡張した地理情報データのためのデータ形式と、それを操作する基礎的な関数群を提供するパッケージです。以下のページの「Reverse depends」や「Reverse imports」を見ればわかるように、数百のパッケージがspに依存しています。

CRAN - Package sp

そんなspパッケージの歴史については、以下のブログ記事の前半で触れられています。見てみましょう。

Simple Features Now on CRAN | R Consortium

In 2003, a group of package developers sat together and decided to adopt a shared understanding of how spatial data should be organized in R. This led to the development of the package sp and its helper packages rgdal and rgeos.

spパッケージがつくられたのは、2003年。実に14年前のことです。rgdalパッケージとrgeosパッケージも同時期に作られたそうです。

After 2003, the rest of the world has broadly settled on adopting a standard for so-called “features”, which can be thought of as “things” in the real world that have a geometry along with other properties.

しかし、そのあと、R以外の世の中では「feature」をベースにしたデータ形式が広く採用され、OGCのsimple feature accessという標準にまとめられるに至ります。

(略) has been adopted widely over the past ten years, not only by spatial databases such as PostGIS, but also more recent standards such as GeoJSON.

このsimple feature accessはたとえば、PostGISやGeoJSONにも採用されています。

しかし、spパッケージは、というよりRのGISまわりの基礎は、この標準以前に作られたもので、simple featureのデータをうまく扱うことができていません。

具体的には、

  • Rに読み込めるが、異なる型に変換されてしまうデータ型がある
  • Rに読み込むことができないデータ型がある
  • 三次元データやM座標を含む二次元データをまったく扱えない

といった問題があるらしいです(Simple features for Rの最後に指摘されている)。

Simple Features for Rプロジェクト

この状況を打開しようと立ち上がったのが、spの作者でもあるEdzer Pebesma氏です。Edzerは、「Simple Features for R」というプロジェクトを立ち上げ、R Consortiumのグラントを獲得します。

Simple Features Access For R
Proposed By: Edzer Pebesma (Institute for Geoinformatics, University of Muenster)
Funded: 10,000USD
Project URL: TBA
Project Summary: Using the “Simple Features” standard supported by the Open Geospatial Consortium and the International Organization for Standardization, this tool will simplify analysis on modern geospatial data.
(Awarded Projects | R Consortium)

その成果として出来上がったのが、sfパッケージです。

CRAN - Package sf

ちなみに、sfパッケージのための試行錯誤は以下のブログにいくつか綴られていて、変遷を追っていくと興味深い気がします(私は追いきれませんでした)。

r-spatial

sfパッケージとは

sfパッケージの説明は、以下のvignetteに詳しいです。simple featuresとはなにか、というところから始まって、sfパッケージの使い方まで丁寧に説明されています。

冒頭で、sfの位置づけについての記述があります。

R has well-supported classes for storing spatial data (sp) and interfacing to the above mentioned environments (rgdal, rgeos), but has so far lacked a complete implementation of simple features, making conversions at times convoluted, inefficient or incomplete. The package sf tries to fill this gap, and aims at succeeding sp in the long term.

ということで、sp + rgdal + rgeosを置き換えていくものであることが控えめに宣言されています。「long term」というのがどれくらいかわかりませんが、すでにR界はsfのサポートに向けて動き出しているように見えます。

上のvignetteでは、sfのデータ構造について詳しく書かれてます。残りのvignetteもかなり本気の文章量です。どういうことかというと、つまり......私はまだ読めてません。内容の紹介はまた今度(誰かお願いします)。

細かいことはいいから使ってみる

simple featureがどうとかjustificationはいろいろあるんですが、sfはついでにつくりがモダンになっています。私がなぜsfパッケージを推すかというと理由は単純で、rgeosパッケージでもmaptoolsパッケージでも読めなかったshapefileが読めたからです。

これがそのデータです。

library(rgdal)
library(curl)

f <- tempfile(fileext = ".zip")
# OSによってdownload.file()だとエラーになるのでcurl_download()で
curl_download('http://nlftp.mlit.go.jp/ksj/gml/data/P12/P12-14/P12-14_05_GML.zip', destfile = f)

# ほんとはディレクトリを掘った方がいいけどめんどくさいのでtempdir直下に展開する
unzip(f, exdir = tempdir())

d <- readOGR(tempdir(), layer = 'P12a-14_05')
#> OGR data source with driver: ESRI Shapefile 
#> Source: "C:\Users\user1\AppData\Local\Temp\RtmpCwXYvr", layer: "P12a-14_05"
#> with 621 features
#> It has 7 fields
#> Error in readOGR(tempdir(), layer = "P12a-14_05") : 
#>   Incompatible geometry: 4
#> In addition: There were 50 or more warnings (use warnings() to see the first 50)

このエラーは、geometry: 4MULTIPOINTなんですが、それをrgdalでは扱えないので出ているエラーらしいです(ただのバグっぽい気もする)。

これを、sfでやってみます。データの読み込みにはst_read()という関数を使います(sfの関数はすべてst_というプレフィックスがついています)。

library(sf)

d <- st_read(tempdir(), layer = 'P12a-14_05')
#> Reading layer `P12a-14_05' from data source `C:\Users\user1\AppData\Local\Temp\RtmpCwXYvr\P12a-14_05.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 621 features and 7 fields
#> geometry type:  MULTIPOINT
#> dimension:      XY
#> bbox:           xmin: 139.7046 ymin: 38.95984 xmax: 140.8945 ymax: 40.46308
#> epsg (SRID):    NA
#> proj4string:    +proj=longlat +ellps=GRS80 +no_defs

あっさり読み込めました。中身を見てみましょう。

knitr::kable(head(d))
P12_001 P12_002 P12_003 P12_004 P12_005 P12_006 P12_007 geometry
10001 秋田竿燈まつり 05 05201 年中行事 秋田市旭北 1 140.10262, 39.71988
10002 秋田市大森山公園 05 05201 秋田市浜田 4 140.07247, 39.67297
10003 秋田市花木観光農園 05 05201 秋田市雄和向野字桧谷地1 6 140.23706, 39.54113
10004 秋田市浜田森林総合公園 05 05201 秋田市浜田 1 140.08242, 39.65611
10005 自然科学学習館 05 05201 秋田市東通仲町4-1 6 140.13152, 39.71584
10006 秋田市文化会館 05 05201 秋田市山王7-3-1 6 140.0980, 39.7185

右端のgeometryというのが座標情報のカラムで、あとはメタデータです。これがspのデータ形式SpatialMultiPointsDataFrameとか)だと、別々に分かれてしまっていました。sfでは、data.frameと同じ感覚で扱えます。

これの何がすごいかというと、dplyrでただのdata.frameのように扱えるということです。

library(dplyr)

# 全体を灰色でplot
plot(d, col = "grey")

# P12_004の条件が当てはまるものを赤色にplot
d %>%
  filter(P12_004 == "05211") %>%
  plot(add = TRUE, col = "red", pch = 19)

f:id:yutannihilation:20170106200711p:plain:w450

あっさりfilter()できてしまいました。spの時はこんなにあっさりいきませんでした。すごい。

filter()に使えるGIS用のpredicate function(例:点がポリゴン内にあるか)とかsummarise()に使える集計関数とかが充実して、うまく使いこなせるようになると便利そうですね。

可視化

書き忘れてましたが、sfのデータは↑のような感じでplot()するだけでプロットできます。

ggplot2は上のHadleyのツイートにあったように、次期バージョンでサポートされるらしいです。

leafletは、issueが見当たらなかったので対応がどうなるか少し不安です。以下の記事では、as(x, "Spatial")でいちどspのデータに変換してからleafletに渡していました(変換の仕方についてはsfのvignetteに書かれています)。

Spatial analysis pipelines with simple features in R · Kyle Walker

追記(2017/1/7):
id:u_riboの記事を見ていて知りましたが、この記事を書いた直後にHadleyがleafletにpull requestを出していました。ということで、心配する必要はなさそうです。

PostGIS

私はPostGISを使ったことないのでわからないんですけど、PostGISとか他のデータベースもソースとして扱えます。DBIに準拠しているようなので、このへんも透過的にdplyrのバックエンドとして扱えるようになるでしょう(まだなってないと思います)。

meuse <- st_read("PG:dbname=postgis", "meuse")

使えるデータソースはvignetteのこのあたりが参考になりそうです。

GeoJSONは?

ところで、rOpenSciもGIS関連のパッケージを多くつくっていますが、GeoJSON推しです。geojsonパッケージでは、これさえあればもうsp要らずでおk、と言わんばかりの挑発的な文句が書かれています。

geojson aims to deal only with geojson data, without requiring any of the sp/rgdal/rgeos stack. That means this package is light weight.
(GitHub - ropensci/geojson: GeoJSON classes for R)

sfパッケージについては以下のブログ記事で言及されているので、存在は認識しているようです。rOpenSciの主張は、やりたいのはただのデータ形式の変換程度であって、フルスタックなGISツールはいらないケースもけっこうあるよね、という感じのスタンスです。

The rOpenSci geospatial suite

たしかにそういうケースは結構ありそうなので、なんだかんだ共存していくような感じはします。(Turf.jsのラッパーパッケージとかはかぶってる気もしますけど)

ただ、sfはGeoJSONも読み書きできるので、ツールを一つにそろえたいという場合はsfでよさそうです。

tmp_json <- tempfile(fileext = ".geojson")
st_write(d, tmp_json)
#> Writing layer `file69429ae5b0e.geojson' to data source `C:\Users\user1\AppData\Local\Temp\RtmpCwXYvr\file69429ae5b0e.geojson' using driver `GeoJSON'
#> features:       621
#> fields:         7
#> geometry type:  MULTIPOINT

readr::read_file(tmp_json)
[1] "{\n\"type\": \"FeatureCollection\",\n\"features\": [\n{ \"type\": \"Feature\", \"properties\": { \"P12_001\": 10001, ... <truncated>

GDALとGEOSのビルドが大変そう

とはいえ、フルスタックなGISツールは大変です。具体的にはGDALとGEOSです。いろんなOSでGDALとGEOSをビルドするのがつらかった話は以下に書かれています。

そしてここでLinuxユーザの人には少し悲しいお知らせです。sfは、GDAL 2.0.0以上、GEOS 3.3.0以上というモダンな環境でしか動きません。これがどれくらいモダンかというと、Ubuntuで言えば16.10以降でなければ動きません。16.04以前はubuntugis/ubuntugis-unstableというPPAを入れないとだめらしいです。unstableと書かれるとちょっとためらいますね...。

ちなみにWindowsは、rwinlibがつくったバイナリを使ってるのでインストールするだけです。お手軽ですね(rwinlibについてはTokyo.RでrwinlibについてLTしてきました。 - Technically, technophobic.あたりに書いたのでご参照ください)。

感想

ということで、Linuxへのインストールにはまだ難があるし、関連パッケージの対応もまだまだですが、sfはRのGISの未来です。まだ少し触っただけですが、spでつらいなーと思ってたこと(データの読み込みとデータの前処理)がsfだとあっさりできて感動しました。

まだGISのこともさっぱりわからない私ですが、とりあえずvignetteを読んでもうちょっといろいろ試そうと思います。