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

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)

# 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

• Equal Area projections can't keep the area perfectly equal at all locations and for all areas, especially at continental scale. – Vince May 10 at 11:51
• @Vince got it! If you can point me to a good source illustrating that issue, I'll accept it as the answer to my question – HAVB May 10 at 12:58
• You might be able to reduce the differences a bit (although the first one is only 0.004% different) by densifying the polygons. – mkennedy May 10 at 17:23
• @mkennedy does that mean that if I had "perfect" polygons, ie with infinite density thus needing no interpolation, the equal area projection would estimate the same areas no matter where the center of projection is set? – HAVB May 10 at 19:36
• I double-checked with our map projections expert. Yes, that is true but because we calculate area (and distance) using discrete points...you get differences. – mkennedy May 10 at 19:57

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!