多边形的表面积加权空间连接

Question

问题陈述

我有两组多边形，我想将其中一组多边形的数量特征合并到另一组多边形中。

例如，考虑Yolo县的多边形，yolo。我想将所有符合城市多边形davis内字段estimate的收缩级别数据聚合起来。

结果应该是一个新的多边形davis，具有一个新的字段estimate，它是yolo中所有落在davis范围内的要素的表面积加权estimate。 我如何在sp或sf中实现这个功能？

---

可重复的示例

从此网站下载的城市多边形davis，文件：CityLimits.zip。

# packages
library(tidycensus)
library(tidyverse)
library(raster)

# get tract level data for yolo county
yolo  <- get_acs(state = "CA", county = "Yolo", geography = "tract", 
                 variables = "B19013_001", geometry = TRUE)


# city of davis shapefile
davis <- raster::shapefile("Davis_DBO_CityLimits.shp")
davis <- davis %>% spTransform(., st_crs(yolo)$`proj4string` %>% crs())
davis <- st_as_sf(davis)
yolo <- yolo %>% st_transform(st_crs(davis)$`proj4string`) 

# plot
ggplot() + 
  geom_sf(data = yolo, aes(fill = estimate)) +
  geom_sf(data = davis, alpha = 0.3, color = "red") +
  coord_sf(xlim=c(-121.6, -121.9), ylim = c(38.5, 38.6))

---

注意: 我看到了这篇SO帖子，但是由于链接失效所以无法重现。

- Rich Pauloo

---

没有尝试过，sf::st_interpolate_aw(yolo, davies) 应该可以实现你想要的功能。 - TimSalabim

那也是我的第一想法。 - Edzer Pebesma

嗨，Tim和Edzer，这正是我正在寻找的函数。谢谢！如果你想将其写成答案，我会接受它，以便未来遇到此问题的用户可以使用。我还发现?st_interpolate_aw中的示例非常有用，比我的更具可重复性。 - Rich Pauloo

Answer 1

我来举个可重复的例子。

数据示例：

library(raster)
p <- shapefile(system.file("external/lux.shp", package="raster"))
p$value <- 1:length(p)
b <- as(extent(6, 6.4, 49.76, 50), 'SpatialPolygons')
b <- SpatialPolygonsDataFrame(b, data.frame(bid = 1))
crs(b) <- crs(p)
plot(p)
plot(b, add=T, border='red', lwd=2)

首先是“手动方式”

i <- intersect(b, p)    
i$AREA <- area(i) / 1000000
aw <- sum(i$AREA * i$value) / sum(i$AREA)
aw
# 5.086891

sp 方法：

a <- aggregate(p['value'], b, FUN=sum, areaWeighted=TRUE)
a$value
# 5.085438

现在有了sf。

library(sf)
pf <- as(p, 'sf')
bf <- as(b, 'sf')

x <- sf::st_interpolate_aw(pf['value'], bf, extensive=F)
x$value
# 5.086891

- Robert Hijmans