Неверный результат при перепроецировании SpatialPolygonsDataFrame с помощью spTransform()

Версия R 4.4.0 (24 апреля 2024 г., ucrt) – "Кубок щенка" Платформа: x86_64-w64-mingw32/x64

У меня есть очень сложный SpatialPolygonsDataFrame, представляющий сайт Natura 2000 из Швеции, SE0820430, показанный черным цветом:

Неверный результат при перепроецировании SpatialPolygonsDataFrame с помощью spTransform()

Сначала я проверил площадь исходной проекции (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,
  • Полигон SE0820430_in_Mollweide представляет собой результат спроецированного SE0820430_in_LAEA в равновеликой проекции Моллвейде и
  • Полигон SE0820430_in_CEA представляет собой результат проецирования SE0820430_in_LAEA в цилиндрической равновеликой проекции.

Я пробовал использовать две разные проекции, но это не сработало. Я изменил проекцию в ArcGIS 10.4, и она работает — площадь отличается всего на 16 квадратных метров и форма правильная.

Здесь я рассказываю о конкретном случае в качестве примера для упрощения объяснения, но ошибка произошла в сотнях различных охраняемых территорий в наборе данных Natura 2000. Исходный шейп-файл содержит 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 ...

1
68
1

Ответ:

Решено

Я подозреваю, что ваша проблема связана с применением raster::area() к объекту пространственного вектора. Кроме того, если вам не нужны объекты raster и/или sp, рекомендуется подготовить свой рабочий процесс к будущему. Это связано с тем, что оба пакета raster и sp устарели в октябре 2023 года и были заменены 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)