sf vignette 6

r sf

Rのsfパッケージの vignette6について、 これまでのvignetteで書けなかったことなどを雑多に書いてあるので、 あまりまとまりがない感じである。

EPSGとは

EPSGは空間座標参照システムとしてよく知られており、 IOGP (International Association of Oil & Gass Producers)によってメンテナンスされている。 下記の関数を利用することでEPSG一覧を見ることができる。

rgdal::make_EPSG()

sfで第二のgeometry columnを扱うには

sfオブジェクトは一つ以上のgeometry list-columnを持つことができるが、有効となるのは一つだけであり、 st_geometry()関数によって取得できる。 printメソッドを利用することでgeometry columnが有効か確認することができる。

library(sf)
demo(nc, ask =F, echo = F)

nc$geom2 = st_centroid(st_geometry(nc))
print(nc, n = 2)

Simple feature collection with 100 features and 14 fields
Attribute-geometry relationship: 0 constant, 8 aggregate, 6 identity, 1 NA's
Active geometry column: geom
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
First 2 features:
   AREA PERIMETER CNTY_ CNTY_ID      NAME  FIPS FIPSNO CRESS_ID BIR74
1 0.114     1.442  1825    1825      Ashe 37009  37009        5  1091
2 0.061     1.231  1827    1827 Alleghany 37005  37005        3   487
  SID74 NWBIR74 BIR79 SID79 NWBIR79                           geom
1     1      10  1364     0      19 MULTIPOLYGON (((-81.47276 3...
2     0      10   542     3      12 MULTIPOLYGON (((-81.23989 3...
                       geom2
1  POINT (-81.49826 36.4314)
2 POINT (-81.12515 36.49101)

geometry columnを変更するにはst_geometry<-またはst_set_geometry関数を利用する。

plot(st_geometry(nc))

f:id:jerrarrdan:20180224104751p:plain

st_geometry(nc) <- "geom2"
plot(st_geometry(nc))

f:id:jerrarrdan:20180224104912p:plain

st_simplifyでのトポロジー関係の維持について

st_simpligyはトポロジー保持する関数であるが、これは個々のfeature geometrieで行われる。 つまり、ポリゴンに対して実行したものはポリゴンとなる。 しかし、2つのfeaturesが境界を共有しているときにst_simplifyを実行した場合には2つのポリゴンが同じ境界を共有しているとは限らない。これは個々のポリゴンごとに単純化が行われるためである。

sfオブジェクトに対してdplyr verbsが使えない

多くの開発者はパッケージをロードせずに、sf::のように関数の属するパッケージを明示して下記のように記述する。

i = sf::st_intersects(sf1, sf2)

しかし、sfオブジェクトに対して、dplyr::selectは利用できません。 このケースではsfパッケージをロードする必要がある。

表示されるerror/warning/messageについて

座標系は緯度経度であるが xxxは平面を仮定しています（although coordinates are longitude/latitude, xxx assumes that they are planar）

sfで利用されるgeometryに関する計算においてはGEOSライブラリを利用している。このライブラリは2次元のフラットなユクリッド距離の空間座標を考慮する。しかし、例えば極を含む緯度経度データを扱い場合にはこれだけでは不十分です。

下記に例を示すが、これは間違った結果を返すこととなる。

polygon = st_sfc(
  st_polygon(
    list(rbind(c(0, 80), c(120, 80),
               c(240, 80), c(0, 80)))), 
  crs =4326)
pole = st_sfc(st_point(c(0, 90)), crs = 4326)
st_intersects(polygon, pole)

Sparse geometry binary predicate list of length 1, where the predicate was `intersects'
 1: (empty)

st_centroidは緯度経度データに対しては正しい値を返しません(

st_centroid does not give correct centroids for longitude/latitude data)

先ほどの例と同じですが、中心は2次元空間を前提としています。

st_centroid(polygon)[[1]]

この場合中心は極である必要があります。

距離は10進数を前提としています(

dist is assumed to be in decimal degrees (arc_degrees).)

このメッセージはsfでは距離の値が度で与えられていることを示しています。このメッセージを避けるには正しいユニット（単位）を与えてやる必要があります。

pt = st_sfc(st_point(c(0,0)), crs = 4326)
buf = st_buffer(polygon, 1)
buf = st_buffer(polygon, units::set_units(1, degree))

両方の式のオペランドがunitオブジェクトもしくはアトミックベクトルである必要があります(

both operands of the expression should be "units" objects or _Error: $ operator is invalid for atomic vectors")

このメッセージはsfパッケージをロードした後にspatstatパッケージをロードした時に現れる。パッケージ単位で作成された単位オブジェクトの場合、unitクラスも実装するspatstatのメソッドが選択される。この問題はspatstatの後にunitパッケージをロードすることで解決できる。なお、バージョン1.53-2.012のspatstatパッケージではこれが解決されています。