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)
我尝试了两种不同的投影,但没有成功。我已经在 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 ...
我怀疑您的问题源于将
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)