3

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

enter image description here

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

enter image description here

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

  • 2
    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
  • 2
    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
  • 1
    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
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!

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.