两个形状文件之间的空间连接导致不匹配

Question

我有两个形状文件，对于世界不同地区具有不同的名称。我需要

countries

中多边形的名称，它们与

countries_GIFT

中多边形的名称相对应。我已经用其他形状文件完成了此操作，并且一直运行没有问题，但由于某种原因，我在两个数据帧之间出现不匹配。

例如，第一个形状文件中的俄罗斯条目在另一个形状文件中被识别为爱沙尼亚。当检查 x y 坐标时，我可以看到该点确实在俄罗斯，但由于某种原因它被提取到其他地方。

我已经检查了两个shapefile的投影以及区域是否相互重叠，所以我不确定问题是什么。

这是我的代码

library(sf)
library(tidyverse)
library(GIFT)

setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

sf_use_s2(FALSE)

countries <- st_read("../../../../locations/Shapefiles/Locations.shp")
countries <- st_transform(countries, crs = 4326)

GIFT_shapes <- GIFT_shapes() # retrieves all shapefiles by default

countries_GIFT <- GIFT_shapes() 
countries_GIFT_crs <- st_transform(countries_GIFT, crs = 4326)

gift_stand_countries <- st_join(countries_GIFT_crs, countries, join = st_intersects)

# Check if regions overlap
ggplot() +
  geom_sf(data = countries_GIFT_crs, fill = "lightblue", color = "black") +
  geom_sf(data = countries, fill = NA, color = "red")

这是两个形状文件的绘图：

以及两个 df 的几何形状

head(countries$geometry)
Geometry set for 6 features 
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -178.8733 ymin: -14.54889 xmax: 74.89291 ymax: 71.33334
Geodetic CRS:  WGS 84
First 5 geometries:
MULTIPOLYGON (((71.00423 38.47541, 71.02091 38....
MULTIPOLYGON (((20.97472 59.96736, 20.98694 59....
MULTIPOLYGON (((-178.8139 51.8375, -178.7424 51...
MULTIPOLYGON (((19.59271 41.63155, 19.59014 41....
MULTIPOLYGON (((-1.630694 35.21958, -1.630972 3...

head(countries_GIFT_crs$geometry)
Geometry set for 6 features 
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 83.6
Geodetic CRS:  WGS 84
First 5 geometries:
MULTIPOLYGON (((-81.81 24.54, -81.77 24.54, -81...
MULTIPOLYGON (((145.99 43.37, 146.01 43.36, 146...
MULTIPOLYGON (((140.87 51.42, 140.88 51.43, 140...
MULTIPOLYGON (((166.4939 -78.21118, 166.6582 -7...
MULTIPOLYGON (((145.99 43.37, 146.01 43.36, 146...

提前非常感谢

Answer 1

如果无法访问这两个来源，就很难指出具体的东西，但您很可能正在处理具有不同分辨率和精度的数据集，后者也可以从对象的几何形状中看出，

71.00423 38.47541

与

-81.81 24.54

。即使具有相同/相似的分辨率和精度，也不太可能获得来自不同数据源的形状的完美匹配，因此预计与邻居相交（+1 @VinceGreg 的评论）。

为了说明，我们可以以爱沙尼亚为

GIFT_shapes()

，检查 GISCO 中的哪些中分辨率形状相交，以及更高分辨率的国家轮廓以供参考：

library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
library(GIFT)
library(ggplot2)

# 2016; resolution: 1:20million
countries <- giscoR::gisco_countries[, "CNTR_ID"]

# Estonia
GIFT_est <- GIFT_shapes(entity_ID = 11456)
#> You are asking for the latest stable version of GIFT which is 3.2.
#> ================================================================================

intersecting_countries <- 
  st_join(
    countries,
    GIFT_est,
    left = FALSE
  ) |> 
  # crop with extended bbox for plot
  st_crop(st_bbox(GIFT_est) + c(-.2, -.2, .2, .2))
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries
intersecting_countries[, 1:2]
#> Simple feature collection with 3 features and 2 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 21.56417 ymin: 57.31555 xmax: 28.40378 ymax: 60.05231
#> Geodetic CRS:  WGS 84
#>     CNTR_ID entity_ID                       geometry
#> 82       EE     11456 MULTIPOLYGON (((27.83301 59...
#> 138      LV     11456 MULTIPOLYGON (((27.75134 57...
#> 185      RU     11456 MULTIPOLYGON (((28.40378 59...

ggplot() +
  geom_sf(data = intersecting_countries, aes(fill = CNTR_ID)) +
  geom_sf(data = GIFT_est, fill = "white", alpha = .6 ) +
  # 1:1million outline for reference
  geom_sf(data = giscoR::gisco_get_countries(resolution = "01", country = "Estonia"), color = "navy", fill = NA) +
  scale_fill_viridis_d() +
  coord_sf(expand = FALSE)

或者您可以查看 https://gift.uni-goettingen.de/map .

---

您可以通过首先从一个数据集的多边形生成st_point_on_surface()点（不保证质心最终位于多边形中）并将这些点连接到另一个数据集的多边形来获得所需的结果：

GIFT_est |> 
  st_point_on_surface() |> 
  st_join(giscoR::gisco_countries)
#> Warning: st_point_on_surface assumes attributes are constant over geometries
#> Warning in st_point_on_surface.sfc(st_geometry(x)): st_point_on_surface may not
#> give correct results for longitude/latitude data
#> Simple feature collection with 1 feature and 17 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 25.50592 ymin: 58.62939 xmax: 25.50592 ymax: 58.62939
#> Geodetic CRS:  WGS 84
#>   entity_ID geo_entity  point_x  point_y     area    x_min    x_max    y_min
#> 1     11456    Estonia 25.54518 58.67211 45483.44 21.76361 28.21478 57.51511
#>      y_max entity_class    entity_type polygon_source CNTR_ID CNTR_NAME
#> 1 59.82264     Mainland Political Unit           GADM      EE     Eesti
#>   ISO3_CODE NAME_ENGL FID                  geometry
#> 1       EST   Estonia  EE POINT (25.50592 58.62939)

对于预计的 CRS，我希望它能够很好地工作，并且具有地理坐标，上面的警告可能不应该太轻视。