使用 spTransform() 重新投影 SpatialPolygonsDataFrame 时结果不正确

Question

R版本4.4.0（2024-04-24 ucrt）——“小狗杯” 平台：x86_64-w64-mingw32/x64

我有一个非常复杂的 SpatialPolygonsDataFrame，代表来自瑞典的 Natura 2000 站点，SE0820430，以黑色显示：

首先，我检查了原始投影（LAEA）的面积，为 1753394516 平方米。根据此 BISE 页面，注册面积为 1760923000 平方米，比较相似（注册面积的 99.57%）。

问题是，当我运行下面的代码来转换投影时，我看到创建了一个新的多边形/形状，它与原始的 SpatialPolygonsDataFrame 没有任何共同点（参见最后的图）。毫不奇怪，面积相差很大（约 13 平方米）。我没有收到警告或错误。

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 是原始 SpatialPolygonsDataFrame 的绘图，
• Polygon SE0820430_in_Mollweide 表示在 Mollweide 等面积投影中投影 SE0820430_in_LAEA 的结果，并且
• Polygon SE0820430_in_CEA 表示 SE0820430_in_LAEA 在圆柱等积投影中的投影结果

我尝试了两种不同的投影，但没有成功。我已经在 ArcGIS 10.4 中更改了投影并且它有效 - 面积仅相差 16 平方米，并且形状是正确的。

我在这里报告的是一个具体案例作为示例以简化说明，但错误发生在 Natura 2000 数据集中的数百个不同的保护区中。原始 shapefile 包含 27,020 个特征，可以从此欧洲环境署链接下载。当我重新投影原始形状文件时，我检查了它，大多数多边形都被正确地重新投影。后来，我将这些区域与所有其他区域进行了比较，并发现了这些错误。

此外，我使用

st_transform()

包中的

sf

函数运行相同的代码，但出现了相同的错误。我想我在应用

spTransform()

函数时做错了什么，但我找不到它是什么。

这是帮助重现它的对象结构（由于空间限制，我无法包含所有内容）：

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

我怀疑您的问题源于将

raster::area()

应用于空间矢量对象。另外，除非您需要

raster

和/或

sp

对象，否则最好让您的工作流程面向未来。这是因为

raster

和

sp

软件包均已于 2023 年 10 月弃用，并已被

terra

和

sf

取代。

此解决方案将产生您想要的结果，如果您需要将数据作为 SpatialPolygonsDataFrame 对象，您可以在最后使用

sf::as_Spatial()

进行转换。

至于注册的区域和R函数返回的区域（对于其原生CRS中的数据）之间的差异，我建议您联系数据的作者，了解其值是如何得出的。

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)