使用 R 从栅格图层中提取可能重叠缓冲区外部区域中的总和

Question

我对栅格数据和使用 R 进行空间数据分析非常陌生，所以如果我错过了一个明显的解决方案或过程，我深表歉意。

我有一个来自 WorldPop 的人口数据的光栅文件，以及一组覆盖在其上的纬度/经度位置点。我试图确定(根据WorldPop估计)在这些兴趣点的给定距离内的人口比例，以及哪一部分不是。

我知道使用 raster：：extract，我应该能够从(例如)每个点周围的 1 公里缓冲区获取人口值的总和。(虽然我的点和栅格数据都在纬度/纬度投影中，所以我收集我需要首先通过将投影更改为utm来更正这一点，如此处所示。

但是，由于这些点中的某些点之间的距离不到 1 公里，我担心此总和会重复计算缓冲区重叠的某些像元的总体。缓冲是否会自动纠正此问题，或者是否有一种有效的方法来确保情况并非如此，并且还可以从缓冲点区域选择的反函数中获取值？

Answer 1

这是一个最小的自包含可重现示例，

library(raster)
r <- raster(system.file("external/rlogo.grd", package="raster"))
d <- matrix(c(48, 48, 48, 53, 50, 46, 54, 70, 84, 85, 74, 84, 95, 85, 
66, 42, 26, 4, 19, 17, 7, 14, 26, 29, 39, 45, 51, 56, 46, 38, 31, 
22, 34, 60, 70, 73, 63, 46, 43, 28), ncol=2)
p <- SpatialPoints(d, proj4string=crs(r))     

具有点p和栅格r的简单工作流将是

b <- buffer(p, 10)
m <- mask(r, b)
ms <- cellStats(m, "sum")
rs <- cellStats(r, "sum")
ms/rs
#[1] 0.4965083

或者你可以使用terra，使它走得更快，像这样

library(terra)
r <- rast(system.file("ex/logo.tif", package="terra")) [[1]]
p <- vect(d, crs=crs(r))
b <- buffer(p, 10)
m <- mask(r, b)
ms <- global(m, "sum", na.rm=TRUE)
rs <- global(r, "sum")
ms/rs

顺便说一下，对于raster包，您关于需要转换 lon/lat 数据的断言对于extract或buffer是不正确的。相比之下，对于terra，您需要这样做(要修复)。

Answer 2

好吧，多亏了通过Twitter提出的建议和这份围绕点创建空间多边形的指南，我已经能够找到答案。这可能不是最有效的方法——事实证明，它在大型多边形上非常慢——但它对我的目的来说是可行的。

下面是示例代码：

library(tabularaster)
library(raster)
library(tidyverse)
library(geos)
# -----------------------
# load point data ---
p <- read_csv("points_of_interest.csv")
p_df <- p %>% rename(x = lat, y = lon)
p_coords <- p_df[, c("y","x")]
p_spdf <- SpatialPointsDataFrame(
coords = pc_coords,
data = p_df,
proj4string = CRS("+init=epsg:4326"))
# convert projection to metric units
p_mrc <- spTransform(
p_spdf,
CRS("+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 
+x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs")
)
# buffer to 1000 meters
p_mrc_1k_mrc <- gBuffer(
p_mrc, byid = TRUE, width = 1000)
# switch back to lat/lon
p_mrc_1k <- spTransform(p_mrc_1k_mrc, CRS("+init=epsg:4326"))
# load raster data -------
r <- raster("pop.tif")
r_tib <- tabularaster::as_tibble(r)
# get intersection of cells and polygons
cell_df_1k <- cellnumbers(r, p_mrc_1k)
# get list of cells where there is intersection
target_cell_1k <- cell_df_1k$cell_
# add cell values to df listing all cells covered by polys
target_cells_extract_1k <- cell_df_1k %>%
rename(cellindex = cell_) %>%
left_join(r_tib)
# calculate the sum of population within 1k radius for each object 
# (this includes overlapping population cells shared between polys)
cell_sum_1k <- target_cells_extract_1k %>%
group_by(object_) %>%
summarize(pop_1k = sum(cellvalue, na.rm = T))
# get only unique cell values for total overlapping coverage of all polys
target_cells_unique_1k <- r_tib %>% filter(cellindex %in% target_cell_1k)
total_coverage_pop <- sum(target_cells_unique_1k$cellvalue, na.rm = T)
outside_coverage_pop <- sum(r_tib$cellvalue) - total_coverage_pop