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
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!