area_and_crs.md

    library(terra)
#> terra 1.7.71
set.seed(2)
query <- buffer(project(vect(geosphere::randomCoordinates(1), crs = "EPSG:4326"),   "EPSG:3857"), 
                c(1e6, 4e6), capstyle = "square")

maps::map()
plot(project(query, "EPSG:4326"), add = TRUE, border= "firebrick", lwd = 4)

expanse(query) ## correct, crs-aware
#> [1] 2.407833e+12
sf::st_area(sf::st_as_sf(query))  ## Mercator
#> 4e+12 [m^2]

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] 4e+12
set.crs(query, "EPSG:3857")  ## restore

sf::st_area(sf::st_as_sf(project(query, "EPSG:4326")))  ## correct, s2
#> 2.406033e+12 [m^2]
sf::sf_use_s2(FALSE)
#> Spherical geometry (s2) switched off
sf::st_area(sf::st_as_sf(project(query, "EPSG:4326")))  ## correct, but different - underlying ellipsoid engine

#> 2.407833e+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")
set.crs(query.cyl, NA)

expanse(query.cyl)
#> Warning: [expanse] unknown CRS. Results can be wrong
#> [1] 2.411209e+12

sf::st_area(sf::st_as_sf(query.cyl))   ## also close enough
#> [1] 2.411209e+12


## compare to local laea (or aea or many lcc or many others)
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)
}

expanse(laea_local(query)) ## correctish but might not be segmentized enough
#> [1] 2.407833e+12
expanse(laea_local(densify(query, 5000)))
#> [1] 2.411209e+12

## compare
plot(laea_local(densify(query, 5000)))

plot(laea_local(query))

^{Created on 2024-04-02 with reprex v2.0.2}