Incorrect result when I reproject a SpatialPolygonsDataFrame using spTransform()

Question

R version 4.4.0 (2024-04-24 ucrt) -- "Puppy Cup" Platform: x86_64-w64-mingw32/x64

I have a very complex SpatialPolygonsDataFrame representing a Natura 2000 site from Sweden, SE0820430, shown in black:

First, I checked the area for the original projection (LAEA), which was 1753394516 square meters. The registered area, according to this BISE page, is 1760923000 square meters, which is relatively similar (99.57% of the registered area).

The problem is that when I run the code below to transform the projection, I see that a new polygon/shape has been created which has nothing in common with the original SpatialPolygonsDataFrame (see figure at the end). Not surprisingly, the area is very different (around 13 square meters). I do not receive either a warning or an error.

rm(list = ls()) 
library(raster)
library(sp)

#We read the spatial information of the Nature2000:
spatial_nature2000_original <- shapefile("Natura2000_end2021_rev1_epsg3035.shp")

#Subset of a feature/protected area which I am using as example, the site "SE0820430"
spatial_nature2000_original_SE0820430 <- subset(spatial_nature2000_original, SITECODE == "SE0820430")
save(spatial_nature2000_original_SE0820430, file = 'spatial_nature2000_original_SE0820430.RData')

spatial_nature2000_original_SE0820430$Area_laea <- raster::area(spatial_nature2000_original_SE0820430) #Take care of units # In sq meters
spatial_nature2000_original_SE0820430@data

#Project and calculate the area into Mollweide  
#Define the target CRS for the Mollweide equal-area projection
equal_area_crs <- CRS("+proj=moll +datum=WGS84 +units=m +no_defs")

#Reproject the SpatialPolygonsDataFrame
spatial_nature2000_original_SE0820430_moll <- spTransform(spatial_nature2000_original_SE0820430, equal_area_crs)

#Calculate the area
spatial_nature2000_original_SE0820430_moll$Area_moll <- raster::area(spatial_nature2000_original_SE0820430_moll) 
spatial_nature2000_original_SE0820430_moll@data

#Project and calculate the area into Cilyndrical Equal Area  
#Define the target CRS for Cylindrical Equal Area 
cea_proj <- "+proj=cea +lon_0=0 +lat_ts=0 +x_0=0 +y_0=0 +ellps=WGS84 +towgs84=0,0,0,-0,-0,-0,0 +units=m +no_defs" 

#Reproject the SpatialPolygonsDataFrame
spatial_nature2000_original_SE0820430_cea <- spTransform(spatial_nature2000_original_SE0820430, cea_proj)

#Calculate the area
spatial_nature2000_original_SE0820430_cea$Area_cea <- raster::area(spatial_nature2000_original_SE0820430_cea)
spatial_nature2000_original_SE0820430_cea@data

#Plot the three projections  
par(mfrow = c(3,1))
plot(spatial_nature2000_original_SE0820430, main = "SE0820430_in_LAEA", col = "blue", border = NA) 
plot(spatial_nature2000_original_SE0820430_moll, main = "Polygon SE0820430_in_Mollweide", col = "blue", border = NA)
plot(spatial_nature2000_original_SE0820430_cea, main = "Polygon SE0820430_in_CEA", col = "blue", border = NA)

• SE0820430_in_LAEA is the plot of the original SpatialPolygonsDataFrame,
• Polygon SE0820430_in_Mollweide represents the result of the projected SE0820430_in_LAEA in Mollweide equal-area projection and
• Polygon SE0820430_in_CEA represents the result of the projected SE0820430_in_LAEA in Cylindrical equal-area projection

I tried with two different projections and it didn't work. I have changed the projection in ArcGIS 10.4 and it works - the area is just 16 square meters different and the shape is correct.

What I am reporting here is a specific case as an example to simplify the explanation, but the error occurred across hundreds of different protected areas in the Natura 2000 dataset. The original shapefile contains 27,020 features, it can be downloaded from this European Environment Agency link. When I reprojected the original shapefile, I checked it and most of the polygons were correctly reprojected. Later, I compared the areas with all other areas and detected these errors.

In addition, I ran the same code using the st_transform() function from the sf package, but I got the same error. I guess that I am doing something wrong applying the spTransform() function but I can't find what it is.

Here is the structure of the objects to help reproduce it (I can't include all due to space limitations):

str(spatial_nature2000_original_SE0820430)#@ Polygons :List of 7591
# Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
#   ..@ data       :'data.frame': 1 obs. of  7 variables:
#   .. ..$ SITECODE  : chr "SE0820430"
#   .. ..$ SITENAME  : chr "Torne och Kalix älvsystem"
#   .. ..$ RELEASE_DA: chr "2021/10/04"
#   .. ..$ MS        : chr "SE"
#   .. ..$ SITETYPE  : chr "B"
#   .. ..$ INSPIRE_ID: chr "SE.SWEPA.SE0820430"
#   .. ..$ Area_laea : num 1.75e+09
#   ..@ polygons   :List of 1
#   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
#   .. .. .. ..@ Polygons :List of 7591
#   .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
#   .. .. .. .. .. .. ..@ labpt  : num [1:2] 4864853 4959222
#   .. .. .. .. .. .. ..@ area   : num 4.3e+09
#   .. .. .. .. .. .. ..@ hole   : logi FALSE
#   .. .. .. .. .. .. ..@ ringDir: int 1
#   .. .. .. .. .. .. ..@ coords : num [1:1202650, 1:2] 4930104 4930080 4930072 4930064 4930057 ...

str(spatial_nature2000_original_SE0820430_moll)#@ polygons   :List of 1
# Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
#   ..@ data       :'data.frame': 1 obs. of  8 variables:
#   .. ..$ SITECODE  : chr "SE0820430"
#   .. ..$ SITENAME  : chr "Torne och Kalix älvsystem"
#   .. ..$ RELEASE_DA: chr "2021/10/04"
#   .. ..$ MS        : chr "SE"
#   .. ..$ SITETYPE  : chr "B"
#   .. ..$ INSPIRE_ID: chr "SE.SWEPA.SE0820430"
#   .. ..$ Area_laea : num 1.75e+09
#   .. ..$ Area_moll : num 12.9
#   ..@ polygons   :List of 1
#   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
#   .. .. .. ..@ Polygons :List of 1
#   .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
#   .. .. .. .. .. .. ..@ labpt  : num [1:2] 974250 7586174
#   .. .. .. .. .. .. ..@ area   : num 12.9
#   .. .. .. .. .. .. ..@ hole   : logi FALSE
#   .. .. .. .. .. .. ..@ ringDir: int 1
#   .. .. .. .. .. .. ..@ coords : num [1:4, 1:2] 974241 974258 974250 974241 7586175 ...

str(spatial_nature2000_original_SE0820430_cea)#@ polygons   :List of 1
# Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
#   ..@ data       :'data.frame': 1 obs. of  8 variables:
#   .. ..$ SITECODE  : chr "SE0820430"
#   .. ..$ SITENAME  : chr "Torne och Kalix älvsystem"
#   .. ..$ RELEASE_DA: chr "2021/10/04"
#   .. ..$ MS        : chr "SE"
#   .. ..$ SITETYPE  : chr "B"
#   .. ..$ INSPIRE_ID: chr "SE.SWEPA.SE0820430"
#   .. ..$ Area_laea : num 1.75e+09
#   .. ..$ Area_cea  : num 13
#   ..@ polygons   :List of 1
#   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
#   .. .. .. ..@ Polygons :List of 1
#   .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
#   .. .. .. .. .. .. ..@ labpt  : num [1:2] 2000291 5887632
#   .. .. .. .. .. .. ..@ area   : num 13
#   .. .. .. .. .. .. ..@ hole   : logi FALSE
#   .. .. .. .. .. .. ..@ ringDir: int 1
#   .. .. .. .. .. .. ..@ coords : num [1:4, 1:2] 2000274 2000307 2000291 2000274 5887633 ...

Answer 1

I suspect your issue stems from applying raster::area() to a spatial vector object. Also, unless you need raster and/or sp objects, it is a good idea to future-proof your workflow. That's because both the raster and sp packages were deprecated in October 2023, and have been replaced by terra and sf.

This solution will produce your desired outcome, and if you need your data as SpatialPolygonsDataFrame objects, you can convert them using sf::as_Spatial() at the end.

As for the discrepancy between the registered area and the area returned by R functions (for the data in its native CRS), I suggest you contact the authors of the data to find out how their value was derived.

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

# Load data from working directory previously downloaded and unzipped from
# https://www.eea.europa.eu/data-and-maps/data/natura-14/natura-2000-spatial-data
spatial_nature2000_original <- st_read("Natura2000_end2021_rev1_epsg3035.shp")

# Filter by SITECODE and add area column
spatial_nature2000_original_SE0820430 <- spatial_nature2000_original %>%
  filter(SITECODE == "SE0820430") %>%
  mutate(Area_laea = st_area(.))

spatial_nature2000_original_SE0820430[, c("SITECODE", "Area_laea")]
# Simple feature collection with 1 feature and 2 fields
# Geometry type: MULTIPOLYGON
# Dimension:     XY
# Bounding box:  xmin: 4655780 ymin: 4803259 xmax: 4967643 ymax: 5134291
# Projected CRS: ETRS89-extended / LAEA Europe
#    SITECODE        Area_laea                       geometry
# 1 SE0820430 1753394506 [m^2] MULTIPOLYGON (((4930104 487...

# Define Mollweide PROJ4 string
equal_area_crs <- "+proj=moll +datum=WGS84 +units=m +no_defs"
# or search for equivalent EPSG/ESRI code
equal_area_crs <- "ESRI:53009"

# Project and calculate area
spatial_nature2000_original_SE0820430_moll <- spatial_nature2000_original_SE0820430 %>%
  st_transform(equal_area_crs) %>%
  mutate(Area_moll = st_area(.))

spatial_nature2000_original_SE0820430_moll[, c("SITECODE", "Area_moll")]
# Simple feature collection with 1 feature and 2 fields
# Geometry type: MULTIPOLYGON
# Dimension:     XY
# Bounding box:  xmin: 962626.4 ymin: 7414242 xmax: 1378628 ymax: 7693381
# Projected CRS: +proj=moll +datum=WGS84 +units=m +no_defs
#    SITECODE        Area_moll                       geometry
# 1 SE0820430 1745038571 [m^2] MULTIPOLYGON (((1328694 746...

# Define Cylindrical Equal Area PROJ4 string
cea_proj <- "+proj=cea +lon_0=0 +lat_ts=0 +x_0=0 +y_0=0 +ellps=WGS84 +towgs84=0,0,0,-0,-0,-0,0 +units=m +no_defs" 

# Project and calculate area
spatial_nature2000_original_SE0820430_cea <- spatial_nature2000_original_SE0820430 %>%
  st_transform(cea_proj) %>%
  mutate(Area_cea = st_area(.))

spatial_nature2000_original_SE0820430_cea[, c("SITECODE", "Area_cea")]
# Simple feature collection with 1 feature and 2 fields
# Geometry type: MULTIPOLYGON
# Dimension:     XY
# Bounding box:  xmin: 1996714 ymin: 5801257 xmax: 2688687 ymax: 5939193
# Projected CRS: +proj=cea +lon_0=0 +lat_ts=0 +x_0=0 +y_0=0 +ellps=WGS84 +towgs84=0,0,0,-0,-0,-0,0 +units=m +no_defs
#    SITECODE         Area_cea                       geometry
# 1 SE0820430 1753397068 [m^2] MULTIPOLYGON (((2629712 582...

# Plot
par(mfrow=c(3,1))
plot(st_geometry(spatial_nature2000_original_SE0820430),
     main = "SE0820430_in_LAEA", col = "blue", border = NA) 
plot(st_geometry(spatial_nature2000_original_SE0820430_moll),
     main = "Polygon SE0820430_in_Mollweide", col = "blue", border = NA)
plot(st_geometry(spatial_nature2000_original_SE0820430_cea),
     main = "Polygon SE0820430_in_CEA", col = "blue", border = NA)