Неверный результат при перепроецировании SpatialPolygonsDataFrame с помощью spTransform()
Версия R 4.4.0 (24 апреля 2024 г., 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,
- Полигон 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) Где/как я получил исходный шейп-файл "Natura2000_end2021_rev1_epsg3035.shp" 2) Код для получения объекта "spatial_nature2000_SE0820430" (теперь называется "spatial_nature2000_original_SE0820430") 3) Структура объектов 4) Комментарий, указывающий, что это не единичный случай, а я использую его в качестве примера, чтобы лучше объяснить и обнаружить проблему.
Извините, вы правы. Я отредактировал и добавил строку cea_proj. Спасибо всем вам за вашу помощь.
Ответы 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)
Привет, Пау Люк! Вы можете помочь воспроизвести проблему, добавив вывод
str(spatial_nature2000_SE0820430)и код, который вы использовали для получения из ресурса EUNIS (или другого) в свой пространственный объектspatial_nature2000_SE0820430.