5

I have a number of raster layers contained in a single stack as well as a SpatialPoints object containing point coordinates. I am trying to change the value of the raster layers to NA if the cell contains a point.

I have provided a reproducible example below. However to put this in context, my real data consists of a 10 layers of 'habitat' (elevation, canopy cover, etc) that are all contained in a stack. The rasters cover an area that is roughly 34 x 26 km. In addition I have ~10,000 GPS locations from animals.

So, for the reproducible example.

library(raster)
library(sp)

Make three rasters and combine them into a 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)

Make the SpatialPoints

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

Plot one of the rasters and the points

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

enter image description here

Now, is it possible to change the value of every cell that contains at least one point to NA? It would be great to perform this operation on the stack.

B. Davis
  • 3,391
  • 5
  • 42
  • 78

1 Answers1

6

I'd use extract(), which can be asked to return the cell number of the raster cell over which each point lies:

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

Alternatively, one could use rasterize(), though it'll likely (?) be slower for very large rasters:

ii <- !is.na(rasterize(Pts, Stack))
Stack[ii] <- NA
Josh O'Brien
  • 159,210
  • 26
  • 366
  • 455
  • Cool. Very helpful. Is there a reason you indexed layer `[[1]]` in the `Stack`? Testing this out, it seems like you could perform the operation on the entire stack at once. e.g. `ii <- unique(extract(Stack, Pts, cellnumbers=TRUE)[,"cells"])` `Stack[ii] <- NA` – B. Davis Feb 12 '16 at 23:41
  • @B.Davis No great reason, and in fact the call to `unique()` isn't needed either. I've edited the post to simplify accordingly. – Josh O'Brien Feb 12 '16 at 23:48