多边形的表面积加权空间连接
Surface Area Weighted Spatial Join of Polygons
问题陈述
我有两组多边形,我想将一组多边形的定量特征连接到另一组多边形中。
例如,考虑 Yolo 县的多面体 yolo
。我想从适合 davis
.
戴维斯市的多边形内的字段 estimate
中的所有要素汇总道级数据
结果应该是一个 davis
的多边形,带有一个新字段 estimate
,即 中所有要素的 estimate
表面积加权 estimate
11=] 属于 davis
。我如何在 sp
或 sf
中执行此操作?
可重现的例子
从 this website, file: CityLimits.zip
下载的戴维斯市多边形 (davis
)。
# 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))
注:看过this SO post。死链接使其不可复制。
这是一个可重现的例子。
示例数据:
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)
第一个'by hand'
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
问题陈述
我有两组多边形,我想将一组多边形的定量特征连接到另一组多边形中。
例如,考虑 Yolo 县的多面体 yolo
。我想从适合 davis
.
estimate
中的所有要素汇总道级数据
结果应该是一个 davis
的多边形,带有一个新字段 estimate
,即 中所有要素的 estimate
表面积加权 estimate
11=] 属于 davis
。我如何在 sp
或 sf
中执行此操作?
可重现的例子
从 this website, file: CityLimits.zip
下载的戴维斯市多边形 (davis
)。
# 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))
注:看过this SO post。死链接使其不可复制。
这是一个可重现的例子。
示例数据:
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)
第一个'by hand'
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