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

Question

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

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

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

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

Answer 1

好吧，感谢 suggestion via Twitter 和 this guide 在点周围创建 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

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