基本的な違い
\(R^2\) 上の幾何学や DE9-IM と比較して、s2 パッケージは根本的に新しい概念をいくつかもたらしています。
\(S^2\) 上の多角形は球を2つに分割する
球面 (\(S^2\)) 上では、どんな多角形も2つの領域を定義します。外側の環を追うとき、内側を定義する必要があり、その定義は囲む辺の左側です。これは、多角形を反転させれば（辺の順序を反転させれば）地球儀の他の部分を得ることができ、空の多角形（空集合）に加えて完全な多角形（地球全体）を得ることができることも意味します。
シンプルフィーチャは、リング方向にしたがっています。外側のリングは反時計回り、内側（穴）のリングは時計回りですが、外側のリングと内側のリングの違いはその位置（外側、その後に0個以上の内側）で定義されるので、ある意味でこれは時代遅れです。sf::read_sf には check_ring_dir という引数があり、リングの方向をチェックして修正することができます。間違ったリング方向でも、多くのものは動作します。
\(S^2\) の場合、リングの方向は重要です。そのため、st_as_s2 には oriented = FALSE という引数があり、すべての外側のリングが地球の半分より小さい面積を占めると仮定して、リングの方向をチェックし、修正することができます。
nc = read_sf(system.file("gpkg/nc.gpkg", package="sf")) # リングの方向が間違っている
s2_area(st_as_s2(nc, oriented = FALSE)[1:3]) # リングの方向を直し、面積が正確
## [1] 1137107793 610916077 1423145355
s2_area(st_as_s2(nc, oriented = TRUE)[1:3]) # 方向が間違っている: Earth's surface minus area
## [1] 5.100649e+14 5.100655e+14 5.100646e+14
nc = read_sf(system.file("gpkg/nc.gpkg", package="sf"), check_ring_dir = TRUE)
s2_area(st_as_s2(nc, oriented = TRUE)[1:3]) # 修正は不要
## [1] 1137107793 610916077 1423145355ここでは、地球全体を表すフルポリゴンとの差分として海を計算した例を紹介します。
g = as_s2_geography(TRUE)
g
## <s2_geography[1]>
## [1] <POLYGON ((0 -90, 0 -90))>さらに、各国を正射（オルソ画像）投影したものです。
co = s2_data_countries()
oc = s2_difference(g, s2_union_agg(co)) # oceans
b = s2_buffer_cells(as_s2_geography("POINT(-30 52)"), 9800000) # visible half
i = s2_intersection(b, oc) # visible ocean
plot(st_transform(st_as_sfc(i), "+proj=ortho +lat_0=52 +lon_0=-30"), col = 'blue')
(フルポリゴン g の WKT 表示は、シンプルフィーチャ基準に従っていないため有効な WKT ではなく、入力として使えないことに注意してください。フルポリゴンの内部表現を反映しているようで、WKTの拡張が必要なようです。）
これで、海によって覆われた地球表面の割合を計算することができます。
s2_area(oc) / s2_area(g)
## [1] 0.711301SEMI-OPEN ポリゴンの境界線
s2geometry のポリゴンは、次のようなものがあります。
- CLOSED: 境界を含み、境界上の点がポリゴンと交差している。
- OPEN: 境界を含まず、境界上の点がポリゴンと交差しない。
- SEMI-OPEN: 境界の一部を含むが、重ならないポリゴンの境界が2つ以上のポリゴンに含まれることはない。
原理的に、DE9-IM モデルは内部、境界、外部を扱い、交差述語はこれに敏感です（contains と covers の違いは全て境界に関するものです）。しかし、DE9-IM は、ポリゴンがポリゴン coverage（重なりはないが、境界は共有）を形成している場合、ポリゴンに点を一意に割り当てることができない。このため、ポリゴン単位で点数をカウントする場合、ある点がポリゴンの境界を共有している場合、その点を見逃すか（ contains ）、二重にカウントする（ covers 、 intersects ）ため、偏りが生じることがあり、後処理が必要になることがあります。SEMI-OPEN の非重複ポリゴンを使用すると、すべての点が交差点内の最大で1つのポリゴンに割り当てられることが保証されます。これは、例えば、グリッド（ラスタ）カバレッジにおいて、各グリッドセルが（通常）左上隅とその上下左右を含む場合に、どのように処理されるかに対応します。
a = as_s2_geography("POINT(0 0)")
b = as_s2_geography("POLYGON((0 0,1 0,1 1,0 1,0 0))")
s2_intersects(a, b, s2_options(model = "open"))
## [1] FALSE
s2_intersects(a, b, s2_options(model = "closed"))
## [1] TRUE
s2_intersects(a, b, s2_options(model = "semi-open")) # a toss
## [1] FALSE
s2_intersects(a, b) # default: semi-open
## [1] FALSE外接キャップ、外接矩形
sf が st_bbox() で行っているような，座標範囲に対する最小値と最大値の計算は，球面座標ではあまり意味がありません。なぜなら、球面空間では、 area covered が必ずしも座標範囲に含まれるとは限らないからです。2つの例を示します。
- 反経（経度±180）をカバーする小さな地域は、大きな経度範囲になってしまい、反経がカバーされていることが明確になりません。
- 極を含む地域は、緯度範囲が+/-90まで広がらないことになります。
S2には、外接キャップと外接矩形という2つの選択肢があります。
fiji = s2_data_countries("Fiji")
aa = s2_data_countries("Antarctica")
s2_bounds_cap(fiji)
## lng lat angle
## 1 178.7459 -17.15444 1.801369
s2_bounds_rect(c(fiji,aa))
## lng_lo lat_lo lng_hi lat_hi
## 1 177.285 -18.28799 -179.7933 -16.02088
## 2 -180.000 -90.00000 180.0000 -63.27066キャップは、外接キャップ（円）を中点（lat, lng）として、この点を中心とした角度を報告する。外接矩形は、 lat と lng 座標の _lo と _hi の境界を報告します。 フィジーの場合、 lng_lo が lng_hi よりも大きいということは、その領域が反経線を覆っている（交差している）ことを示していることに注意してください。
