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

suppressPackageStartupMessages({
    #' @title get administritive outlines for a country
    #' @description using the geoBoundaires API,
    #' download the geojson outline of a country
    #' @param country character vector of country names
    #' @param admin_level character vector of admin levels to download
    #' @return sf object of the outlines
    #' @details check out the documentation for the geoboundaries API at:
    #' geoBoundaries.org
    #'
    geo_bounds <- function(country, admin_level = c("ADM0", "ADM1", "ADM2")) {
        country <- countrycode::countrycode(country,
            origin = "country.name",
            destination = "iso3c"
        )
        url <- paste(
            "https://www.geoboundaries.org/api/current/gbOpen/",
            country, admin_level[1],
            sep = "/"
        )
        get <- httr::GET(url)
        cont <- httr::content(get, as = "parsed")
        area <- sf::read_sf(cont$gjDownloadURL)
        return(area)
    }

    #' @title get the area of a polygon in multiple projections and from packages
    #'
    #' @param v sf object
    #' @param add_crs character or numeric vector of crs to add to the sf object
    test_proj_area <- function(v, add_crs = c(4326, 3857)) {
        crs_opts <- suppressMessages(crsuggest::suggest_crs(v))
        auto_epsg <- as.numeric(crs_opts$crs_code)

        area_tab <- function(x) {
            if (!suppressWarnings(is.na((as.numeric(x))))) {
                c_crs <- paste0("EPSG:", x)
            } else {
                c_crs <- x
            }

            vp <- sf::st_transform(v, c_crs)

            sf_a <- function(.s2, .crs) {
                suppressMessages(sf::sf_use_s2(.s2))
                vp |>
                    sf::st_transform(.crs) |>
                    sf::st_area() |>
                    units::set_units("km^2")
            }

            terra_a <- function(.transform, .crs) {
                tv <- terra::vect(vp)
                # tv <- terra::project(tv, .crs)
                terra::expanse(tv, transform = .transform) |>
                    units::as_units("m^2") |>
                    units::set_units("km^2")
            }

            tibble::tibble(
                "crs" = c_crs,
                "lon-lat" = sf::st_is_longlat(vp),
                "sf s2" = sf_a(TRUE, c_crs),
                "sf no s2" = sf_a(FALSE, c_crs),
                "terra trans" = terra_a(TRUE, c_crs),
                "terra no trans" = terra_a(FALSE, c_crs),
            )
        }

        purrr::map(c(add_crs, auto_epsg), area_tab) |>
            dplyr::bind_rows()
    }

    bolivia <- geo_bounds("Bolivia")

    plot(sf::st_geometry(bolivia), axes = TRUE)

    tab <- test_proj_area(bolivia,
        add_crs = c(
            4326,
            4674,
            9005,
            3857,
            "+proj=laea +lon_0=-41 +lat_0=72", # inappropriate laea for bolivia
            "+proj=laea +lon_0=-64 +lat_0=-16", # appropriate laea for bolivia
            "+proj=aeqd +lon_0=-2 +lat_0=40", # inappropriate aeqd for bolivia
            "+proj=aeqd +lon_0=-64 +lat_0=-16" # appropriate aeqd for bolivia
        )
    )

    options(width = 100)
    print(dplyr::mutate(
        tab,
        across(where(is.numeric), ~ tibble::num(., digits = 3))
    ))
})

#> # A tibble: 18 × 6
#>    crs                              `lon-lat`     `sf s2`  `sf no s2` `terra trans` `terra no trans`
#>    <chr>                            <lgl>       <num:.3!>   <num:.3!>     <num:.3!>        <num:.3!>
#>  1 EPSG:4326                        TRUE      1072375.122 1068784.978   1068784.978      1068784.978
#>  2 EPSG:4674                        TRUE      1072375.122 1068784.978   1068784.978      1068784.978
#>  3 EPSG:9005                        TRUE      1072375.122 1068784.978   1068784.978      1068784.978
#>  4 EPSG:3857                        FALSE     1174659.446 1174659.446   1068784.978      1174659.446
#>  5 +proj=laea +lon_0=-41 +lat_0=72  FALSE     1068779.063 1068779.063   1068784.978      1068779.063
#>  6 +proj=laea +lon_0=-64 +lat_0=-16 FALSE     1068784.760 1068784.760   1068784.978      1068784.760
#>  7 +proj=aeqd +lon_0=-2 +lat_0=40   FALSE     1534126.450 1534126.450   1068784.978      1534126.450
#>  8 +proj=aeqd +lon_0=-64 +lat_0=-16 FALSE     1069692.568 1069692.568   1068784.978      1069692.568
#>  9 EPSG:5355                        FALSE     1071080.602 1071080.602   1068784.978      1071080.602
#> 10 EPSG:24880                       FALSE     1071106.590 1071106.590   1068785.100      1071106.590
#> 11 EPSG:5356                        FALSE     1076083.973 1076083.973   1068784.978      1076083.973
#> 12 EPSG:5357                        FALSE     1088049.640 1088049.640   1068784.978      1088049.640
#> 13 EPSG:24881                       FALSE     1088071.455 1088071.455   1068785.002      1088071.455
#> 14 EPSG:5389                        FALSE     1076083.973 1076083.973   1068784.978      1076083.973
#> 15 EPSG:24893                       FALSE     1080647.120 1080647.120   1068785.803      1080647.120
#> 16 EPSG:5876                        FALSE     1076086.453 1076086.453   1068784.978      1076086.453
#> 17 EPSG:5536                        FALSE     1088110.007 1088110.007   1068784.978      1088110.007
#> 18 EPSG:22521                       FALSE     1088050.048 1088050.048   1068784.978      1088050.048

^{Created on 2023-07-06 with reprex v2.0.2}

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!