## helper functions for later 

## define a local projection from a terra object
laea_local <- function(x) {
  centre <- terra::geom(terra::project(terra::centroids(x), "EPSG:4326"))[1, c("x", "y")]
  prj <- sprintf("+proj=laea +lon_0=%f +lat_0=%f", centre[1], centre[2])
  terra::project(x, prj)
}
## make a polygon ring from an extent (xmin,xmax,ymin,ymax)
make_poly <- function(extent) {
  matrix(ex[c(1, 1, 2, 2, 1, 
              4, 3, 3, 4, 4)], ncol = 2)
}
library(terra)

maps::map()

ex <- c(-120, -100, 35, 55)

rect(ex[1], ex[3], ex[2], ex[4], border = "firebrick", lwd = 3)

What is the area

What is the area of this box? There are a few numbers.

Value Method
3.486716e+12 Use map projection correctly.
3.469788e+12 Use ellipsoidal methods (hardcore mode, GeographicLib)
3.462492e+12 Use spherical method (hardcore mode, S2)
7.119961e+12 Use a specific map projection incorrectly (Mercator)

GeoBox

GeoBox is a really cool helper for geographic extent in Python, it came out of OpenDataCube.

We are specifying a raster abstractly, this is four numbers for extent, two numbers for dimension (resolution), and a string for CRS. When we ask for area we have to specifically specify the projection to use (the authors know this scene well). We are subject to the resolution of the grid and will get less accurate results for coarser grids, even though we’re only using the extent at the start, the extent gets materalized into a geospatial polygon in the nominated CRS, with the segmentizing around the boundary depending on the grid itself.

library(reticulate)
GeoBox <- import("odc.geo.geobox")$GeoBox

g <- GeoBox$from_bbox(ex[c(1, 3, 2, 4)], crs = "EPSG:4326", shape = c(1024, 1024))
g$footprint("+proj=laea +lon_0=-110 +lat_0=45")$area
## [1] 3.486716e+12

GeographicLib and S2

Ellipsoidal and spherical (hardcore mode).

geosphere::areaPolygon(poly <- make_poly(ex))
## [1] 3.469788e+12
s2_poly <- s2::s2_make_polygon(poly[,1], poly[,2])
s2::s2_area(s2_poly)
## [1] 3.462492e+12

terra and sf

query <- project(as.polygons(ext(ex), crs = "EPSG:4326"), "EPSG:3857")
expanse(query) ## ellipsoidal
## [1] 3.469788e+12
sf::st_area(sf::st_as_sf(query))  ## Mercator
## 7.119961e+12 [m^2]
query <- set.crs(query, NA)  ## beware doing this or to a copy, it's "in-place"
expanse(query)  ## also wrong
## Warning: [expanse] unknown CRS. Results can be wrong
## [1] 7.119961e+12
query <- set.crs(query, "EPSG:3857")  ## restore

suppressMessages(sf::sf_use_s2(TRUE))
sf::st_area(sf::st_as_sf(project(query, "EPSG:4326")))  ## correct, s2
## 3.462492e+12 [m^2]
suppressMessages(sf::sf_use_s2(FALSE))
sf::st_area(sf::st_as_sf(project(query, "EPSG:4326")))  ## correct, but different - underlying ellipsoid engine
## 3.469788e+12 [m^2]
## now compare non-crs-aware, but a good projection for this isomorphic case (Merc -> longlat -> cea)
query.cyl <- project(query, "+proj=cea")
query.cyl <- set.crs(query.cyl, NA)
expanse(query.cyl)  ## closest to the correct map projection way
## Warning: [expanse] unknown CRS. Results can be wrong
## [1] 3.486722e+12
sf::st_area(sf::st_as_sf(query.cyl))   ## also correct unsurprising
## [1] 3.486722e+12
expanse(laea_local(query)) ## correctish but might not be segmentized enough
## [1] 3.469788e+12
expanse(laea_local(densify(query, 5000)))
## [1] 3.486722e+12
## compare
plot(laea_local(densify(query, 5000)))
plot(laea_local(query), add = TRUE, lty = 2)