correct_area.md

Here is some code to explore the variety available in calculating area across various packages.

x <- rnaturalearth::ne_countries(country = c("Greenland", "Spain"), returnclass = "sf")
as.numeric(sf::st_area(x))/1e6

## Greenland 2.166e6 km²  -41, 72
## Spain     505990 km²   -2, 40

## the numbers are different in different projections
as.numeric(sf::st_area(sf::st_transform(x, "+proj=laea")))/1e6

=## especially different in this one
as.numeric(sf::st_area(sf::st_transform(x, "EPSG:3857")))/1e6


## what did sp do? 
spx <- sf::as_Spatial(x)
## area in native units (not useful, unless you want the metrics that you see when you plot)
rgeos::gArea(spx, byid = T)

## this is the same as sf now, so sf pivots on longlat vs non-longlat 
rgeos::gArea(sf::as_Spatial(sf::st_transform(x, "EPSG:3857")), byid = TRUE)

## raster was the same, corrected for longlat 
raster::area(spx)
## just gives native for projected
raster::area(sf::as_Spatial(sf::st_transform(x, "EPSG:3857")))
## but assumes it is longlat if it looks like longlat 
spx@proj4string@projargs <- NA_character_
raster::area(spx)

## what does terra do
terra::expanse(terra::vect(x))
## gives the correct answer no matter what
terra::expanse(terra::project(terra::vect(x), "EPSG:3857"))

nocrs <- terra::project(terra::vect(x), "EPSG:3857")
terra::set.crs(nocrs, NA)
## this is the actual native area again (because terra ignores geography because it doesn't have a crs)
terra::expanse(nocrs)

nocrs_longlat <- terra::vect(x)
terra::set.crs(nocrs_longlat, NA)
terra::expanse(nocrs_longlat)  ## different to raster which assumed longlat when it looked like it

## slightly different for the longlat case, because this is ellipsoidal (same as terra)
sf::sf_use_s2(FALSE)
as.numeric(sf::st_area(x))/1e6
as.numeric(sf::st_area(sf::st_transform(x, "+proj=laea")))/1e6
as.numeric(sf::st_area(sf::st_transform(x, "EPSG:3857")))/1e6




## ok but maybe *this* laea is not suitable, they stay the same because this projection is _equal_area_
as.numeric(sf::st_area(sf::st_transform(x, "+proj=laea +lon_0=-41 +lat_0=72")))/1e6
as.numeric(sf::st_area(sf::st_transform(x, "+proj=laea +lon_0=-2 +lat_0=40")))/1e6

## if we use an equidistant projection our numbers for area are better when we're centred on the object
as.numeric(sf::st_area(sf::st_transform(x, "+proj=aeqd +lon_0=-41 +lat_0=72")))/1e6  ## good for Greenland
as.numeric(sf::st_area(sf::st_transform(x, "+proj=aeqd +lon_0=-2 +lat_0=40")))/1e6   ## good for Spain

So, the above is my current exploration of this process. I'm not totally sure how helpful this is but the thing here that I'm trying to understand is which approach is better more generally? So when not using S2 in sf, can we assume terra and sf use the same ellipsoid for calculating area when the crs is geogrpahic?

What's your take on s2 for measuring distances/areas - is it an improvement here?

thanks! this is a great way to illustrate 👌

yes to same when crs is longlat

for s2 I don't think it's an improvement, but it's probably fine for lots of scenarios, many packages in R do spherical calcs for cases when spherical accuracy is enough.

but, I need to flesh out examples a lot more to ensure i understand what's different and where. your layout is really helpful!