Как выполнить расчет растра (например, аспект) на подмножестве растра на основе пересечения точек в R

Question

Я работаю с некоторыми растровыми данными в R, используя пакет raster.Я хочу рассчитать и извлечь некоторую географическую информацию (например, уклон, аспект) из растра, но только в определенных точках (у меня также есть некоторые данные в виде SpatialPointsDataFrame, по которым я хочу вычислить уклон / аспект / и т. Д.).Я делаю это для нескольких растров с высоким разрешением, и кажется неэффективным использование ресурсов для вычисления этого для каждой растровой ячейки, когда мне нужно всего лишь 5-10% из них.

Я подумал, что, возможно, функция raster::stackApply могла бы работать, но, похоже, она выполняет вычисления на подмножествах растрового кирпича, а не вычисления на подмножествах одного растра на основе точечных местоположений (пожалуйста, исправьте меня, если я ошибаюсь).Я также подумал, что мог бы сделать цикл for, где я извлекаю окружающие ячейки, ближайшие к каждой точке интереса, и таким образом итеративно вычисляю наклон / аспект.Это кажется неуклюжим, и я надеялся на более элегантное или встроенное решение, но оно должно работать.

Это мои мысли о цикле for, но я не уверен, как лучше это сделать.

# Attach packages
library(rgdal)
library(raster)

# Generate example raster data
r = raster()
set.seed(0)
values(r) = runif(ncell(r), min = 0, max = 1000)

# Generate example point data
df.sp = SpatialPoints(
  coords = cbind(runif(25, min = -100, max = 100),  
                 runif(25, min = -50, max = 50)), 
  proj4string = crs(r))

# Iterate on each row of SpatialPoints
for (i in 1:nrow(df.sp)) {

  # Find cell index of current SpatialPoint
  cell.idx = raster::extract(r, df.sp[i,], cellnumbers = TRUE)[1]

  # Find indices of cells surrounding point of interest
  neighbors.idx = raster::adjacent(r, cell.idx, directions = 16)

  # Get DEM values for cell and surrounding cells
  vals.local = r[c(cell.idx, neighbors.idx[,2])]

  # Somehow convert this back to an appropriate georeferenced matrix
  #r.local = ...

  # Perform geometric calculations on local raster
  #r.stack = terrain(r.local, opt = c('slope', 'aspect'))

  # Remaining data extraction, etc. (I can take it from here...)
}

Итак, мне нужен метод для вычисления наклона и аспекта из растра DEM только в определенных точках, заданных объектом SpatialPoints.Если вы знаете о готовом или более элегантном решении, то отлично!Если нет, то некоторая помощь в завершении цикла for (как наилучшим образом извлечь окрестность окружающих ячеек и выполнить вычисления на этом уровне) также была бы наиболее ценной.

Answer 1

Интересный вопрос. Вот возможный подход.

library(raster)

r <- raster()
set.seed(0)
values(r) <- runif(ncell(r), min = 0, max = 1000)
coords <- cbind(runif(25, min = -100, max = 100),  
               runif(25, min = -50, max = 50))

x <- rasterize(coords, r)
f <- focal(x, w=matrix(1, nc=3, nr=3), na.rm=TRUE)             
rr <- mask(r, f)
slope <- terrain(rr, "slope")

extract(slope, coords)
# [1] 0.0019366236 0.0020670699 0.0006305257 0.0025334280 0.0023480935 0.0007527267 0.0002699272 0.0004699626
# [9] 0.0004869054 0.0025651333 0.0010415805 0.0008574920 0.0010664869 0.0017700297 0.0001666226 0.0008405391
#[17] 0.0017682167 0.0009854172 0.0015350466 0.0017714466 0.0012994945 0.0016563132 0.0003276584 0.0020499529
#[25] 0.0006582073

Возможно, не так много прироста эффективности, поскольку он все еще обрабатывает все NA значения

Так что, может быть, так, больше по вашей линии мышления:

cells <- cellFromXY(r, coords)
ngbs <- raster::adjacent(r, cells, pairs=TRUE)

slope <- rep(NA, length(cells))
for (i in 1:length(cells)) {
   ci <- ngbs[ngbs[,1] == cells[i], 2]
   e <- extentFromCells(r, ci)
   x <- crop(r, e)
   slope[i] <- terrain(x, "slope")[5]
}
slope
#[1] 0.0019366236 0.0020670699 0.0006305257 0.0025334280 0.0023480935 0.0007527267 0.0002699272 0.0004699626
#[9] 0.0004869054 0.0025651333 0.0010415805 0.0008574920 0.0010664869 0.0017700297 0.0001666226 0.0008405391
#[17] 0.0017682167 0.0009854172 0.0015350466 0.0017714466 0.0012994945 0.0016563132 0.0003276584 0.0020499529
#[25] 0.0006582073

Но я нахожу эту грубую силу

slope <- terrain(r, "slope")
extract(slope, coords)

быстрее, в 9 раз быстрее, чем моя первая альтернатива, и в 4 раза быстрее, чем вторая альтернатива