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

Question

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

我有一个来自 WorldPop 的人口数据栅格文件，以及一组叠加在上面的纬度/经度位置点。我正在尝试确定在这些兴趣点的给定距离内有多少人口（根据 WorldPop 估计），还有哪些人口不在。

我知道使用 raster::extract，我应该能够从（例如）每个点周围 1 公里的缓冲区中获得人口值的总和。（虽然我的点和栅格数据都在 lat/lon 投影中，所以我想我需要首先通过将投影更改为 utm 来纠正这一点 here。）

但是，由于这些点中的某些点相距不到 1 公里，我担心这个总和会重复计算缓冲区重叠的某些单元格的人口。缓冲是否会自动纠正这一点，或者是否有一种有效的方法来确保不是这种情况，并且还可以从缓冲点区域选择的倒数中获取值？

Answer 1

好吧，感谢 suggestion via Twitter and 围绕点创建 SpatialPolygons，我已经找到了这个问题的答案。这可能不是最有效的方法 - 它在大多边形上被证明非常慢 - 但它对我的目的是可行的。

示例代码如下：

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

Answer 2

这是一个最小的独立可复制示例，

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 你需要这样做（待修复）。

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

Using R to extract sums from a raster layer, in areas outside potentially overlapping buffers

gis

r

raster