3. シンプルフィーチャの形状を操作する

Edzer Pebesma

For a better version of the sf vignettes see https://r-spatial.github.io/sf/articles/

この vignette では、シンプルフィーチャ（地物）の形状を操作する方法を説明します。

ジオメトリタイプの変換

このセクションでは、あるタイプのシンプルフィーチャの形状を別のタイプに変換する方法について説明します。線をポリゴンに変換する方法については、以下の st_polygonize も参照してください。

単一ジオメトリの場合

単一ジオメトリの場合 st_cast は

  1. XX から MULTIXX に変換します。例えば LINESTRING から MULTILINESTRING に変換します。
  2. MULTIXX が長さ1の場合、MULTIXX から XX に変換する（それ以外の場合、変換は行われるが情報の損失について警告が出される）。
  3. MULTIXX が長さ1でない場合、MULTIXX から XX に変換するが、情報の損失について警告します。
  4. 長さ1の GEOMETRYCOLLECTION をその構成要素に変換します。

最初の3種類の例としては:

library(sf)
## Linking to GEOS 3.10.3, GDAL 3.5.1, PROJ 7.2.1; sf_use_s2() is TRUE
suppressPackageStartupMessages(library(dplyr))
st_point(c(1,1)) %>% st_cast("MULTIPOINT")
## MULTIPOINT ((1 1))
st_multipoint(rbind(c(1,1))) %>% st_cast("POINT")
## Warning in st_cast.MULTIPOINT(., "POINT"): point from first coordinate only
## POINT (1 1)
st_multipoint(rbind(c(1,1),c(2,2))) %>% st_cast("POINT")
## Warning in st_cast.MULTIPOINT(., "POINT"): point from first coordinate only
## POINT (1 1)

第4のタイプの例としては:

st_geometrycollection(list(st_point(c(1,1)))) %>% st_cast("POINT")
## POINT (1 1)

ジオメトリコレクション（sfc）、シンプルフィーチャコレクション（sf）用

ここで注意すべきは、 st_read を使用してジオメトリを読み込む際に、type 引数を使用して返されるジオメトリのクラスを制御することができることです。

shp = system.file("shape/nc.shp", package="sf")
class(st_geometry(st_read(shp, quiet = TRUE)))
## [1] "sfc_MULTIPOLYGON" "sfc"
class(st_geometry(st_read(shp, quiet = TRUE, type = 3)))
## [1] "sfc_POLYGON" "sfc"
class(st_geometry(st_read(shp, quiet = TRUE, type = 1)))
## [1] "sfc_GEOMETRY" "sfc"

このオプションは GDAL ライブラリによって処理されます。ターゲット型への変換に失敗する場合には、元の型を返します。この場合、POLYGON と MULTIPOLYGON の幾何学が混在し、スーパークラスとして GEOMETRY を導くことになります。マルチポリゴンをポリゴンとして読もうとすると、マルチポリゴンのすべての二次リングが失われてしまいます。

関数がジオメトリタイプが混在しているオブジェクト (GEOMETRY) を返す場合、st_write などの下流の関数では扱いが難しくなることがあります。このような場合、st_cast を使用すると、オブジェクトの型を変更することができます。 ジオメトリオブジェクトのセット (sfc) やシンプルフィーチャセット (sf) に対して、st_cast は対象の型を指定しても、指定しなくても使用することができます。

ls <- st_linestring(rbind(c(0,0),c(1,1),c(2,1)))
mls <- st_multilinestring(list(rbind(c(2,2),c(1,3)), rbind(c(0,0),c(1,1),c(2,1))))
(sfc <- st_sfc(ls,mls))
## Geometry set for 2 features 
## Geometry type: GEOMETRY
## Dimension:     XY
## Bounding box:  xmin: 0 ymin: 0 xmax: 2 ymax: 3
## CRS:           NA
## LINESTRING (0 0, 1 1, 2 1)
## MULTILINESTRING ((2 2, 1 3), (0 0, 1 1, 2 1))
st_cast(sfc, "MULTILINESTRING")
## Geometry set for 2 features 
## Geometry type: MULTILINESTRING
## Dimension:     XY
## Bounding box:  xmin: 0 ymin: 0 xmax: 2 ymax: 3
## CRS:           NA
## MULTILINESTRING ((0 0, 1 1, 2 1))
## MULTILINESTRING ((2 2, 1 3), (0 0, 1 1, 2 1))
sf <- st_sf(a = 5:4, geom = sfc)
st_cast(sf, "MULTILINESTRING")
## Simple feature collection with 2 features and 1 field
## Geometry type: MULTILINESTRING
## Dimension:     XY
## Bounding box:  xmin: 0 ymin: 0 xmax: 2 ymax: 3
## CRS:           NA
##   a                           geom
## 1 5 MULTILINESTRING ((0 0, 1 1,...
## 2 4 MULTILINESTRING ((2 2, 1 3)...

対象となる型が与えられていない場合、 st_cast は 2 つのケースに対して賢くなろうとします。

  1. オブジェクトのクラスが GEOMETRY で、すべての要素が同じ型である場合、および
  2. すべての要素が長さ1の GEOMETRYCOLLECTION オブジェクトである場合、GEOMETRYCOLLECTION オブジェクトはその内容（再び GEOMETRY の混合物になる可能性もある）に置き換えられます。

例としては、以下のようなものがあります。

ls <- st_linestring(rbind(c(0,0),c(1,1),c(2,1)))
mls1 <- st_multilinestring(list(rbind(c(2,2),c(1,3)), rbind(c(0,0),c(1,1),c(2,1))))
mls2 <- st_multilinestring(list(rbind(c(4,4),c(4,3)), rbind(c(2,2),c(2,1),c(3,1))))
(sfc <- st_sfc(ls,mls1,mls2))
## Geometry set for 3 features 
## Geometry type: GEOMETRY
## Dimension:     XY
## Bounding box:  xmin: 0 ymin: 0 xmax: 4 ymax: 4
## CRS:           NA
## LINESTRING (0 0, 1 1, 2 1)
## MULTILINESTRING ((2 2, 1 3), (0 0, 1 1, 2 1))
## MULTILINESTRING ((4 4, 4 3), (2 2, 2 1, 3 1))
class(sfc[2:3])
## [1] "sfc_MULTILINESTRING" "sfc"
class(st_cast(sfc[2:3]))
## [1] "sfc_MULTILINESTRING" "sfc"

gc1 <- st_geometrycollection(list(st_linestring(rbind(c(0,0),c(1,1),c(2,1)))))
gc2 <- st_geometrycollection(list(st_multilinestring(list(rbind(c(2,2),c(1,3)), rbind(c(0,0),c(1,1),c(2,1))))))
gc3 <- st_geometrycollection(list(st_multilinestring(list(rbind(c(4,4),c(4,3)), rbind(c(2,2),c(2,1),c(3,1))))))
(sfc <- st_sfc(gc1,gc2,gc3))
## Geometry set for 3 features 
## Geometry type: GEOMETRYCOLLECTION
## Dimension:     XY
## Bounding box:  xmin: 0 ymin: 0 xmax: 4 ymax: 4
## CRS:           NA
## GEOMETRYCOLLECTION (LINESTRING (0 0, 1 1, 2 1))
## GEOMETRYCOLLECTION (MULTILINESTRING ((2 2, 1 3)...
## GEOMETRYCOLLECTION (MULTILINESTRING ((4 4, 4 3)...
class(st_cast(sfc))
## [1] "sfc_GEOMETRY" "sfc"
class(st_cast(st_cast(sfc), "MULTILINESTRING"))
## [1] "sfc_MULTILINESTRING" "sfc"

アフィン変換

アフィン変換とは、\(f(x) = xA + b\)という型の変換で、行列 \(A\) は平坦化、拡大縮小、回転に使われ、\(b\) は \(x\) を平行移動させます。低レベルの例としては

(p = st_point(c(0,2)))
## POINT (0 2)
p + 1
## POINT (1 3)
p + c(1,2)
## POINT (1 4)
p + p
## POINT (0 4)
p * p
## POINT (0 4)
rot = function(a) matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2)
p * rot(pi/4)
## POINT (1.414214 1.414214)
p * rot(pi/2)
## POINT (2 1.224647e-16)
p * rot(pi)
## POINT (2.449294e-16 -2)

例えば、ノースカロライナ州の郡を重心を中心に時計回りに90度回転させ、元のサイズの75%に縮小することができます。

nc = st_read(system.file("shape/nc.shp", package="sf"), quiet = TRUE)
ncg = st_geometry(nc)
plot(ncg, border = 'grey')
cntrd = st_centroid(ncg)
ncg2 = (ncg - cntrd) * rot(pi/2) * .75 + cntrd
plot(ncg2, add = TRUE)
plot(cntrd, col = 'red', add = TRUE, cex = .5)

座標参照系の変換と変形

sf オブジェクトの座標参照系の取得と設定

クラス sf または sfc のオブジェクトの座標参照系は次のとおりです。 は st_crs によって取得され、st_crs<- によって置き換えられます。

library(sf)
geom = st_sfc(st_point(c(0,1)), st_point(c(11,12)))
s = st_sf(a = 15:16, geometry = geom)
st_crs(s)
## Coordinate Reference System: NA
s1 = s
st_crs(s1) <- 4326
st_crs(s1)
## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]
s2 = s
st_crs(s2) <- "+proj=longlat +datum=WGS84"
all.equal(s1, s2)
## [1] "Component \"geometry\": Attributes: < Component \"crs\": Component \"input\": 1 string mismatch >"
## [2] "Component \"geometry\": Attributes: < Component \"crs\": Component \"wkt\": 1 string mismatch >"

よりパイプに優しい代替の st_crs<- は以下の通りです。

s1 %>% st_set_crs(4326)
## Simple feature collection with 2 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 0 ymin: 1 xmax: 11 ymax: 12
## Geodetic CRS:  WGS 84
##    a      geometry
## 1 15   POINT (0 1)
## 2 16 POINT (11 12)

座標参照系の変換

座標参照系をある非欠測値から別の非欠測値に変更すると、座標は変更されずにCRSが変更されますが、値を再投影しなかいているという警告が出されます。

s3 <- s1 %>% st_set_crs(4326) %>% st_set_crs(3857)
## Warning: st_crs<- : replacing crs does not reproject data; use st_transform for
## that

このような警告が発生しない、より適切な方法は、まずCRSに欠損値を割り当てて消去し、その後、意図している値に設定することです。

s3 <- s1  %>% st_set_crs(NA) %>% st_set_crs(3857)

座標変換を行うには、 st_transform を使用します。

s3 <- s1 %>% st_transform(3857)
s3
## Simple feature collection with 2 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 0 ymin: 111325.1 xmax: 1224514 ymax: 1345708
## Projected CRS: WGS 84 / Pseudo-Mercator
##    a                geometry
## 1 15      POINT (0 111325.1)
## 2 16 POINT (1224514 1345708)

この場合、実際に座標が変更（投影）されることがわかります。

幾何学的な操作

すべての幾何演算 st_op(x) または st_op2(x,y) は sf オブジェクトと sfc オブジェクト x と y の両方で動作します。この演算はジオメトリに対して行われるので、sf オブジェクトの非ジオメトリ部分は単に破棄されます。また、 st_op2(x,y) のような単一の引数で呼ばれるすべての二項演算は st_op2(x,x) として扱われます。

ここでは、非常に簡単なデータセットを用いて、幾何学的な操作を説明します。

b0 = st_polygon(list(rbind(c(-1,-1), c(1,-1), c(1,1), c(-1,1), c(-1,-1))))
b1 = b0 + 2
b2 = b0 + c(-0.2, 2)
x = st_sfc(b0, b1, b2)
a0 = b0 * 0.8
a1 = a0 * 0.5 + c(2, 0.7)
a2 = a0 + 1
a3 = b0 * 0.5 + c(2, -0.5)
y = st_sfc(a0,a1,a2,a3)
plot(x, border = 'red')
plot(y, border = 'green', add = TRUE)

単項演算

st_is_valid はポリゴン形状が位相的に有効であるかどうかを返す。

b0 = st_polygon(list(rbind(c(-1,-1), c(1,-1), c(1,1), c(-1,1), c(-1,-1))))
b1 = st_polygon(list(rbind(c(-1,-1), c(1,-1), c(1,1), c(0,-1), c(-1,-1))))
st_is_valid(st_sfc(b0,b1))
## [1]  TRUE FALSE

と st_is_simple は、線状のジオメトリが単純であるかどうかを示します。

s = st_sfc(st_linestring(rbind(c(0,0), c(1,1))), 
    st_linestring(rbind(c(0,0), c(1,1),c(0,1),c(1,0))))
st_is_simple(s)
## [1]  TRUE FALSE

また、 st_area はポリゴン形状の面積を、st_length は線状形状の長さを返す。

st_area(x)
## [1] 4 4 4
st_area(st_sfc(st_point(c(0,0))))
## [1] 0
st_length(st_sfc(st_linestring(rbind(c(0,0),c(1,1),c(1,2))), st_linestring(rbind(c(0,0),c(1,0)))))
## [1] 2.414214 1.000000
st_length(st_sfc(st_multilinestring(list(rbind(c(0,0),c(1,1),c(1,2))),rbind(c(0,0),c(1,0))))) # ignores 2nd part!
## [1] 2.414214

二項演算：距離と関係

st_distance は、ジオメトリ間の最短距離行列を求めます。これは密な行列です。

st_distance(x,y)
##           [,1] [,2] [,3]     [,4]
## [1,] 0.0000000  0.6    0 0.500000
## [2,] 0.2828427  0.0    0 1.000000
## [3,] 0.2000000  0.8    0 1.220656

st_relate は、ジオメトリの各ペア間の DE9-IM 関係を表す密な文字行列を返します。

st_relate(x,y)
##      [,1]        [,2]        [,3]        [,4]       
## [1,] "212FF1FF2" "FF2FF1212" "212101212" "FF2FF1212"
## [2,] "FF2FF1212" "212101212" "212101212" "FF2FF1212"
## [3,] "FF2FF1212" "FF2FF1212" "212101212" "FF2FF1212"

この行列の要素 [i,j] は9文字で、x[i] と y[j] の関係を表し、\(I_xI_y,I_xB_y,I_xE_y, B_xI_y,B_xB_y,B_xE_y,E_xI_y,E_xB_y,E_xE_y\) （\(I\)が内部、\(B\)が境界、\(E\)が外部）とコード化されています。 g. \(B_xI_y\) は x[i] の境界 \(B_x\) と y[j] の内部 \(I_y\) の交点の次元数で、それぞれ {0,1,2,F} のいずれかで、0次元、1次元、2次元の交差、および (F) 交差なし、を示します。

2項論理演算

二項論理演算は、疎な行列を返す

st_intersects(x,y)
## Sparse geometry binary predicate list of length 3, where the predicate
## was `intersects'
##  1: 1, 3
##  2: 2, 3
##  3: 3

or a dense matrix

st_intersects(x, x, sparse = FALSE)
##      [,1]  [,2]  [,3]
## [1,] TRUE  TRUE  TRUE
## [2,] TRUE  TRUE FALSE
## [3,] TRUE FALSE  TRUE
st_intersects(x, y, sparse = FALSE)
##       [,1]  [,2] [,3]  [,4]
## [1,]  TRUE FALSE TRUE FALSE
## [2,] FALSE  TRUE TRUE FALSE
## [3,] FALSE FALSE TRUE FALSE

ここで、疎な行列のリスト要素 i は、密な行列の行 i に含まれる TRUE 要素のインデックスを表します。大きなジオメトリ集合では、密な行列は多くのメモリを消費し、そのほとんどが FALSE 値で埋め尽くされます。しているがって、デフォルトでは、疎な行列を返します。

st_intersects は、ジオメトリのペアごとに、それらが交差しているかどうか（密な行列）、またはどの要素が交差しているか（疎な行列）を返します。このパッケージの関数 st_intersection は、st_intersects のような論理値ではなく、交差するジオメトリを返すことに注意してください（このビネットの次のセクションを参照してください）。

他の二項述語は以下の通りである（可読性のためにsparseを使用する）。

st_disjoint(x, y, sparse = FALSE)
##       [,1]  [,2]  [,3] [,4]
## [1,] FALSE  TRUE FALSE TRUE
## [2,]  TRUE FALSE FALSE TRUE
## [3,]  TRUE  TRUE FALSE TRUE
st_touches(x, y, sparse = FALSE)
##       [,1]  [,2]  [,3]  [,4]
## [1,] FALSE FALSE FALSE FALSE
## [2,] FALSE FALSE FALSE FALSE
## [3,] FALSE FALSE FALSE FALSE
st_crosses(s, s, sparse = FALSE)
##       [,1]  [,2]
## [1,] FALSE FALSE
## [2,] FALSE FALSE
st_within(x, y, sparse = FALSE)
##       [,1]  [,2]  [,3]  [,4]
## [1,] FALSE FALSE FALSE FALSE
## [2,] FALSE FALSE FALSE FALSE
## [3,] FALSE FALSE FALSE FALSE
st_contains(x, y, sparse = FALSE)
##       [,1]  [,2]  [,3]  [,4]
## [1,]  TRUE FALSE FALSE FALSE
## [2,] FALSE FALSE FALSE FALSE
## [3,] FALSE FALSE FALSE FALSE
st_overlaps(x, y, sparse = FALSE)
##       [,1]  [,2] [,3]  [,4]
## [1,] FALSE FALSE TRUE FALSE
## [2,] FALSE  TRUE TRUE FALSE
## [3,] FALSE FALSE TRUE FALSE
st_equals(x, y, sparse = FALSE)
##       [,1]  [,2]  [,3]  [,4]
## [1,] FALSE FALSE FALSE FALSE
## [2,] FALSE FALSE FALSE FALSE
## [3,] FALSE FALSE FALSE FALSE
st_covers(x, y, sparse = FALSE)
##       [,1]  [,2]  [,3]  [,4]
## [1,]  TRUE FALSE FALSE FALSE
## [2,] FALSE FALSE FALSE FALSE
## [3,] FALSE FALSE FALSE FALSE
st_covered_by(x, y, sparse = FALSE)
##       [,1]  [,2]  [,3]  [,4]
## [1,] FALSE FALSE FALSE FALSE
## [2,] FALSE FALSE FALSE FALSE
## [3,] FALSE FALSE FALSE FALSE
st_covered_by(y, y, sparse = FALSE)
##       [,1]  [,2]  [,3]  [,4]
## [1,]  TRUE FALSE FALSE FALSE
## [2,] FALSE  TRUE FALSE FALSE
## [3,] FALSE FALSE  TRUE FALSE
## [4,] FALSE FALSE FALSE  TRUE
st_equals_exact(x, y,0.001, sparse = FALSE)
##       [,1]  [,2]  [,3]  [,4]
## [1,] FALSE FALSE FALSE FALSE
## [2,] FALSE FALSE FALSE FALSE
## [3,] FALSE FALSE FALSE FALSE

ジオメトリを返す演算

u = st_union(x)
plot(u)

par(mfrow=c(1,2), mar = rep(0,4))
plot(st_buffer(u, 0.2))
plot(u, border = 'red', add = TRUE)
plot(st_buffer(u, 0.2), border = 'grey')
plot(u, border = 'red', add = TRUE)
plot(st_buffer(u, -0.2), add = TRUE)

plot(st_boundary(x))

par(mfrow = c(1:2))
plot(st_convex_hull(x))
plot(st_convex_hull(u))

par(mfrow = c(1,1))
par(mfrow=c(1,2))
plot(x)
plot(st_centroid(x), add = TRUE, col = 'red')
plot(x)
plot(st_centroid(u), add = TRUE, col = 'red')

2つのジオメトリの交差点は、両者によってカバーされるジオメトリであり、 st_intersection によって得られます。

plot(x)
plot(y, add = TRUE)
plot(st_intersection(st_union(x),st_union(y)), add = TRUE, col = 'red')

関数 st_intersects は、各ジオメトリペアが交差しているかどうかを示す論理行列を返すことに注意してください（このビネットの前のセクションを参照してください）。

交差以外のすべての情報を取得するには、st_difference または st_sym_difference を使用してください。

par(mfrow=c(2,2), mar = c(0,0,1,0))
plot(x, col = '#ff333388'); 
plot(y, add=TRUE, col='#33ff3388')
title("x: red, y: green")
plot(x, border = 'grey')
plot(st_difference(st_union(x),st_union(y)), col = 'lightblue', add = TRUE)
title("difference(x,y)")
plot(x, border = 'grey')
plot(st_difference(st_union(y),st_union(x)), col = 'lightblue', add = TRUE)
title("difference(y,x)")
plot(x, border = 'grey')
plot(st_sym_difference(st_union(y),st_union(x)), col = 'lightblue', add = TRUE)
title("sym_difference(x,y)")

関数 st_segmentize は直線やポリゴンの直線部分に点を追加します。

par(mfrow=c(1,3),mar=c(1,1,0,0))
pts = rbind(c(0,0),c(1,0),c(2,1),c(3,1))
ls = st_linestring(pts)
plot(ls)
points(pts)
ls.seg = st_segmentize(ls, 0.3)
plot(ls.seg)
pts = ls.seg
points(pts)
pol = st_polygon(list(rbind(c(0,0),c(1,0),c(1,1),c(0,1),c(0,0))))
pol.seg = st_segmentize(pol, 0.3)
plot(pol.seg, col = 'grey')
points(pol.seg[[1]])

関数 st_polygonize は、ポイントが閉じた多角形を形成する限り、マルチリンストリングをポリゴン化します。

par(mfrow=c(1,2),mar=c(0,0,1,0))
mls = st_multilinestring(list(matrix(c(0,0,0,1,1,1,0,0),,2,byrow=TRUE)))
x = st_polygonize(mls)
plot(mls, col = 'grey')
title("multilinestring")
plot(x, col = 'grey')
title("polygon")