4. シンプルフィーチャを操作する
Edzer Pebesma
For a better version of the sf vignettes see https://r-spatial.github.io/sf/articles/
この vignette では、シンプルフィーチャ（ジオメトリに付随するレコード）を操作する方法について、ジオメトリを含む場合を例に説明します。操作の内容は以下の通りです。
- フィーチャの集計
- フィーチャの要約
- フィーチャジオメトリに基づく 2 つのフィーチャセットの結合
フィーチャは sf オブジェクトのレコードとして表現され、フィーチャ属性（ジオメトリ以外のすべてのフィールド）とフィーチャジオメトリを持つ。sfオブジェクトはdata.frameやtbl_dfのサブクラスなので、フィーチャ属性の操作はdata.frame` に対して行うのと同じように行うことができる（例.
library(sf)
## Linking to GEOS 3.10.3, GDAL 3.5.1, PROJ 7.2.1; sf_use_s2() is TRUE
nc <- st_read(system.file("shape/nc.shp", package="sf"))
## Reading layer `nc' from data source
## `/Users/baba/Library/R/x86_64/4.2/library/sf/shape/nc.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 100 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## Geodetic CRS: NAD27
nc <- st_transform(nc, 2264)
nc[1,]
## Simple feature collection with 1 feature and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 1193283 ymin: 913326.4 xmax: 1340553 ymax: 1044143
## Projected CRS: NAD83 / North Carolina (ftUS)
## AREA PERIMETER CNTY_ CNTY_ID NAME FIPS FIPSNO CRESS_ID BIR74 SID74 NWBIR74
## 1 0.114 1.442 1825 1825 Ashe 37009 37009 5 1091 1 10
## BIR79 SID79 NWBIR79 geometry
## 1 1364 0 19 MULTIPOLYGON (((1270813 913...
は最初のレコードを表示します。
tidyverse/dplyr の動詞の多くは、sfオブジェクトのためのメソッドを持っています。これは、 sf と dplyr の両方がロードされている場合、1つの属性を選択するような操作は sf オブジェクトを返すことを意味します。
library(dplyr)
##
## 次のパッケージを付け加えます: 'dplyr'
## 以下のオブジェクトは 'package:stats' からマスクされています:
##
## filter, lag
## 以下のオブジェクトは 'package:base' からマスクされています:
##
## intersect, setdiff, setequal, union
nc %>% select(NWBIR74) %>% head(2)
## Simple feature collection with 2 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 1193283 ymin: 913326.4 xmax: 1441000 ymax: 1044143
## Projected CRS: NAD83 / North Carolina (ftUS)
## NWBIR74 geometry
## 1 10 MULTIPOLYGON (((1270813 913...
## 2 10 MULTIPOLYGON (((1340553 959...
これは、ジオメトリがスティッキーであり、自動的に追加されることを意味します。ジオメトリを削除しているい場合は、まず data.frame を作成し、ジオメトリのリストカラムを削除します。
nc %>% as.data.frame %>% select(NWBIR74) %>% head(2)
## NWBIR74
## 1 10
## 2 10
フィーチャセットのサブセット化
角括弧の表記でフィーチャセットをサブセットすることができる
nc[1, "NWBIR74"]
## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 1193283 ymin: 913326.4 xmax: 1340553 ymax: 1044143
## Projected CRS: NAD83 / North Carolina (ftUS)
## NWBIR74 geometry
## 1 10 MULTIPOLYGON (((1270813 913...
で、ジオメトリをドロップするには drop 引数を使用します。
nc[1, "NWBIR74", drop = TRUE]
## [1] 10
## attr(,"class")
## [1] "numeric"
しかし、空間オブジェクトを行セレクタとして使用し、別の空間フィーチャと交差するフィーチャを選択することも可能です。
Ashe = nc[nc$NAME == "Ashe",]
class(Ashe)
## [1] "sf" "data.frame"
nc[Ashe,]
## Simple feature collection with 4 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 1142157 ymin: 823077.4 xmax: 1448917 ymax: 1044143
## Projected CRS: NAD83 / North Carolina (ftUS)
## AREA PERIMETER CNTY_ CNTY_ID NAME FIPS FIPSNO CRESS_ID BIR74 SID74
## 1 0.114 1.442 1825 1825 Ashe 37009 37009 5 1091 1
## 2 0.061 1.231 1827 1827 Alleghany 37005 37005 3 487 0
## 18 0.199 1.984 1874 1874 Wilkes 37193 37193 97 3146 4
## 19 0.081 1.288 1880 1880 Watauga 37189 37189 95 1323 1
## NWBIR74 BIR79 SID79 NWBIR79 geometry
## 1 10 1364 0 19 MULTIPOLYGON (((1270813 913...
## 2 10 542 3 12 MULTIPOLYGON (((1340553 959...
## 18 200 3725 7 222 MULTIPOLYGON (((1402673 837...
## 19 17 1775 1 33 MULTIPOLYGON (((1171157 868...
結果セットには Ashe が含まれています。[.sf の引数 op のデフォルト値は st_intersects で、Ashe は自分自身と交差していることがわかります。述語 st_touches を用いることで、自己交差を除外することができます（重なり合うフィーチャは触れません）。
Ashe = nc[nc$NAME == "Ashe",]
nc[Ashe, op = st_touches]
## Simple feature collection with 3 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 1142157 ymin: 823077.4 xmax: 1448917 ymax: 1035625
## Projected CRS: NAD83 / North Carolina (ftUS)
## AREA PERIMETER CNTY_ CNTY_ID NAME FIPS FIPSNO CRESS_ID BIR74 SID74
## 2 0.061 1.231 1827 1827 Alleghany 37005 37005 3 487 0
## 18 0.199 1.984 1874 1874 Wilkes 37193 37193 97 3146 4
## 19 0.081 1.288 1880 1880 Watauga 37189 37189 95 1323 1
## NWBIR74 BIR79 SID79 NWBIR79 geometry
## 2 10 542 3 12 MULTIPOLYGON (((1340553 959...
## 18 200 3725 7 222 MULTIPOLYGON (((1402673 837...
## 19 17 1775 1 33 MULTIPOLYGON (((1171157 868...
dplyr を使えば、述語を直接呼び出すことで同じことができます。
nc %>% filter(lengths(st_touches(., Ashe)) > 0)
## Simple feature collection with 3 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 1142157 ymin: 823077.4 xmax: 1448917 ymax: 1035625
## Projected CRS: NAD83 / North Carolina (ftUS)
## AREA PERIMETER CNTY_ CNTY_ID NAME FIPS FIPSNO CRESS_ID BIR74 SID74
## 1 0.061 1.231 1827 1827 Alleghany 37005 37005 3 487 0
## 2 0.199 1.984 1874 1874 Wilkes 37193 37193 97 3146 4
## 3 0.081 1.288 1880 1880 Watauga 37189 37189 95 1323 1
## NWBIR74 BIR79 SID79 NWBIR79 geometry
## 1 10 542 3 12 MULTIPOLYGON (((1340553 959...
## 2 200 3725 7 222 MULTIPOLYGON (((1402673 837...
## 3 17 1775 1 33 MULTIPOLYGON (((1171157 868...
フィーチャセットの集計または要約
例えば、Asheと交差する郡の SID (sudden infant death) の1974年の割合を、残りの郡と比較しているいとします。これは次のようにして行うことができます。
a <- aggregate(nc[, c("SID74", "BIR74")], list(Ashe_nb = lengths(st_intersects(nc, Ashe)) > 0), sum)
(a <- a %>% mutate(frac74 = SID74 / BIR74) %>% select(frac74))
## Simple feature collection with 2 features and 1 field
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: 406265 ymin: 48359.7 xmax: 3052877 ymax: 1044143
## Projected CRS: NAD83 / North Carolina (ftUS)
## frac74 geometry
## 1 0.0020406588 MULTIPOLYGON (((454155.5 58...
## 2 0.0009922276 POLYGON ((1372051 837036.9,...
plot(a[2], col = c(grey(.8), grey(.5)))
plot(st_geometry(Ashe), border = '#ff8888', add = TRUE, lwd = 2)
属性に基づく 2 つのフィーチャセットの結合
ベース R (merge) や dplyr (left_join など) の通常の結合動詞は、 sf オブジェクトに対しても同様に動作します。ジオメトリがマッチしない場合は、空のジオメトリが代用されます。2 番目の引数は sf オブジェクトではなく、 data.frame (または類似のもの) であるべきです。
x = st_sf(a = 1:2, geom = st_sfc(st_point(c(0,0)), st_point(c(1,1))))
y = data.frame(a = 2:3)
merge(x, y)
## Simple feature collection with 1 feature and 1 field
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 1 ymin: 1 xmax: 1 ymax: 1
## CRS: NA
## a geometry
## 1 2 POINT (1 1)
merge(x, y, all = TRUE)
## Simple feature collection with 3 features and 1 field (with 1 geometry empty)
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: 0 ymin: 0 xmax: 1 ymax: 1
## CRS: NA
## a geometry
## 1 1 POINT (0 0)
## 2 2 POINT (1 1)
## 3 3 GEOMETRYCOLLECTION EMPTY
right_join(x, y)
## Joining, by = "a"
## Simple feature collection with 2 features and 1 field (with 1 geometry empty)
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 1 ymin: 1 xmax: 1 ymax: 1
## CRS: NA
## a geom
## 1 2 POINT (1 1)
## 2 3 POINT EMPTY
幾何学に基づく 2 つのフィーチャセットの結合
空間的な交点に基づく結合には、（あらゆる種類の） st_join が使用されます。
x = st_sf(a = 1:3, geom = st_sfc(st_point(c(1,1)), st_point(c(2,2)), st_point(c(3,3))))
y = st_buffer(x, 0.1)
x = x[1:2,]
y = y[2:3,]
plot(st_geometry(x), xlim = c(.5, 3.5))
plot(st_geometry(y), add = TRUE)
結合方法は左結合で、最初の属性のすべてのレコードを保持します。
st_join(x, y)
## Simple feature collection with 2 features and 2 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 1 ymin: 1 xmax: 2 ymax: 2
## CRS: NA
## a.x a.y geom
## 1 1 NA POINT (1 1)
## 2 2 2 POINT (2 2)
st_join(y, x)
## Simple feature collection with 2 features and 2 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 1.9 ymin: 1.9 xmax: 3.1 ymax: 3.1
## CRS: NA
## a.x a.y geom
## 2 2 2 POLYGON ((2.1 2, 2.099863 1...
## 3 3 NA POLYGON ((3.1 3, 3.099863 2...
であり、保持されるジオメトリは第1引数のものです。
空間結合述語は st_intersects （デフォルト）と互換性のある任意の関数で制御することができます。
st_join(x, y, join = st_covers) # 一致するYレコードがない: 点が円をカバーしない
## Simple feature collection with 2 features and 2 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 1 ymin: 1 xmax: 2 ymax: 2
## CRS: NA
## a.x a.y geom
## 1 1 NA POINT (1 1)
## 2 2 NA POINT (2 2)
st_join(y, x, join = st_covers) # 点を含む円にマッチ
## Simple feature collection with 2 features and 2 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 1.9 ymin: 1.9 xmax: 3.1 ymax: 3.1
## CRS: NA
## a.x a.y geom
## 2 2 2 POLYGON ((2.1 2, 2.099863 1...
## 3 3 NA POLYGON ((3.1 3, 3.099863 2...