04-Spherical.qmd

# Spherical Geometries {#sec-spherical}

"_There are too many false conclusions drawn and stupid measurements
made when geographic software, built for projected Cartesian
coordinates in a local setting, is applied at the global scale_"
[@chrisman]

The previous chapter discussed geometries defined on the plane,
$R^2$.  This chapter discusses what changes when we consider
geometries not on the plane, but on the sphere ($S^2$).

Although we learned in @sec-cs that the shape of the Earth
is usually approximated by an ellipsoid, none of the libraries shown
in green in @fig-gdal-fig-nodetails provide access
to a comprehensive set of functions that compute on an ellipsoid.
Only the s2geometry [@R-s2; @s2geometry] library does provide it
using a sphere rather than an ellipsoid.  However, when compared
to using a flat (projected) space as we did in the previous chapter,
a sphere is a _much_ better approximation to an ellipsoid.

\index{coordinates!spherical or planar}

## Straight lines {#sec-straight}

The basic premise of _simple features_ of @sec-geometries
is that geometries are represented by sequences of points _connected
by straight lines_.  On $R^2$ (or any Cartesian space), this is
trivial, but on a sphere straight lines do not exist. The shortest
line connecting two points is an arc of the circle through both
points and the centre of the sphere, also called a _great circle
segment_. A consequence is that "the" shortest distance line
connecting two points on opposing sides of the sphere does not exist,
as any great circle segment connecting them has equal length.
Note that the GeoJSON standard [@geojson] has its own definition
of straight lines in geodetic coordinates (see Exercise 1 at the
end of this chapter).

\index{coordinates!straight line}
\index{coordinates!great circle}
\index{great circle segment}

## Ring direction and full polygon

Any polygon on the sphere divides the sphere surface in two parts
with finite area: the inside and the outside. Using the
"counter-clockwise rule" as was done for $R^2$ will not work
because the direction interpretation depends on what is defined
as inside. A convention here is to define the inside as the left
(or right) side of the polygon boundary when traversing its points
in sequence. Reversal of the node order then switches inside and
outside.

\index{polygon!inside or outside on sphere}
\index{polygon!ring direction}

In addition to empty polygons, one can define the _full polygon_
on a sphere, which comprises its entire surface. This is useful,
for instance for computing the oceans as the geometric difference
between the full polygon and the union of the land mass (see
@fig-world and @fig-srsglobe).

\index{polygon!full vs. empty}
\index{full polygon}

## Bounding box, rectangle, and cap

\index{coordinates!bounding box}
\index{bounding cap}
\index{bounding box}
\index{bounding rectangle}
\index{antimeridian}

Where in $R^2$ one can easily define bounding boxes as the range
of the $x$ and $y$ coordinates, for ellipsoidal coordinates these
ranges are not of much use when geometries cross the antimeridian
(longitude +/- 180) or one of the poles. The assumption in $R^2$
that lower $x$ values are Westwards of higher ones does not hold
when crossing the antimeridian. An alternative to delineating
an area on a sphere that is more natural is the _bounding cap_,
defined by its centre coordinates and a radius. For Antarctica,
as depicted in Figures -@fig-antarctica (a) and (c), the
bounding box formed by coordinate ranges is
```{r fig-caprect, echo=!knitr::is_latex_output()}
#| code-fold: true
#| collapse: false
#| out.width: 100%
library(sf) |> suppressPackageStartupMessages()
library(maps) |> suppressPackageStartupMessages()
library(dplyr) |> suppressPackageStartupMessages()
map(fill = TRUE, plot = FALSE) |>
  st_as_sf() |>
  filter(ID == "Antarctica") -> a
st_bbox(a)
```
which clearly does not contain the region (`ymin` being -90 and `xmax` 180).
Two geometries that do contain the region are the bounding cap:
```{r echo=!knitr::is_latex_output()}
#| code-fold: true
#| collapse: false
library(s2)
s2_bounds_cap(a)
```
and the bounding _rectangle_:
```{r echo=!knitr::is_latex_output()}
#| code-fold: true
#| collapse: false
s2_bounds_rect(a)
```
For an area spanning the antimeridian, here the Fiji island country,
the bounding box:
```{r echo=!knitr::is_latex_output()}
#| code-fold: true
#| collapse: false
map(fill = TRUE, plot = FALSE) |>
  st_as_sf() |>
  filter(ID == "Fiji") -> Fiji
st_bbox(Fiji)
```
seems to span most of the Earth, as opposed to the bounding rectangle:
```{r echo=!knitr::is_latex_output()}
#| code-fold: true
#| collapse: false
s2_bounds_rect(Fiji)
```
where a value `lng_lo` _larger_ than `lng_hi` indicates that the
bounding rectangle spans the antimeridian. This property could not
be inferred from the coordinate ranges.

## Validity on the sphere

\index{geometry!valid on the sphere}
\index{polygon!valid on the sphere}

Many global datasets are given in ellipsoidal coordinates but are
prepared in a way that they "work" when interpreted on the $R^2$
space [-180,180] $\times$ [-90,90]. This means that:

* geometries crossing the antimeridian (longitude +/- 180) are cut in
half, such that they no longer cross it (but nearly touch each other)
* geometries including a pole, like Antarctica, are cut at +/- 180 and
make an excursion through -180,-90 and 180,-90 (both representing the
Geographic South Pole)

@fig-antarctica  shows two different representations of
Antarctica, plotted with ellipsoidal coordinates taken as $R^2$
(top) and in a Polar Stereographic projection (bottom), without
(left) and with (right) an excursion through the Geographic South
Pole. In the projections as plotted, polygons (b) and
(c) are valid, polygon (a) is not valid as it self-intersects, and
polygon (d) is not valid because it traverses the same edge to the
South Pole twice.  On the sphere ($S^2$), polygon (a) is valid but
(b) is not, for the same reason as (d) is not valid.

```{r fig-antarctica, echo=!knitr::is_latex_output()}
#| code-fold: true
#| out.height: 70%
#| fig.cap: "Antarctica polygon, (a, c): _not_ passing through `POINT(-180 -90)`; (b, d): passing through `POINT(-180 -90)` and `POINT(180 -90)`"
# maps:
par(mfrow = c(2,2))
par(mar = c(1,1.2,1,1))
m <- st_as_sf(map(fill=TRUE, plot=FALSE))
m <- m[m$ID == "Antarctica", ]
plot(st_geometry(m), asp = 2)
title("a (not valid)")
# ne:
library(rnaturalearth)
ne <- ne_countries(returnclass = "sf")
ne <- ne[ne$region_un == "Antarctica", "region_un"]
plot(st_geometry(ne), asp = 2)
title("b (valid)")
# 3031
m |>
  st_geometry() |>
  st_transform(3031) |>
  plot()
title("c (valid)")
ne |>
  st_geometry() |>
  st_transform(3031) |>
  plot()
title("d (not valid)")
```

## Exercises

For the following exercises, use R where possible or relevant.

1. How does the [GeoJSON](https://tools.ietf.org/html/rfc7946) format [@geojson] define
   "straight" lines between ellipsoidal coordinates (Section 3.1.1)?
   Using this definition of straight, how does `LINESTRING(0 85,180 85)` 
   look like in an Arctic polar projection? How could this geometry be modified to
   have it cross the North Pole?
2. For a typical polygon on $S^2$, how can you find out ring direction?
3. Are there advantages of using bounding caps over using bounding boxes? If so, list them.
4. Why is, for small areas, the orthographic projection centred
   at the area a good approximation of the geometry as handled on $S^2$?
5. For `rnaturalearth::ne_countries(country = "Fiji", returnclass = "sf")`,
   check whether the geometry is valid on $R^2$,
   on an orthographic projection centred on the country, and on $S^2$.
   How can the geometry be made valid on $S^2$? Plot the resulting
   geometry back on $R^2$.  Compare the centroid of the country,
   as computed on $R^2$ and on $S^2$, and the distance between the two.
6. Consider dataset `gisco_countries` in R package **giscoR**, and select
   the country with `NAME_ENGL == "Fiji"`. Does it have a valid geometry
   on the sphere? If so, how was this accomplished?