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
