lwgeomパッケージでRからliblwgeomの関数を使う

R GIS sf lwgeom

lwgeomは、liblwgeomというライブラリをRから使えるようにするパッケージです。

次のような関数が用意されています。

  • st_make_valid(): 不正な地物を整形
  • st_geohash(): geohashを計算
  • st_geod_*(): 面積や方位を計算
  • st_transform_proj: 投影法の変更
  • st_minimum_bounding_circle(): minimum bounding circle（境界円？）を計算

sfとの関係

次のようにsfとlwgeomを読み込んでみると分かりますが、いくつかの関数は名前が被っています。また、st_transform_proj()にはst_transform()、st_geod_*()にはst_*()という似たような関数があります。なぜこんなことになっているのでしょう？

library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.0, proj.4 4.9.3
library(lwgeom)
#> Linking to liblwgeom 2.5.0dev r16016
#> 
#> Attaching package: 'lwgeom'
#> The following objects are masked from 'package:sf':
#> 
#>     st_geohash, st_make_valid, st_split

実は、同名の関数はliblwgeomというライブラリに依存していて、それを使ってsfをビルドする設定にしていないと使えない関数です。 私の手元のWindowsではエラーになりました。

p1_valid <- sf::st_make_valid(p1)
#> Error in CPL_make_valid(x): st_make_valid requires compilation against liblwgeom

どうも、liblwgeomをsfに含めるのは大変だから、ということで別個のパッケージとしてlwgeomがつくられた、という経緯のようです。

st_make_valid()

上に書いたように、sfにはすでにst_make_valid()があります。これに加えて、地物が不正な形ではないか判定するst_is_valid()という関数があります（こっちは動く）。

まず、仮に不正なPOLYGONを作ってみましょう。

p1 <- st_as_sfc("POLYGON((0 0, 0 10, 10 0, 10 10, 0 0))")

p1
#> Geometry set for 1 feature 
#> geometry type:  POLYGON
#> dimension:      XY
#> bbox:           xmin: 0 ymin: 0 xmax: 10 ymax: 10
#> epsg (SRID):    NA
#> proj4string:    NA
#> POLYGON ((0 0, 0 10, 10 0, 10 10, 0 0))

plot()してみると分かりますが、線が交差してしまっています。

plot(p1)

このため、st_is_valid()で判定すると不正であることがわかります。

st_is_valid(p1)
#> Warning in evalq((function (..., call. = TRUE, immediate. = FALSE,
#> noBreaks. = FALSE, : Self-intersection at or near point 5 5
#> [1] FALSE

ここでst_make_valid()の出番です。sfのはエラーになりましたが、lwgeomの方の関数はちゃんと動きます。

p1_valid <- lwgeom::st_make_valid(p1)

p1_valid
#> Geometry set for 1 feature 
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: 0 ymin: 0 xmax: 10 ymax: 10
#> epsg (SRID):    NA
#> proj4string:    NA
#> MULTIPOLYGON (((5 5, 10 10, 10 0, 5 5)), ((0 0,...

見た目的には変わっていませんが、交差しているところで分割されて、正しいMULTIPOLYGONになりました。

plot(p1_valid)

st_geod_*()

さて、同名の関数は上に書いたようなことなんですが、名前が違う関数はよくわかりませんでした。

例えばこれはst_geod_area()のExamplesに乗っている例なんですが、st_area()と結果が違います。地表面を平面とみなすか球面とみなすかの違い、とか...？ 詳しい人教えてください。

nc <- st_read(system.file("gpkg/nc.gpkg", package="sf"))
#> Reading layer `nc.gpkg' from data source `C:\Users\hiroaki-yutani\Documents\R\win-library\3.4\sf\gpkg\nc.gpkg' using driver `GPKG'
#> Simple feature collection with 100 features and 14 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> epsg (SRID):    4267
#> proj4string:    +proj=longlat +datum=NAD27 +no_defs
st_geod_area(nc[1:3,])
#> Units: m^2
#> [1] 1137389166  611077451 1423490699
st_area(nc[1:3,])
#> Units: m^2
#> [1] 1137388604  611077263 1423489919