This vignette describes what spherical geometry implies, and how package sf
uses the s2geometry library (https://s2geometry.io) for geometrical measures, predicates and transformations.
Spatial coordinates either refer to projected (or Cartesian) coordinates, meaning that they are associated to points on a flat space, or to unprojected or geographic coordinates, when they refer to angles (latitude, longitude) pointing to locations on a sphere (or ellipsoid). The flat space is also referred to as R2, the sphere as S2.
Package sf
implements simple features, a standard for point, line, and polygon geometries where geometries are built from points (nodes) connected by straight lines (edges). The simple feature standard does not say much about its suitability for dealing with geographic coordinates, but the topological relational system it builds upon (DE9-IM) refer to R2, the two-dimensional flat space.
Yet, more and more data are routinely served or exchanged using geographic coordinates. Using software that assumes an R2, flat space may work for some problems, and although sf
up to version 0.9-x had some functions in place for spherical/ellipsoidal computations (from package lwgeom
, for computing area, length, distance, and for segmentizing), it has also happily warned the user that it is doing R2, flat computations with such coordinates with messages like
although coordinates are longitude/latitude, st_intersects assumes that they are planar
hinting to the responsibility of the user to take care of potential problems. Doing this however leaves ambiguities, e.g. whether LINESTRING(-179 0,179 0)
POINT(0 0)
, orPOINT(180 0)
and whether it is * a straight line, cutting through the Earth's surface, or * a curved line following the Earth's surface
Starting with sf
version 1.0, sf
uses the new package s2
(Dunnington, Pebesma, Rubak 2020) for spherical geometry, which has functions for computing pretty much all measures, predicates and transformations on the sphere. This means:
The s2
package is really a wrapper around the C++ s2geometry library which was written by Google, and which is used in many of its products (e.g. Google Maps, Google Earth Engine, Bigquery GIS) and has been translated in several programming other languages.
Compared to geometry on R2, and DE9-IM, the s2
package brings a few fundamentally new concepts, which are discussed first.
On the sphere (S2), 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 S2, 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:
library(sf)
nc = read_sf(system.file("gpkg/nc.gpkg", package="sf")) # wrong ring directions
library(s2)
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 1423145355
Here is an example where the oceans are computed as the difference from the full polygon,
as_s2_geography(TRUE)
## <s2_geography[1]>
## [1] <POLYGON ((0 -90, 0 -90))>
and the countries, and shown in an orthographic projection:
co = s2_data_countries()
oc = s2_difference(as_s2_geography(TRUE), 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')
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 mostly common 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); this leads to bias or need for 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 = 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
Computing the minimum and maximum values over coordinate ranges, as sf
does with st_bbox()
, is of limited value for spherical coordinates because
S2 has two alternatives: the cap and the lat_lng_rect
:
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
The cap reports a bounding cap (circle) as a mid point (lat, lng) and an angle around this point. The rect reports the _lo
and _hi
bounds of lat
and lng
, as well as its center (_cnt
) values. Note that for Fiji, lng_lo
being higher than lng_hi
indicates that the region covers (crosses) the antimeridian; the lng_cnt
values is not the mean of lng_lo
and lng_hi
.
The two-dimensional R2 library that was formerly used by sf
is GEOS, and sf
can be instrumented to use GEOS or sf
. First we will ask if s2
is being used by default:
sf_use_s2()
## [1] TRUE
then we can switch it of (and use GEOS) by
sf_use_s2(FALSE)
and switch it on (and use S2) by
sf_use_s2(TRUE)
library(sf)
library(units)
## udunits system database from /usr/share/xml/udunits
nc = read_sf(system.file("gpkg/nc.gpkg", package="sf"))
sf_use_s2(TRUE)
a1 = st_area(nc)
sf_use_s2(FALSE)
a2 = st_area(nc)
plot(a1, a2)
abline(0, 1)
summary((a1 - a2)/a1)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.638e-04 -1.650e-04 -7.133e-05 -6.448e-05 1.598e-05 2.817e-04
nc_ls = st_cast(nc, "MULTILINESTRING")
l1 = st_length(nc_ls)
l2 = st_length(nc_ls, use_lwgeom = TRUE)
plot(l1 , l2)
abline(0, 1)
summary((l1-l2)/l1)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0
d1 = st_distance(nc, nc[1:10,])
d2 = st_distance(nc, nc[1:10,], use_lwgeom = TRUE)
dim(d1)
## [1] 100 10
dim(d2)
## [1] 100 10
plot(as.vector(d1), as.vector(d2))
abline(0, 1)
summary(as.vector(d1)-as.vector(d2))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0