How to subset (classify ) raster based on another raster grid cells values?

Question

How to reclassify (subset) a raster r1 ( of the same dimension and extent as r2) based on following conditions in r2 in the given example.

Conditions:

• If the grid cells values of r2 are >0.5, retain corresponding value and adjacent 2 grids cells next to 0.5 values (ie buffer the grid cells with values >0.5 in r2 to the surrounding two grids in all directions) in r1 and change other values to 0.

ie. how can i change grid cells values in r1 such that it gives those values which correspond to >0.5 value grid cells and its buffering (surrounding) two grid cells in each direction in r2.

If I only had to get grid cells >0.5 I would have easily obtained by the following code, however, I want to extract the >0.5 value as well as the value of the 2 surrounding gridcells as well.

Sample example calculation code is :

set.seed(150)
r1 <- raster(ncol=10, nrow=5) #Create rasters
values(r1) = round(runif(ncell(r1),5,25))
r2 <- raster(ncol=10, nrow=5)
values(r2) = round(runif(ncell(r2),0.1,1))

selfun <- function(x,y) {
  ifelse( x >0.5, y,0)
}  # It works only for >0.5 gridcells, i need this gridcells and its adjacent
  #two gridcells in each direction.
# it's like buffering the >0.5 grid cells with adjacent two grids and retaining corresponding grid cells values.

r3<-overlay(r2,r1,fun=selfun)
plot(r3)

Thank you.

Answer 1

We can use the focal function to create a mask showing the pixels of interest, and use the mask function to retrieve the values.

I am going to create my examples because the example rasters you created are too small for demonstration.

# Create example raster r1
set.seed(150)
r1 <- raster(ncol = 100, nrow = 50) 
values(r1) <- round(runif(ncell(r1), 5, 25))

r1 <- crop(r1, extent(-60, 60, -30, 30))

plot(r1)

# Create example raster r2
r2 <- raster(ncol = 100, nrow = 50)
values(r2) <- rnorm(ncell(r2), mean = -1)

r2 <- crop(r2, extent(-60, 60, -30, 30))

plot(r2)

The first step is to process r2 by replacing any values larger than 0.5 to be 1 and others to be NA.

# Replace values > 0.5 to be 1, others to be NA
r2[r2 > 0.5] <- 1
r2[r2 <= 0.5] <- NA

plot(r2)

Before using the focal function, we need to define a matrix representing the window and a function.

# Define a matrix as the window
w <-  matrix(c(NA, NA, 1, NA, NA, 
               NA, 1, 1, 1, NA,
               1, 1, 1, 1, 1, 
               NA, 1, 1, 1, NA,
               NA, NA, 1, NA, NA), ncol = 5)
# Define a function to populate the maximum values when w = 1
max_fun <- function(x, na.rm = TRUE) if (all(is.na(x))) NA else max(x, na.rm = na.rm)

We can then apply the focal function.

r3 <- focal(r2, w = w, fun = max_fun, pad = TRUE)

plot(r3)

r3 is a layer showing the pixels we want values from r1. We can now use the mask function for this.

# Use the mask function extract values in r1 based on r3
r4 <- mask(r1, r3)
# Replace NA with 0
r4[is.na(r4)] <- 0

plot(r4)

r4 is the final output.