计算形状文件的交集

Question

我合并了两个形状文件：一个邮政编码 SF 和一个德国投票区的 SF。许多邮政编码位于一个投票区，有些则位于两个（或更多）投票区的边界。我需要找出投票区（占总邮政编码区域）的比例。

我用 sf_intersect 尝试过：

overlap_areas <- st_intersection( shape_wahl,shape_PLZ) %>%
  mutate(overlap_area2 = st_area(.))

但问题是，这段代码计算的是整个交点，而不仅仅是那些非 100% 位于投票区的交点（灰色和黑色轮廓区域相交的区域）

我绘制了“交叉点”，如图所示，整个地图都是红色的

pi <- st_intersection(shape_wahl,shape_PLZ)

ggplot() +
  geom_sf(data = shape_PLZ, aes(fill = category), color = "lightgrey") +
  geom_sf(data = shape_wahl, aes(), color = "black", fill = NA)+
  geom_sf(data = pi, aes(), fill = "red")

最后我想要这样的东西：

邮编	投票区	百分比
12345	1	0.34
12345	2	0.66

我将非常感谢一些提示:)

Answer 1

您可以向您的ZIP对象添加区域，经过

st_intersecation()

之后，您可以执行相同的操作来获取新区域，从那里您可以获得比率。

通用示例，在本例中基于

sf

包中的示例数据集，可能如下所示：

library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
library(dplyr, warn.conflicts = FALSE)
library(ggplot2)

sf_1 <- 
  read_sf(system.file("shape/nc.shp", package="sf"))[, "NAME"] |>
  mutate(area_orig = st_area(geometry))
print(sf_1, n = 3)
#> Simple feature collection with 100 features and 2 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> Geodetic CRS:  NAD27
#> # A tibble: 100 × 3
#>   NAME                                                        geometry area_orig
#> * <chr>                                             <MULTIPOLYGON [°]>     [m^2]
#> 1 Ashe      (((-81.47276 36.23436, -81.54084 36.27251, -81.56198 36.2…    1.14e9
#> 2 Alleghany (((-81.23989 36.36536, -81.24069 36.37942, -81.26284 36.4…    6.11e8
#> 3 Surry     (((-80.45634 36.24256, -80.47639 36.25473, -80.53688 36.2…    1.42e9
#> # ℹ 97 more rows

sf_2 <- 
  st_sf(geometry = st_make_grid(sf_1, n = c(9,3))) |> 
  mutate(id = row_number())
print(sf_2, n = 3)
#> Simple feature collection with 27 features and 1 field
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> Geodetic CRS:  NAD27
#> First 3 features:
#>                         geometry id
#> 1 POLYGON ((-84.32385 33.8819...  1
#> 2 POLYGON ((-83.33864 33.8819...  2
#> 3 POLYGON ((-82.35344 33.8819...  3

sf_3 <- 
  st_intersection(sf_2, sf_1) |>
  mutate(percentage = as.numeric(st_area(geometry) / area_orig)) 
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries
print(sf_3, n = 3)
#> Simple feature collection with 194 features and 4 fields
#> Geometry type: GEOMETRY
#> Dimension:     XY
#> Bounding box:  xmin: -84.32385 ymin: 33.88252 xmax: -75.45698 ymax: 36.58965
#> Geodetic CRS:  NAD27
#> First 3 features:
#>      id      NAME        area_orig                       geometry percentage
#> 21   21      Ashe 1137107793 [m^2] POLYGON ((-81.45289 36.2395...  0.8582486
#> 22   22      Ashe 1137107793 [m^2] POLYGON ((-81.36823 36.5740...  0.1417514
#> 22.1 22 Alleghany  610916077 [m^2] POLYGON ((-81.17667 36.4154...  1.0000000

ggplot() +
  geom_sf(data = sf_3, aes(fill = percentage)) +
  geom_sf(data = sf_1, color = "red", fill = NA) +
  geom_sf(data = sf_2, color = "gray80", linewidth = 1, fill = NA) +
  scale_fill_viridis_b() +
  theme_minimal()

^{创建于 2024-06-04，使用 reprex v2.1.0}