Skip to content

Instantly share code, notes, and snippets.

@mdsumner
Last active July 6, 2023 20:43
Show Gist options
  • Save mdsumner/6caeefaeb5200ddf5e77d2aeced3bebd to your computer and use it in GitHub Desktop.
Save mdsumner/6caeefaeb5200ddf5e77d2aeced3bebd to your computer and use it in GitHub Desktop.

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
@mdsumner
Copy link
Author

mdsumner commented Jul 6, 2023

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!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment