如何根据另一个栅格的网格单元值对栅格进行子集(分类)?

4
如何根据给定示例中r2中的以下条件重新分类(子集)与r2相同维度和范围的栅格r1
条件:
- 如果r2的网格单元值为>0.5,则保留相应的值并在r1中将邻近的2个网格单元保持在0.5值之上(即将r2中数值大于0.5的网格单元缓冲到所有方向的周围两个网格单元),并将其他值更改为0。
即:如何更改r1中的网格单元值以使其给出与r2中大于0.5的值对应的值以及每个方向中的其两个缓冲(周围)网格单元的值。
如果我只需要获取大于0.5的网格单元,我可以通过以下代码轻松获得,但是我想提取大于0.5的值以及相邻两个网格单元的值。
样例计算代码如下:
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.

1个回答

7
我们可以使用focal函数创建一个掩码,显示感兴趣的像素,并使用mask函数检索值。
我将创建我的示例,因为您创建的示例栅格太小,无法进行演示。
# 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)

enter image description here

# 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)

enter image description here

第一步是通过将大于0.5的任何值替换为1并将其余值替换为NA来处理r2

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

plot(r2)

enter image description here

在使用focal函数之前,我们需要定义一个代表窗口和一个函数的矩阵。
# 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)

我们可以应用focal函数。
r3 <- focal(r2, w = w, fun = max_fun, pad = TRUE)

plot(r3)

enter image description here

r3是显示我们想要从r1获取值的像素层。现在我们可以使用mask函数来实现。

# 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是最终输出结果。

enter image description here


谢谢,但是你是如何定义窗口矩阵的? - Lily Nature
您可以打印矩阵w并查看其结构。您说您想要相邻的两个网格单元。我认为这可能意味着所有方向(上,下,左,右),但我不确定对角线的条件。我的示例矩阵假定对角线有一个单元格,但您可以根据需要修改矩阵。 - www

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接