R のシンプルフィーチャはどのように構成されるか
パッケージ sf は、シンプルフィーチャを R のネイティブオブジェクトとして表現します。 PostGIS と同様に、空間データを操作する sf のすべての関数とメソッドには、spatial type を意味する st_ がプレフィックスとして付加されています。 シンプルフィーチャは、単純なデータ構造 (S3 クラス、リスト、行列、ベクトル) を用いて、R ネイティブデータとして実装されています。 典型的な使用方法は、属性とジオメトリを持つフィーチャのセットの読み込み、操作、書き込みです。
属性は通常 data.frame オブジェクト（または非常に類似している tbl_df ）に格納されるので、フィーチャのジオメトリも data.frame カラムに格納することにします。ジオメトリは単一値ではないので、リストカラムに格納されます。リストカラムは data.frame のレコード数と同じ長さで、各リスト要素にそのフィーチャのシンプルフィーチャのジオメトリが格納されます。 シンプルフィーチャの表現に使用されるクラスは以下の3つです。
sfは、フィーチャの属性とフィーチャジオメトリを含むテーブル (data.frame) であり、以下のものを含みます： *sfは、フィーチャの属性とフィーチャジオメトリを含むテーブル (data.frame) です。- 各フィーチャ（レコード）の形状を表すリストカラムで、以下の要素から構成されます。
sfgは、個々のシンプルフィーチャのフィーチャ形状です。
以下、これら 3 つのクラスについてそれぞれ説明します。
sf: シンプルフィーチャを持つオブジェクト
通常、単一のシンプルフィーチャからなるジオメトリを扱うことはなく、属性付きのフィーチャの集合からなるデータセットを扱うため、この2つは sf (シンプルフィーチャ) オブジェクトにまとめられます。 次のコマンドは nc データセットを読み込みます。 sf パッケージを使用します。
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(ユーザーは system.file を使わず、直接 filename を指定します。また、シェープファイルは同じディレクトリにある、同じベースネームの複数のファイルから構成されていることに注意してください)。 表示されている短いレポートには、ファイル名、ドライバ (ESRI Shapefile)、100 個のフィーチャ (行で表されるレコード) と 14 個のフィールド (列で表される属性) があることが記載されています。このオブジェクトはクラス
class(nc)
## [1] "sf" "data.frame"は、data.frameを拡張している（そして「ある」）ことを意味しますが、以下の名前のカラムに、ジオメトリのリスト列を伴います。
attr(nc, "sf_column")
## [1] "geometry"最初の3つのフィーチャを表示すると、その属性値とジオメトリの簡略版をみることができます。
print(nc[9:15], n = 3)以下のように出力されます:
出力では、以下のようになります。
- 緑色の部分はシンプルフィーチャ：属性とジオメトリからなる 1 つのレコード、つまり
data.frame行です。 - 青色はシンプルフィーチャジオメトリ（クラス
sfgのオブジェクト）。 - 赤字は単純なリストカラム (
sfcクラスのオブジェクト。data.frameのカラム) - ジオメトリは R のネイティブオブジェクトですが、よく知られたテキストとして表示されている
訳注: wkt では？wkbは、「よく知られたバイナリ」。wkt については、下に解説あり。
sf オブジェクトのメソッドは、
methods(class = "sf")
## [1] [ [[<-
## [3] $<- aggregate
## [5] as.data.frame cbind
## [7] coerce dbDataType
## [9] dbWriteTable filter
## [11] identify initialize
## [13] merge plot
## [15] print rbind
## [17] show slotsFromS3
## [19] st_agr st_agr<-
## [21] st_area st_as_s2
## [23] st_as_sf st_as_sfc
## [25] st_bbox st_boundary
## [27] st_buffer st_cast
## [29] st_centroid st_collection_extract
## [31] st_convex_hull st_coordinates
## [33] st_crop st_crs
## [35] st_crs<- st_difference
## [37] st_drop_geometry st_filter
## [39] st_geometry st_geometry<-
## [41] st_inscribed_circle st_interpolate_aw
## [43] st_intersection st_intersects
## [45] st_is_valid st_is
## [47] st_join st_line_merge
## [49] st_m_range st_make_valid
## [51] st_minimum_rotated_rectangle st_nearest_points
## [53] st_node st_normalize
## [55] st_point_on_surface st_polygonize
## [57] st_precision st_reverse
## [59] st_sample st_segmentize
## [61] st_set_precision st_shift_longitude
## [63] st_simplify st_snap
## [65] st_sym_difference st_transform
## [67] st_triangulate st_union
## [69] st_voronoi st_wrap_dateline
## [71] st_write st_z_range
## [73] st_zm transform
## see '?methods' for accessing help and source codeまた、クラス sf ではないジオメトリのリストカラムを持つ data.frame オブジェクトを作成することも可能です。
nc.no_sf <- as.data.frame(nc)
class(nc.no_sf)
## [1] "data.frame"ただし、そのようなオブジェクトは
- どの列がジオメトリのリストカラムであるかを登録しない。
- plot メソッドを持たなくなる。
- 上記の
sfクラスの専用メソッドがすべてない。
sfc: シンプルフィーチャジオメトリリストカラム
ジオメトリを含む sf data.frame のカラムはリストであり、クラスは sfc です。 この場合のジオメトリのリストカラムは、 nc$geom または nc[[15]] によって取得することができますが、次のような方法もあります。 より一般的な方法では st_geometry を使用します。
(nc_geom <- st_geometry(nc))
## Geometry set for 100 features
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## Geodetic CRS: NAD27
## First 5 geometries:
## MULTIPOLYGON (((-81.47276 36.23436, -81.54084 3...
## MULTIPOLYGON (((-81.23989 36.36536, -81.24069 3...
## MULTIPOLYGON (((-80.45634 36.24256, -80.47639 3...
## MULTIPOLYGON (((-76.00897 36.3196, -76.01735 36...
## MULTIPOLYGON (((-77.21767 36.24098, -77.23461 3...ジオメトリは省略されている形で表示されますが、選択することで完全なジオメトリを見ることができます。例えば、最初のジオメトリは
nc_geom[[1]]
## MULTIPOLYGON (((-81.47276 36.23436, -81.54084 36.27251, -81.56198 36.27359, -81.63306 36.34069, -81.74107 36.39178, -81.69828 36.47178, -81.7028 36.51934, -81.67 36.58965, -81.3453 36.57286, -81.34754 36.53791, -81.32478 36.51368, -81.31332 36.4807, -81.26624 36.43721, -81.26284 36.40504, -81.24069 36.37942, -81.23989 36.36536, -81.26424 36.35241, -81.32899 36.3635, -81.36137 36.35316, -81.36569 36.33905, -81.35413 36.29972, -81.36745 36.2787, -81.40639 36.28505, -81.41233 36.26729, -81.43104 36.26072, -81.45289 36.23959, -81.47276 36.23436)))この表示方法は well-known text と呼ばれ、規格の一部となっています。MULTIPOLYGON という単語の後には3つの括弧がありますが、これは MULTIPOLYGON(POL1,POL2) のように複数のポリゴンで構成することができるからです。ここで POL1 は (EXT1,HOLE1,HOLE2) のように外部のリングと0個以上の内部の穴から成るかもしれません。座標のセットは括弧でくくられますので、((crds_ext)(crds_hole1)(crds_hole2)) となり、 crds_ はカンマで区切られたリングの座標のセットとなります。つまり、上の例では MULTIPOLYGON(((crds_ext))) は最初の多角形 (3) の外側の穴 (1), 穴なし (2) したがって、3つの括弧があります。
穴を持たない1つの多角形があることがわかります。
par(mar = c(0,0,1,0))
plot(nc[1], reset = FALSE) # reset = FALSE: we want to add to a plot with a legend
plot(nc[1,1], col = 'grey', add = TRUE)
しかし、このデータセットのポリゴンの中には、複数の外周リングを持つものがあり、そのようなものは
par(mar = c(0,0,1,0))
(w <- which(sapply(nc_geom, length) > 1))
## [1] 4 56 57 87 91 95
plot(nc[w,1], col = 2:7)
データ構造 MULTIPOLYGON に従って、R では行列のリストのリストのリストがあります。例えば、フィーチャ 4 のジオメトリについて、2 番目の外接リング（最初のリングは常に外接）の最初の 3 つの座標ペアを次のようにして取得します。
nc_geom[[4]][[2]][[1]][1:3,]
## [,1] [,2]
## [1,] -76.02717 36.55672
## [2,] -75.99866 36.55665
## [3,] -75.91192 36.54253ジオメトリの列は、独自のクラスを持ちます。
class(nc_geom)
## [1] "sfc_MULTIPOLYGON" "sfc"ジオメトリ・リスト・カラムのメソッドには以下のものがあります。
methods(class = 'sfc')
## [1] [ [<-
## [3] as.data.frame c
## [5] coerce format
## [7] identify initialize
## [9] Ops print
## [11] rep show
## [13] slotsFromS3 st_area
## [15] st_as_binary st_as_grob
## [17] st_as_s2 st_as_sf
## [19] st_as_text st_bbox
## [21] st_boundary st_buffer
## [23] st_cast st_centroid
## [25] st_collection_extract st_convex_hull
## [27] st_coordinates st_crop
## [29] st_crs st_crs<-
## [31] st_difference st_geometry
## [33] st_inscribed_circle st_intersection
## [35] st_intersects st_is_valid
## [37] st_is st_line_merge
## [39] st_m_range st_make_valid
## [41] st_minimum_rotated_rectangle st_nearest_points
## [43] st_node st_normalize
## [45] st_point_on_surface st_polygonize
## [47] st_precision st_reverse
## [49] st_sample st_segmentize
## [51] st_set_precision st_shift_longitude
## [53] st_simplify st_snap
## [55] st_sym_difference st_transform
## [57] st_triangulate st_union
## [59] st_voronoi st_wrap_dateline
## [61] st_write st_z_range
## [63] st_zm str
## [65] summary vec_cast.sfc
## [67] vec_ptype2.sfc
## see '?methods' for accessing help and source code座標参照系 (st_crs と st_transform) については、座標参照系 のセクションで説明します。 また、st_as_wkb と st_as_text はジオメトリのリストカラムを well-known-binary または well-known-text に変換しますが、これについては 以下 で説明します。 st_bbox は座標のバウンディングボックスを取得します。
属性は以下の通りです。
attributes(nc_geom)
## $n_empty
## [1] 0
##
## $crs
## Coordinate Reference System:
## User input: NAD27
## wkt:
## GEOGCRS["NAD27",
## DATUM["North American Datum 1927",
## ELLIPSOID["Clarke 1866",6378206.4,294.978698213898,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4267]]
##
## $class
## [1] "sfc_MULTIPOLYGON" "sfc"
##
## $precision
## [1] 0
##
## $bbox
## xmin ymin xmax ymax
## -84.32385 33.88199 -75.45698 36.58965ジオメトリタイプの混在
nc_geom のクラスは c("sfc_MULTIPOLYGON", "sfc"): sfc はすべてのジオメトリタイプで共有され、 sfc_TYPE は手元の特定のジオメトリタイプを示す TYPE と共に使用します。
2 つの “特別な” タイプがあります。GEOMETRYCOLLECTION と GEOMETRY です。GEOMETRYCOLLECTION は、各ジオメトリにさまざまなタイプのジオメトリが含まれることを示す。
(mix <- st_sfc(st_geometrycollection(list(st_point(1:2))),
st_geometrycollection(list(st_linestring(matrix(1:4,2))))))
## Geometry set for 2 features
## Geometry type: GEOMETRYCOLLECTION
## Dimension: XY
## Bounding box: xmin: 1 ymin: 2 xmax: 2 ymax: 4
## CRS: NA
## GEOMETRYCOLLECTION (POINT (1 2))
## GEOMETRYCOLLECTION (LINESTRING (1 3, 2 4))
class(mix)
## [1] "sfc_GEOMETRYCOLLECTION" "sfc"それでも、ここではジオメトリは単一のタイプです。
2つ目の GEOMETRY は、ジオメトリのリストカラムのジオメトリがさまざまなタイプであることを示しています。
(mix <- st_sfc(st_point(1:2), st_linestring(matrix(1:4,2))))
## Geometry set for 2 features
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: 1 ymin: 2 xmax: 2 ymax: 4
## CRS: NA
## POINT (1 2)
## LINESTRING (1 3, 2 4)
class(mix)
## [1] "sfc_GEOMETRY" "sfc"この 2 つは根本的に異なります。GEOMETRY はインスタンスを持たないスーパークラスで、GEOMETRYCOLLECTION はジオメトリのインスタンスです。GEOMETRY のリストカラムは、ジオメトリタイプが混在しているデータソースを読み込む際に発生します。GEOMETRYCOLLECTION は 1 つのフィーチャのジオメトリです。2 つのフィーチャのポリゴンの交点は、ポイント、ライン、ポリゴンから構成されます。
sfg: シンプルフィーチャ・ジオメトリ
シンプルフィーチャ・ジオメトリ (sfg) オブジェクトは、点、線分、多角形など、シンプルフィーチャ・ジオメトリを保持します。
シンプルフィーチャ・ジオメトリは、R のネイティブデータとして実装されており、以下のルールに従っています。
- 単一の POINT は数値ベクトルです。
- 点の集合、例えばLINESTRINGやPOLYGONのリングは
matrixで、それぞれの行に点が含まれます。 - その他の集合は
リストです。
空間データの読み書きを一括して行うことが一般的であるため、クリエイター関数は実際にはほとんど使用されません。しかし、説明のために使うには便利です。
(x <- st_point(c(1,2)))
## POINT (1 2)
str(x)
## 'XY' num [1:2] 1 2
(x <- st_point(c(1,2,3)))
## POINT Z (1 2 3)
str(x)
## 'XYZ' num [1:3] 1 2 3
(x <- st_point(c(1,2,3), "XYM"))
## POINT M (1 2 3)
str(x)
## 'XYM' num [1:3] 1 2 3
(x <- st_point(c(1,2,3,4)))
## POINT ZM (1 2 3 4)
str(x)
## 'XYZM' num [1:4] 1 2 3 4
st_zm(x, drop = TRUE, what = "ZM")
## POINT (1 2)つまり、2次元、3次元、4次元の座標を表現することができるのです。すべてのジオメトリオブジェクトは sfg (シンプルフィーチャ・ジオメトリ) を継承していますが、タイプ (例: POINT) と次元 (例: XYM) のクラス名も持っています。図に、最も一般的な7つのタイプのうち6つを示します。
ジオメトリとして単一の点を持つ POINT を除き、単一のフィーチャ（単一のレコード、または data.frame の行）に対応する残りの 6 つの一般的な単一のシンプルフィーチャ・ジオメトリタイプは、次のように作成されます。
p <- rbind(c(3.2,4), c(3,4.6), c(3.8,4.4), c(3.5,3.8), c(3.4,3.6), c(3.9,4.5))
(mp <- st_multipoint(p))
## MULTIPOINT ((3.2 4), (3 4.6), (3.8 4.4), (3.5 3.8), (3.4 3.6), (3.9 4.5))
s1 <- rbind(c(0,3),c(0,4),c(1,5),c(2,5))
(ls <- st_linestring(s1))
## LINESTRING (0 3, 0 4, 1 5, 2 5)
s2 <- rbind(c(0.2,3), c(0.2,4), c(1,4.8), c(2,4.8))
s3 <- rbind(c(0,4.4), c(0.6,5))
(mls <- st_multilinestring(list(s1,s2,s3)))
## MULTILINESTRING ((0 3, 0 4, 1 5, 2 5), (0.2 3, 0.2 4, 1 4.8, 2 4.8), (0 4.4, 0.6 5))
p1 <- rbind(c(0,0), c(1,0), c(3,2), c(2,4), c(1,4), c(0,0))
p2 <- rbind(c(1,1), c(1,2), c(2,2), c(1,1))
pol <-st_polygon(list(p1,p2))
p3 <- rbind(c(3,0), c(4,0), c(4,1), c(3,1), c(3,0))
p4 <- rbind(c(3.3,0.3), c(3.8,0.3), c(3.8,0.8), c(3.3,0.8), c(3.3,0.3))[5:1,]
p5 <- rbind(c(3,3), c(4,2), c(4,3), c(3,3))
(mpol <- st_multipolygon(list(list(p1,p2), list(p3,p4), list(p5))))
## MULTIPOLYGON (((0 0, 1 0, 3 2, 2 4, 1 4, 0 0), (1 1, 1 2, 2 2, 1 1)), ((3 0, 4 0, 4 1, 3 1, 3 0), (3.3 0.3, 3.3 0.8, 3.8 0.8, 3.8 0.3, 3.3 0.3)), ((3 3, 4 2, 4 3, 3 3)))
(gc <- st_geometrycollection(list(mp, mpol, ls)))
## GEOMETRYCOLLECTION (MULTIPOINT ((3.2 4), (3 4.6), (3.8 4.4), (3.5 3.8), (3.4 3.6), (3.9 4.5)), MULTIPOLYGON (((0 0, 1 0, 3 2, 2 4, 1 4, 0 0), (1 1, 1 2, 2 2, 1 1)), ((3 0, 4 0, 4 1, 3 1, 3 0), (3.3 0.3, 3.3 0.8, 3.8 0.8, 3.8 0.3, 3.3 0.3)), ((3 3, 4 2, 4 3, 3 3))), LINESTRING (0 3, 0 4, 1 5, 2 5))作成されているオブジェクトは
ジオメトリは空の場合もあります。
(x <- st_geometrycollection())
## GEOMETRYCOLLECTION EMPTY
length(x)
## [1] 0Well-known text, well-known binary, precision
WKT and WKB
Well-known text (WKT) と Well-known binary (WKB) は、シンプルフィーチャ形状を表す2つの符号化方式です。よく知られたテキストは、例えば、次のようなものとなります。
x <- st_linestring(matrix(10:1,5))
st_as_text(x)
## [1] "LINESTRING (10 5, 9 4, 8 3, 7 2, 6 1)"これは、人間が読むことができる形式です（ただし、先頭の ##[1] とダブルクォートは除きます）。座標は通常浮動小数点数であり、大量の情報をテキストとして移動させるのは時間がかかり、不正確です。そのため、WKB（Well-known binary）エンコーディングを使用します。
st_as_binary(x)
## [1] 01 02 00 00 00 05 00 00 00 00 00 00 00 00 00 24 40 00 00 00 00 00 00 14 40
## [26] 00 00 00 00 00 00 22 40 00 00 00 00 00 00 10 40 00 00 00 00 00 00 20 40 00
## [51] 00 00 00 00 00 08 40 00 00 00 00 00 00 1c 40 00 00 00 00 00 00 00 40 00 00
## [76] 00 00 00 00 18 40 00 00 00 00 00 00 f0 3fWKTとWKBは、どちらも以下の方法でRネイティブオブジェクトに戻すことができます。
st_as_sfc("LINESTRING(10 5, 9 4, 8 3, 7 2, 6 1)")[[1]]
## LINESTRING (10 5, 9 4, 8 3, 7 2, 6 1)
st_as_sfc(structure(list(st_as_binary(x)), class = "WKB"))[[1]]
## LINESTRING (10 5, 9 4, 8 3, 7 2, 6 1)GDAL、GEOS、空間データベース、GISなどでは、高速かつ高精度なWKBの読み書きが可能です。RネイティブオブジェクトとWKB間の変換は、コンパイル済み（C++/Rcpp）コード内のパッケージ sf によって行われ、Rのシンプルフィーチャジオメトリの入出力に再利用可能で高速なルートとなります。
精度 (Precision)
ジオメトリのリストカラム (sfc) の属性の 1 つに precision があります。これは 2 つの数値で、0 でない場合は WKB への変換時に丸められ、浮動小数点表現が原因で失敗するようなジオメトリ操作が成功する可能性があります。モデルはGEOSのもので、Java Topology Suite (JTS) からコピーしているもので、次のように動作します。
- 精度がゼロの場合(デフォルト、未指定)、何も修正されない。
- 負の値はfloat（4バイトの実数）精度に変換されます。
- 正の値は
round(x*precision)/precisionに変換されます。
精度モデルについては、こちらも参照してください。ここには、次のように書かれています。“… 小数点以下3桁の精度を指定するには、1000のスケールファクタを使用します。小数点以下3桁の精度を指定する場合(つまり1000番台に丸める場合)は、スケールファクター0.001を使用します。” と書かれています。すべての座標、つまり Z や M の値 (もしあれば) も影響を受けることに注意してください。精度`の値を選択するには、いくつかの実験が必要かもしれません。
読み書き
ここまで見てきたように、外部ファイルから空間データを読み込むには、次のような方法があります。
filename <- system.file("shape/nc.shp", package="sf")
nc <- st_read(filename)
## 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引数 quiet=TRUE を追加することで、出力を抑制することができます。また、ほぼ同じですが、より静かな
nc <- read_sf(filename)書き込みは、同様に st_write を使って:
st_write(nc, "nc.shp")
## Writing layer `nc' to data source `nc.shp' using driver `ESRI Shapefile'
## Writing 100 features with 14 fields and geometry type Multi Polygon.これを繰り返すと、ファイルがすでに存在するというエラーメッセージが表示されるので、次のようにして上書きすることができます。
st_write(nc, "nc.shp", delete_layer = TRUE)
## Deleting layer `nc' using driver `ESRI Shapefile'
## Writing layer `nc' to data source `nc.shp' using driver `ESRI Shapefile'
## Writing 100 features with 14 fields and geometry type Multi Polygon.またはその静かな代替品で、デフォルトでこれを行います。
write_sf(nc, "nc.shp") # silently overwritesドライバ固有のオプション
st_read と st_write の dsn と layer 引数は、データソース名とオプションでレイヤ名を指定することができます。 これらの正確な解釈やサポートするオプションはドライバごとに異なるので、 GDAL driver documentation を参照するのが一番よいでしょう。 例えば、データベース postgis にある PostGIS テーブルは、次のようにして読み込まれます。
meuse <- st_read("PG:dbname=postgis", "meuse")ここで、PG: という文字列は PostGIS ドライバに関するものであることを示し、その後にデータベース名、そして場合によってはポート番号とユーザ認証が続きます。引数 layer と driver が指定されていない場合、st_read はデータソースからそれらを推測して読み込むか、あるいは単に最初のレイヤを読み込んで、それ以上ある場合は警告を表示します。
また、st_readは通常、座標参照系を proj4string として読み込むが、EPSG (SRID) は読み込まない。 GDAL は proj4string の文字列から SRID (EPSG コード) を取得することができないので、必要な場合はユーザが設定しなければなりません。座標参照系](#crs)のセクションも参照してください。
st_drivers() は、利用可能なドライバとそのメタデータ (ドライバ名、書き込み可能かどうか、ラスタドライバかベクタドライバか) をリストアップしている data.frame を返す。すべてのドライバは読み込みが可能です。一般的なデータフォーマットの読み込みを以下に示します。
st_layers(dsn) はデータソース dsn に存在するレイヤをリストアップし、各レイヤのフィールド数、フィーチャ、ジオメトリタイプを提供します。
st_layers(system.file("osm/overpass.osm", package="sf"))
## Driver: OSM
## Available layers:
## layer_name geometry_type features fields crs_name
## 1 points Point NA 10 WGS 84
## 2 lines Line String NA 10 WGS 84
## 3 multilinestrings Multi Line String NA 4 WGS 84
## 4 multipolygons Multi Polygon NA 25 WGS 84
## 5 other_relations Geometry Collection NA 4 WGS 84というのは、このxmlファイルではファイル全体を読み込む必要があり、大きなファイルではコストがかかる可能性があるからです。このような場合、以下のように強制的にカウントすることができます。
Sys.setenv(OSM_USE_CUSTOM_INDEXING="NO")
st_layers(system.file("osm/overpass.osm", package="sf"), do_count = TRUE)
## Driver: OSM
## Available layers:
## layer_name geometry_type features fields crs_name
## 1 points Point 1 10 WGS 84
## 2 lines Line String 0 10 WGS 84
## 3 multilinestrings Multi Line String 0 4 WGS 84
## 4 multipolygons Multi Polygon 13 25 WGS 84
## 5 other_relations Geometry Collection 0 4 WGS 84また、kml ファイルや kmz ファイルを読み込む例もあります。
# Download .shp data
u_shp <- "http://coagisweb.cabq.gov/datadownload/biketrails.zip"
download.file(u_shp, "biketrails.zip")
unzip("biketrails.zip")
u_kmz <- "http://coagisweb.cabq.gov/datadownload/BikePaths.kmz"
download.file(u_kmz, "BikePaths.kmz")
# Read file formats
biketrails_shp <- st_read("biketrails.shp")
if(Sys.info()[1] == "Linux") # may not work if not Linux
biketrails_kmz <- st_read("BikePaths.kmz")
u_kml = "http://www.northeastraces.com/oxonraces.com/nearme/safe/6.kml"
download.file(u_kml, "bikeraces.kml")
bikraces <- st_read("bikeraces.kml")Create, read, update and delete
GDAL は、保存に対して crud (create, read, update, delete) 関数を提供します。読み込みには st_read (または read_sf) が使用されます。st_write (または write_sf) は作成に使用され、以下の引数によって更新と削除を制御します。
- このオプションは、すべてのデータベースドライバでデフォルトで
TRUEに設定されており、テーブルを追加することでデータベースを更新することができます。 - データソースが存在しない場合や、レイヤがデータソースに存在しない場合、エラーは発生しません。 データソースが存在しない場合や、データソースにレイヤが存在しない場合、エラーは発生しません。
delete_dsn=TRUEの場合、st_writeは新しく作成しているデータソースにレイヤを書き込む前に、データソースが存在すればそれを削除します。データソースが存在しない場合、エラーは発生しません。このオプションは、ディレクトリやデータベースを完全に消去してしまう可能性があるため、注意して取り扱う必要があります。
空間データベースへの接続
読み込み関数と書き込み関数である st_read() と st_write() は、空間データベースへの接続を処理し、GDAL を使用せずに直接 WKB や WKT を読み込むことができます。 DBI インターフェースを使用することを意図しているが、現在のところ、これらの関数の使用およびテストは PostGIS に限定されています。
座標参照系と変換
座標参照系（CRS）は、座標の測定単位のようなもので、特定の座標ペアが地球上のどの場所を参照しているかを指定するものです。上述しているように、sfc オブジェクト (ジオメトリのリストカラム) は CRS を格納するための 2 つの属性、epsg と proj4string を持っています。 これは、ジオメトリリストカラムに含まれるすべてのジオメトリが同じ CRS を持つ必要があることを意味します。例えば、CRS が不明な場合や、ローカルな座標系（例えば、建物内、身体、抽象的な空間など）を扱う場合には、どちらも NA にすることができます。
proj4string は、PROJ ライブラリによって理解される、CRS の一般的な文字列ベースの記述です。これは、投影の種類を定義し、（多くの場合）特定の投影のためのパラメータ値を定義し、それゆえ、無限の異なる投影をカバーすることができます。 このライブラリ(GDALも使用)は、異なるCRS間の変換や変形を行う関数を提供します。 epsgは、特定の既知のCRSの整数IDで、proj4stringに解決することができます。いくつかのproj4stringの値は、対応するepsg` ID に解決して戻すことができますが、これは常に動作するわけではありません。
proj4string 以外のデータで epsg 値を保存することの重要性は、epsg が特定のよく知られた CRS を参照しており、そのパラメータが時間とともに変化（改善）する可能性があるからです。proj4stringだけを固定すると、その改善の恩恵を受けられなくなり、データセットの出自がある程度限定されますが、再現性には貢献するかもしれません。
座標参照系の変換は st_transform を用いて行うことができます。例えば、NAD27 の経度/緯度を Web Mercator (EPSG:3857) に変換するには、以下のようにします。
nc.web_mercator <- st_transform(nc, 3857)
st_geometry(nc.web_mercator)[[4]][[2]][[1]][1:3,]
## [,1] [,2]
## [1,] -8463267 4377507
## [2,] -8460094 4377498
## [3,] -8450437 4375541Conversion, including to and from sp
変換、 sp 含む
sf オブジェクトと Spatial (パッケージ sp) から派生しているオブジェクトは、両方の方法で強制することができます。
showMethods("coerce", classes = "sf")
## Function: coerce (package methods)
## from="sf", to="Spatial"
## from="Spatial", to="sf"
methods(st_as_sf)
## [1] st_as_sf.data.frame* st_as_sf.lpp* st_as_sf.map*
## [4] st_as_sf.owin* st_as_sf.ppp* st_as_sf.ppplist*
## [7] st_as_sf.psp* st_as_sf.s2_geography* st_as_sf.sf*
## [10] st_as_sf.sfc* st_as_sf.Spatial* st_as_sf.SpatVector*
## see '?methods' for accessing help and source code
methods(st_as_sfc)
## [1] st_as_sfc.bbox* st_as_sfc.blob*
## [3] st_as_sfc.character* st_as_sfc.dimensions*
## [5] st_as_sfc.factor* st_as_sfc.list*
## [7] st_as_sfc.map* st_as_sfc.owin*
## [9] st_as_sfc.pq_geometry* st_as_sfc.psp*
## [11] st_as_sfc.raw* st_as_sfc.s2_geography*
## [13] st_as_sfc.sf* st_as_sfc.SpatialLines*
## [15] st_as_sfc.SpatialMultiPoints* st_as_sfc.SpatialPixels*
## [17] st_as_sfc.SpatialPoints* st_as_sfc.SpatialPolygons*
## [19] st_as_sfc.tess* st_as_sfc.WKB*
## see '?methods' for accessing help and source code
# anticipate that sp::CRS will expand proj4strings:
p4s <- "+proj=longlat +datum=NAD27 +no_defs +ellps=clrk66 +nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat"
st_crs(nc) <- p4s
# anticipate geometry column name changes:
names(nc)[15] = "geometry"
attr(nc, "sf_column") = "geometry"
nc.sp <- as(nc, "Spatial")
class(nc.sp)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
nc2 <- st_as_sf(nc.sp)
all.equal(nc, nc2)
## [1] "Attributes: < Component \"class\": Lengths (4, 2) differ (string compare on first 2) >"
## [2] "Attributes: < Component \"class\": 1 string mismatch >"
## [3] "Component \"geometry\": Attributes: < Component \"crs\": Component \"input\": 1 string mismatch >"
## [4] "Component \"geometry\": Attributes: < Component \"crs\": Component \"wkt\": 1 string mismatch >"Spatial* オブジェクトは MULTILINESTRING と MULTIPOLYGON のみをサポートしているので、 LINESTRING と POLYGON ジオメトリは自動的に MULTI 形式に強制されます。Spatial* から sf に変換する際に、すべてのジオメトリが1つの POLYGON (穴がある場合もある) で構成されている場合は POLYGON を、それ以外はすべて MULTIPOLYGON として返します: POLYGON と MULTIPOLYGON を混ぜたような形状は作成しません (シェイプファイルによく含まれるような形状)。forceMulti=TRUEを指定すると、これを上書きして、すべてのケースで MULTIPOLYGON を作成します。LINESの場合も状況は同じです。
幾何学的な操作
シンプルフィーチャアクセスのための標準は、いくつかの幾何学的操作を定義しています。
st_is_valid と st_is_simple は、ジオメトリが有効か単純かを示す真偽値を返します。
st_is_valid(nc[1:2,])
## [1] TRUE TRUEst_distance は、ジオメトリ間の距離を表す密な数値行列を返します。st_relate はジオメトリの各ペアに対する DE9-IM の値を持つ文字行列を返します。
x = st_transform(nc, 32119)
st_distance(x[c(1,4,22),], x[c(1, 33,55,56),])
## Units: [m]
## [,1] [,2] [,3] [,4]
## [1,] 0.00 312176.2 128338.51 475608.8
## [2,] 440548.35 114938.1 590417.79 0.0
## [3,] 18943.74 352708.6 78754.75 517511.6
st_relate(nc[1:5,], nc[1:4,])
## although coordinates are longitude/latitude, st_relate assumes that they are planar
## [,1] [,2] [,3] [,4]
## [1,] "2FFF1FFF2" "FF2F11212" "FF2FF1212" "FF2FF1212"
## [2,] "FF2F11212" "2FFF1FFF2" "FF2F11212" "FF2FF1212"
## [3,] "FF2FF1212" "FF2F11212" "2FFF1FFF2" "FF2FF1212"
## [4,] "FF2FF1212" "FF2FF1212" "FF2FF1212" "2FFF1FFF2"
## [5,] "FF2FF1212" "FF2FF1212" "FF2FF1212" "FF2FF1212"コマンド st_intersects、st_disjoint、st_touches、st_crosses、st_within、st_contains、st_overlaps、st_equals、st_covers、st_covered_by、st_equals_exact および st_is_within_distance はマッチングの結果が TRUE の疎行列か、論理上の完全行列のいずれかとなるように返します。
st_intersects(nc[1:5,], nc[1:4,])
## Sparse geometry binary predicate list of length 5, where the predicate
## was `intersects'
## 1: 1, 2
## 2: 1, 2, 3
## 3: 2, 3
## 4: 4
## 5: (empty)
st_intersects(nc[1:5,], nc[1:4,], sparse = FALSE)
## [,1] [,2] [,3] [,4]
## [1,] TRUE TRUE FALSE FALSE
## [2,] TRUE TRUE TRUE FALSE
## [3,] FALSE TRUE TRUE FALSE
## [4,] FALSE FALSE FALSE TRUE
## [5,] FALSE FALSE FALSE FALSEコマンド st_buffer、st_boundary、st_convexhull、st_union_cascaded、st_simplify、st_triangulate、st_polygonize、st_centroid、st_segmentize および st_union は新しい幾何学的構造を返します。例えば:
sel <- c(1,5,14)
geom = st_geometry(nc.web_mercator[sel,])
buf <- st_buffer(geom, dist = 30000)
plot(buf, border = 'red')
plot(geom, add = TRUE)
plot(st_buffer(geom, -5000), add = TRUE, border = 'blue')
コマンド st_intersection, st_union, st_difference, st_sym_difference は、ジオメトリのペアの関数である新しいジオメトリを返す。
par(mar = rep(0,4))
u <- st_union(nc)
plot(u)
次のコードは、2つのポリゴンの交点を計算することで、点、線、ポリゴンを含む GEOMETRYCOLLECTION を生成する方法を示しています。
opar <- par(mfrow = c(1, 2))
a <- st_polygon(list(cbind(c(0,0,7.5,7.5,0),c(0,-1,-1,0,0))))
b <- st_polygon(list(cbind(c(0,1,2,3,4,5,6,7,7,0),c(1,0,.5,0,0,0.5,-0.5,-0.5,1,1))))
plot(a, ylim = c(-1,1))
title("intersecting two polygons:")
plot(b, add = TRUE, border = 'red')
(i <- st_intersection(a,b))
## GEOMETRYCOLLECTION (POLYGON ((7 0, 7 -0.5, 6 -0.5, 5.5 0, 7 0)), LINESTRING (4 0, 3 0), POINT (1 0))
plot(a, ylim = c(-1,1))
title("GEOMETRYCOLLECTION")
plot(b, add = TRUE, border = 'red')
plot(i, add = TRUE, col = 'green', lwd = 2)
par(opar)非単純・非有効なジオメトリー
非単純な形状とは、例えば自己交差する直線（左）、非有効な形状とは、例えば細長い多角形（中）、自己交差（右）などがあります。
library(sf)
x1 <- st_linestring(cbind(c(0,1,0,1),c(0,1,1,0)))
x2 <- st_polygon(list(cbind(c(0,1,1,1,0,0),c(0,0,1,0.6,1,0))))
x3 <- st_polygon(list(cbind(c(0,1,0,1,0),c(0,1,1,0,0))))
st_is_simple(st_sfc(x1))
## [1] FALSE
st_is_valid(st_sfc(x2,x3))
## [1] FALSE FALSE