st_transform 似乎改变了 shapefile 的几何区域

Posted 2023-03-17

【中文标题】st_transform 似乎改变了 shapefile 的几何区域

【英文标题】：st_transform appears to change the area of the geometry of a shapefile

【发布时间】：2022-01-19 17:56:45

【问题描述】：

转换 shapefile 投影时，多边形的面积不应该保持不变吗？我目前正在为国家 GIN 使用 shapefile (https://data.humdata.org/dataset/f814a950-4d4e-4f46-a880-4da5522f14c4/resource/731e11cb-be02-46cf-8347-db0a07abff4e/download/gin_admbnda_adm4_2021_ocha.zip)，在我使用 R 中的 st_read 函数将其读取后，我将其称为 gin_shp，如下所示：

gin_shp <- st_read(dsn = "INPUT FOLDER", layer = "sous_prefectures") 
gin_shp ##this shows the crs is WGS84 obviously in degrees

##compute the area of the multipolygon geometries
gin_shp$area <- sf::st_area(gin_shp) ##computes areas in m^2
gin_shp$area <- units::set_units(gin_shp$area, "km^2") ##converting the area to square kms

sum(gin_shp$area) ##compute total area of the entire country (google seems to agree!)
245084.2 [km^2]

我需要创建 1sqkm 网格（使用 sf::st_make_grid），所以我尝试通过调用 st_transform 函数将 crs 投影从度数转换为 UTM：

crs_dt <- rgdal::make_EPSG() ##first load the dataframe of CRS projections available in R

gin_shp <- sf::st_transform(gin_shp, crs = crs_dt$prj4[crs_dt$code == 4328]) ##select one and assign it to gin_shp

#now I try to compute total area under new projection regime by running the same code as before
gin_shp$area <- sf::st_area(gin_shp) ##computes areas in m^2
gin_shp$area <- units::set_units(gin_shp$area, "km^2") #converting the area to square kms

#compute total area of the entire country 
sum(gin_shp$area) 
44363.83 [km^2] (way way off)

您知道为什么会发生这种情况吗？任何想法如何解决它？

【问题讨论】：

你在st_transform时使用的CRS是什么？

我试了一个号码。我尝试了 crs = 4328，即“+proj=geocent +datum=WGS84 +units=m +no_defs +type=crs”。我也试过 9191 是 "+proj=aea +lat_0=-40 +lon_0=175 +lat_1=-30 +lat_2=-50 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs + type=crs" 并且都没有接近保留该区域

【参考方案1】：

不，多边形的面积不能保证在不同 CRS 之间转换时保持不变。地球是圆的，地图和电脑屏幕是平的——必须付出一些；无论是面积还是形状都必须有所扭曲。

有一些区域保留投影 - 例如Mollweide - 但这些更多的是例外而不是规则。

举个夸张的例子，考虑世界数据集，取自giscoR 包。格陵兰岛（靠近极地）和刚果（位于赤道）在球面上的面积大致相同（在 WGS84 上使用球面几何工具计算得出），但在投影时却大不相同。尤其是投影到墨卡托及其衍生产品（例如谷歌地图使用的 Web 墨卡托）时。

library(sf)
library(dplyr)
library(giscoR)
library(ggplot2)

# the world in, 1: 20M
svet <- gisco_get_countries(resolution = "20")

# Greenland & Congo - to be drawn in red
glmd <- svet %>% 
   filter(CNTR_ID %in% c('GL', 'CD'))

# web mercator = default for google maps and other web based tools
ggplot() +
   geom_sf(data = glmd, fill = "red", color = NA) +
   geom_sf(data = svet, fill = NA, color = "gray45") +
   coord_sf(crs = st_crs("EPSG:3857"),
            ylim = c(-20e6, 20e6))

# Mollweide - equal area projection
ggplot() +
   geom_sf(data = glmd, fill = "red", color = NA) +
   geom_sf(data = svet, fill = NA, color = "gray45") +
   coord_sf(crs = st_crs("ESRI:53009"))

【讨论】：

谢谢。但是，如果 shapefile 以度为单位，那么我如何将 shapefile 细分为 1sqkm 网格。我似乎无法从 R 中应该保留区域的 crs 投影集中找到保留 gin_shp 总面积的投影？

您感兴趣的领域是什么？公里中的网格可能是全球数据集的一个问题（但你不太可能去 1 平方公里）。对于较小的（国家规模及以下）地区，通常有一个官方（或半官方）预测，作为公认的最佳实践。