在潜在重叠的缓冲区之外的区域中，使用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"))
p <- SpatialPoints(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))
crs(p) <- crs(r)

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

library(raster)
b <- buffer(p, 1000)
m <- mask(r, b)
ms <- cellStats(m, "sum")
rs <- cellStats(r, "sum")
ms/rs