Использование 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, чтобы сделать это go быстрее, вот так

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 ваше утверждение о необходимости преобразования данных долгота / долгота неверно для extract или buffer. Напротив, с terra вам нужно это сделать (чтобы исправить).

Answer 2

Что ж, благодаря предложению через Twitter и этого руководства по созданию 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