Fundamental differences
Compared to geometry on \(R^2\), and
DE9-IM, the s2 package brings a few fundamentally new
concepts, which are discussed first.
Polygons on \(S^2\) divide the sphere in two parts
On the sphere (\(S^2\)), any polygon defines two areas; when following the exterior ring, we need to define what is inside, and the definition is the left side of the enclosing edges. This also means that we can flip a polygon (by inverting the edge order) to obtain the other part of the globe, and that in addition to an empty polygon (the empty set) we can have the full polygon (the entire globe).
Simple feature geometries should obey a ring direction too: exterior
rings should be counter clockwise, interior (hole) rings should be
clockwise, but in some sense this is obsolete as the difference between
exterior ring and interior rings is defined by their position (exterior,
followed by zero or more interior). sf::read_sf() has an
argument check_ring_dir that checks, and corrects, ring
directions and many (legacy) datasets have wrong ring directions. With
wrong ring directions, many things still work.
For \(S^2\), ring direction is
essential. For that reason, st_as_s2 has an argument
oriented = FALSE, which will check and correct ring
directions, assuming that all exterior rings occupy an area smaller than
half the globe:
nc = read_sf(system.file("gpkg/nc.gpkg", package="sf")) # wrong ring directions
s2_area(st_as_s2(nc, oriented = FALSE)[1:3]) # corrects ring direction, correct area:
## [1] 1137107793 610916077 1423145355
s2_area(st_as_s2(nc, oriented = TRUE)[1:3]) # wrong direction: 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]) # no second correction needed here:
## [1] 1137107793 610916077 1423145355The default conversion from sf to s2 uses
oriented = FALSE, so that we get
Here is an example where the oceans are computed as the difference from the full polygon representing the entire globe,
g = st_as_sfc("POLYGON FULL", crs = 'EPSG:4326')
g
## Geometry set for 1 feature
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 90
## Geodetic CRS: WGS 84
## POLYGON FULLand the countries, and shown in an orthographic projection:
options(s2_oriented = TRUE) # don't change orientation from here on
co = st_as_sf(s2_data_countries())
oc = st_difference(g, st_union(co)) # oceans
b = st_buffer(st_as_sfc("POINT(-30 52)", crs = 'EPSG:4326'), 9800000) # visible half
i = st_intersection(b, oc) # visible ocean
plot(st_transform(i, "+proj=ortho +lat_0=52 +lon_0=-30"), col = 'blue')
(Note that the printing of POLYGON FULL is not valid WKT
according to the simple feature standard, which does not include
this.)
We can now calculate the proportion of the Earth’s surface covered by oceans:
Semi-open polygon boundaries
Polygons in s2geometry can be
- CLOSED: they contain their boundaries, and a point on the boundary intersects with the polygon
- OPEN: they do not contain their boundaries, points on the boundary do not intersect with the polygon
- SEMI-OPEN: they contain part of their boundaries, but no boundary of non-overlapping polygons is contained by more than one polygon.
In principle the DE9-IM model deals with interior, boundary and exterior, and intersection predicates are sensitive to this (the difference between contains and covers is all about boundaries). DE9-IM however cannot uniquely assign points to polygons when polygons form a polygon coverage (no overlaps, but shared boundaries). This means that if we would count points by polygon, and some points fall on shared polygon boundaries, we either miss them (contains) or we count them double (covers, intersects); this might lead to bias and require post-processing. Using SEMI-OPEN non-overlapping polygons guarantees that every point is assigned to maximally one polygon in an intersection. This corresponds to e.g. how this would be handled in a grid (raster) coverage, where every grid cell (typically) only contains its upper-left corner and its upper and left sides.
a = st_as_sfc("POINT(0 0)", crs = 'EPSG:4326')
b = st_as_sfc("POLYGON((0 0,1 0,1 1,0 1,0 0))", crs = 'EPSG:4326')
st_intersects(a, b, model = "open")
## Sparse geometry binary predicate list of length 1, where the predicate
## was `intersects'
## 1: (empty)
st_intersects(a, b, model = "closed")
## Sparse geometry binary predicate list of length 1, where the predicate
## was `intersects'
## 1: 1
st_intersects(a, b, model = "semi-open") # a toss
## Sparse geometry binary predicate list of length 1, where the predicate
## was `intersects'
## 1: (empty)
st_intersects(a, b) # default: closed
## Sparse geometry binary predicate list of length 1, where the predicate
## was `intersects'
## 1: 1Bounding cap, bounding rectangle
Computing the minimum and maximum values over coordinate ranges, as
sf does with st_bbox(), is of limited value
for spherical coordinates because due the the spherical space, the
area covered is not necessarily covered by the coordinate
range. Two examples:
- small regions covering the antimeridian (longitude +/- 180) end up with a huge longitude range, which doesn’t make clear the antimeridian is spanned
- regions including a pole will end up with a latitude range not extending to +/- 90
S2 has two alternatives: the bounding cap and the bounding rectangle:
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.27066The cap reports a bounding cap (circle) as a mid point (lat, lng) and
an angle around this point. The bounding rectangle reports the
_lo and _hi bounds of lat and
lng coordinates. Note that for Fiji, lng_lo
being higher than lng_hi indicates that the region covers
(crosses) the antimeridian.
