Same equal area projection results in different areas when shifting the center of projection

Question

I found that using an equal area projection results in different areas for the same polygons if you shift the center of projection (ie to minimize the distortion of your area of interest when you project it).

How is it possible? This defies my understanding of an equal area projection keeping areas intact.

Example:

library(sf)
options(scipen = 20)

south_america <- st_read("https://raw.githubusercontent.com/codeforamerica/click_that_hood/master/public/data/south-america.geojson")


# Transform to Lambert azimuthal equal-area projection, default center (0° N 0° E)
south_america_laea <- st_transform(south_america, "+proj=laea +nodefs")   
plot(st_geometry(south_america_laea))

st_area(south_america_laea)

Units: [m^2]

83375755616 8467140414780 737655036685 2782483050075 1086627258107 1134458903620 254527695525

211024005330 1290996983204 399946127945 144935413395 177467272516 910435563567 11429672354

# Same Lambert azimuthal equal-area projection, just shift the center of projection to center it on the continent 

south_america_laea_2 <- st_transform(south_america, "+proj=laea +lat_0=-32 +lon_0=-60")

plot(st_geometry(south_america_laea_2))

st_area(south_america_laea_2)

Units: [m^2]

83372539973 8467111444577 737677003441 2782466797426 1086627986050 1134453525235 254528299278

211021681369 1291025934349 399948128115 144931686414 177467118512 910429865771 11429872564

Now the areas are different! Not even a linear transformation, just slightly different numbers

Answer 1

Densifying the geometry - adding more vertices - makes the area estimates much more similar. The difference in area between equal-area projections with a different centre varies with the square of the number of vertices.

Create a large lat-lon WGS84 square - in South America with only 4 vertices; the difference in areas is one in 7,366:

sq <- data.frame(id = 1, wkt = "Polygon ((-64 -9, -59 -9, -59 -14, -64 -14, -64 -9))")
sq <- st_as_sf(sq, wkt = 'wkt', crs = 4326)

sq.laea <- st_transform(sq, "+proj=laea +nodefs")
sq.laea2 <- st_transform(sq, "+proj=laea +lat_0=-32 +lon_0=-60")

a1 <- st_area(sq.laea)
a2 <- st_area(sq.laea2)

1 / (as.numeric(a2/a1) - 1) # 7366

With 16 vertices the difference drops to 1 in 117,740:

sq2 <- data.frame(id = 1, wkt = "Polygon ((-64 -9, -62.75 -9, -61.5 -9, -60.25 -9, -59 -9, -59 -10.25, -59 -11.5, -59 -12.75, -59 -14, -60.25 -14, -61.5 -14, -62.75 -14, -64 -14, -64 -12.75, -64 -11.5, -64 -10.25, -64 -9))")
sq2 <- st_as_sf(sq2, wkt = 'wkt', crs = 4326)

sq2.laea <- st_transform(sq2, "+proj=laea +nodefs")
sq2.laea2 <- st_transform(sq2, "+proj=laea +lat_0=-32 +lon_0=-60")

a1 <- st_area(sq2.laea)
a2 <- st_area(sq2.laea2)

1 / (as.numeric(a2/a1) - 1) # 117740

With 64 vertices the difference drops to 1 in 1,883,725:

sq3 <- data.frame(id = 1, wkt = "Polygon ((-64 -9, -63.6875 -9, -63.375 -9, -63.0625 -9, -62.75 -9, -62.4375 -9, -62.125 -9, -61.8125 -9, -61.5 -9, -61.1875 -9, -60.875 -9, -60.5625 -9, -60.25 -9, -59.9375 -9, -59.625 -9, -59.3125 -9, -59 -9, -59 -9.3125, -59 -9.625, -59 -9.9375, -59 -10.25, -59 -10.5625, -59 -10.875, -59 -11.1875, -59 -11.5, -59 -11.8125, -59 -12.125, -59 -12.4375, -59 -12.75, -59 -13.0625, -59 -13.375, -59 -13.6875, -59 -14, -59.3125 -14, -59.625 -14, -59.9375 -14, -60.25 -14, -60.5625 -14, -60.875 -14, -61.1875 -14, -61.5 -14, -61.8125 -14, -62.125 -14, -62.4375 -14, -62.75 -14, -63.0625 -14, -63.375 -14, -63.6875 -14, -64 -14, -64 -13.6875, -64 -13.375, -64 -13.0625, -64 -12.75, -64 -12.4375, -64 -12.125, -64 -11.8125, -64 -11.5, -64 -11.1875, -64 -10.875, -64 -10.5625, -64 -10.25, -64 -9.9375, -64 -9.625, -64 -9.3125, -64 -9))")
sq3 <- st_as_sf(sq3, wkt = 'wkt', crs = 4326)

sq3.laea <- st_transform(sq3, "+proj=laea +nodefs")
sq3.laea2 <- st_transform(sq3, "+proj=laea +lat_0=-32 +lon_0=-60")

a1 <- st_area(sq3.laea)
a2 <- st_area(sq3.laea2)

1 / (as.numeric(a2/a1) - 1) # 1883725

For the WKT geometries I used QGIS to densify the geometry of the basic square (copying and pasting WKT data) - couldn't find a 'densify' function anywhere in R packages!