2016-02-12 6 views
3

У меня есть ряд растровых слоев, содержащихся в одном stack, а также объект SpatialPoints, содержащий координаты точки. Я пытаюсь изменить значение растровых слоев на NA, если ячейка содержит точку.Измените значения растровых ячеек на NA, если ячейка содержит точку

Приведенный ниже пример приведен ниже. Однако, чтобы поместить это в контекст, мои настоящие данные состоят из 10 слоев «среды обитания» (высота, крышка купола и т. Д.), Которые все содержатся в стеке. Площадь rasters составляет примерно 34 х 26 км. Кроме того, у меня есть ~ 10 000 мест GPS от животных.

Итак, для воспроизводимого примера.

library(raster) 
library(sp) 

сделать три rasters и объединить их в stack

set.seed(123) 
r1 <- raster(nrows=10, ncols=10) 
r1 <- setValues(r1, sample(c(1:50), 100, replace = T)) 

r2 <- raster(nrows=10, ncols=10) 
r2 <- setValues(r2, sample(c(40:50), 100, replace = T)) 

r3 <- raster(nrows=10, ncols=10) 
r3 <- setValues(r3, sample(c(50:55), 100, replace = T)) 

Stack <- stack(r1, r2, r3) 
nlayers(Stack) 

Сделать SpatialPoints

Pts <- SpatialPoints(data.frame(x = sample(-175:175, 50, replace = T), 
           y = sample(-75:75, 50, replace = T))) 

Участок один из растров и точки

plot(Stack[[2]]) 
plot(Pts, add = T) 

enter image description here

Теперь, возможно ли изменить значение каждой ячейки, содержащей хотя бы одну точку, на NA? Было бы здорово выполнить эту операцию на stack.

ответ

4

Я хотел бы использовать extract(), который может попросить вернуть номер ячейки растровой ячейки, по которой каждая точка лежит:

ii <- extract(Stack, Pts, cellnumbers=TRUE)[,"cells"] 
Stack[ii] <- NA 

## Check any one of the layers to see that this worked: 
plot(Stack[[2]]) 
plot(Pts, add=TRUE) 

enter image description here

В качестве альтернативы можно использовать rasterize(), хотя это Вероятно (?) будет медленнее для очень больших растров:

ii <- !is.na(rasterize(Pts, Stack)) 
Stack[ii] <- NA 
+0

Прохладный. Очень полезно. Есть ли причина, по которой вы проиндексировали слой '[[1]]' в 'Stack'? Проверяя это, кажется, что вы можете выполнить операцию на весь стек сразу. например 'ii <- unique (extract (Stack, Pts, cellnumbers = TRUE) [," cells "])' 'Stack [ii] <- NA' –

+0

@ B.Davis Нет большой причины, и на самом деле призыв к' unique() 'также не требуется. Я отредактировал сообщение, чтобы упростить соответственно. –

 Смежные вопросы

  • Нет связанных вопросов^_^