st_transform appears to change the area of the geometry of a shapefile

Question

Shouldn't the area of polygons remain the same when you transform the shapefile projection? I am currently working with the shapefile here (https://data.humdata.org/dataset/f814a950-4d4e-4f46-a880-4da5522f14c4/resource/731e11cb-be02-46cf-8347-db0a07abff4e/download/gin_admbnda_adm4_2021_ocha.zip) for the country GIN which I call gin_shp after I read it in with the st_read function in R as follows:

gin_shp <- st_read(dsn = "INPUT FOLDER", layer = "sous_prefectures") 
gin_shp ##this shows the crs is WGS84 obviously in degrees

##compute the area of the multipolygon geometries
gin_shp$area <- sf::st_area(gin_shp) ##computes areas in m^2
gin_shp$area <- units::set_units(gin_shp$area, "km^2") ##converting the area to square kms

sum(gin_shp$area) ##compute total area of the entire country (google seems to agree!)
245084.2 [km^2]

I need to create 1sqkm grids (using sf::st_make_grid) so I tried to transform crs projection from degrees to UTM by calling the st_transform function as follows:

crs_dt <- rgdal::make_EPSG() ##first load the dataframe of CRS projections available in R

gin_shp <- sf::st_transform(gin_shp, crs = crs_dt$prj4[crs_dt$code == 4328]) ##select one and assign it to gin_shp

#now I try to compute total area under new projection regime by running the same code as before
gin_shp$area <- sf::st_area(gin_shp) ##computes areas in m^2
gin_shp$area <- units::set_units(gin_shp$area, "km^2") #converting the area to square kms

#compute total area of the entire country 
sum(gin_shp$area) 
44363.83 [km^2] (way way off)

Do you know why this might be happening? Any ideas how to fix it?

Answer 1

No, the area of polygons is not guaranteed to remain constant when transforming between different CRSes. The Earth is round, maps & computer screens are flat - something has to give; either area or shape has to be distorted somewhat.

There are some area preserving projections - such as the Mollweide - but these are more of exception than rule.

For an exaggerated example consider the world dataset, taken from {giscoR} package. Greenland (close to the pole) and Congo (on equator) have roughly the same area on sphere (calculated on WGS84 using spherical geometry tools) but wildly different one when projected. Especially when projected to Mercator and its derivatives (e.g. Web Mercator as used by Google maps).

library(sf)
library(dplyr)
library(giscoR)
library(ggplot2)

# the world in, 1: 20M
svet <- gisco_get_countries(resolution = "20")

# Greenland & Congo - to be drawn in red
glmd <- svet %>% 
   filter(CNTR_ID %in% c('GL', 'CD'))

# web mercator = default for google maps and other web based tools
ggplot() +
   geom_sf(data = glmd, fill = "red", color = NA) +
   geom_sf(data = svet, fill = NA, color = "gray45") +
   coord_sf(crs = st_crs("EPSG:3857"),
            ylim = c(-20e6, 20e6))

# Mollweide - equal area projection
ggplot() +
   geom_sf(data = glmd, fill = "red", color = NA) +
   geom_sf(data = svet, fill = NA, color = "gray45") +
   coord_sf(crs = st_crs("ESRI:53009"))