使用spTransform()重投影SpatialPolygonsDataFrame结果异常

环境信息 #

• R版本：4.4.0 (2024-04-24 ucrt) -- "Puppy Cup"
• 平台：x86_64-w64-mingw32/x64

问题背景 #

我有一个表示瑞典Natura2000站点SE0820430的复杂SpatialPolygonsDataFrame，原图形如下：

原投影为LAEA，计算面积为1753394516平方米，与登记面积1760923000平方米接近（为登记面积的99.57%）。

运行以下代码进行投影转换后，生成的新多边形形状与原数据完全不符，面积仅约13平方米，且无任何警告或错误。

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

#读取Natura2000空间信息：
spatial_nature2000_original <- shapefile("Natura2000_end2021_rev1_epsg3035.shp")

#提取示例保护区站点"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) #注意单位：平方米
spatial_nature2000_original_SE0820430@data

#投影至Mollweide并计算面积
#定义Mollweide等面积投影的目标CRS
equal_area_crs <- CRS("+proj=moll +datum=WGS84 +units=m +no_defs")

#重投影SpatialPolygonsDataFrame
spatial_nature2000_original_SE0820430_moll <- spTransform(spatial_nature2000_original_SE0820430, equal_area_crs)

#计算面积
spatial_nature2000_original_SE0820430_moll$Area_moll <- raster::area(spatial_nature2000_original_SE0820430_moll) 
spatial_nature2000_original_SE0820430_moll@data

#投影至圆柱等面积投影并计算面积
#定义圆柱等面积投影的目标CRS
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" 

#重投影SpatialPolygonsDataFrame
spatial_nature2000_original_SE0820430_cea <- spTransform(spatial_nature2000_original_SE0820430, cea_proj)

#计算面积
spatial_nature2000_original_SE0820430_cea$Area_cea <- raster::area(spatial_nature2000_original_SE0820430_cea)
spatial_nature2000_original_SE0820430_cea@data

#绘制三种投影结果
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等面积投影的结果
• Polygon SE0820430_in_CEA：投影至圆柱等面积投影的结果

两种投影转换均失败，但在ArcGIS 10.4中转换正常，面积仅相差16平方米且形状正确。该异常并非个例，实际出现在Natura2000数据集的数百个保护区中。原始shapefile包含27020个要素，重投影后多数多边形正常，对比面积时发现异常。

此外，使用sf包的st_transform()函数运行相同操作，也出现相同错误。

对象结构（部分展示） #

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

内容的提问来源于stack exchange，提问作者Pau Luc