当我使用 sp 包中的 spTransform 重新投影 SpatialPolygonsDataFrame 时，获得了错误的结果

问题描述 投票：0回答：1

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

我有一个SpatialPolygonsDataFrame“spatial_nature2000_SE0820430”，它非常复杂，它代表来自瑞典的Nature2000站点，SE0820430从波罗的海扩展到与挪威的边境，它以黑色表示：

SE0820430

首先，我检查了原始投影（LAEA）上的面积，得到了一个完全不正确的值 1753394516 平方米（注册面积，基于 https://biodiversity.europa.eu/sites/natura2000 /SE0820430 为 1760923000 平方米），但相对相似（注册面积的 99.57%）。问题是，当我运行下面的代码来转换投影时，我没有收到任何警告，也没有任何错误，当我检查新投影时，我看到已经创建了一个新的多边形/形状，它与原始的没有任何共同点SpatialPolygonsDataFrame，很简单，就像一个三角形（见最后图），不出意外面积差别很大（大约13平方米）。

  library(raster)
  library(sf)
  library(sp)
  
  #We read the spatial information of the Nature2000, downloaded from https://www.eea.europa.eu/data-and-maps/data/natura-14/natura-2000-spatial-data :
  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  
    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 在圆柱等积投影中的投影结果

three projections

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

我在这里报告的是一个具体案例（一个要素代表一个保护区），我用它作为例子来简化说明，但是当数百个要素（数百个保护区）出现错误投影的错误时我重新投影了 Nature2000（欧洲保护区网）的原始形状文件。原始形状文件包含 27,020 个特征（是多边形），可以在 https://www.eea.europa.eu/data-and-maps/data/natura-14/natura-2000-spatial-data 下载。当我重新投影原始形状文件时，我检查了它，大多数多边形都被正确地重新投影，但后来我将这些区域与所有这些区域进行了比较并检测到了这些错误。

此外，我使用“sf”包中的“st_transform”函数运行相同的代码，投影也错误，具有相同的错误。我想我在应用 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 ...
r maps gis map-projections
1个回答
0
投票

raster
和 
sp
软件包均已于 2023 年 10 月弃用，并已被 
terra
和 
sf
取代。该解决方案将产生您想要的结果，如果您需要将数据作为 SpatialPolygonsDataFrame 对象，您可以使用 
sf::as_Spatial()
进行转换。

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("data/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)

1

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.