满足特定条件时用另一个栅格替换栅格值

Question

我有三个光栅。

library(raster)

a <- raster(ncol=100, nrow=100)
set.seed(2)     
values(a) = runif(10000, min=-1, max=1)    # define the range between -1 and 1

b <- raster(ncol=100, nrow=100)
set.seed(2)     
values(b) = runif(10000, min=0.97, max=0.99)

c <- raster(ncol=100, nrow=100)
set.seed(2)     
values(c) = runif(10000, min=0, max=1)

现在我想根据多个条件创建一个新的栅格。条件如

当 a < 0.2 时，新光栅像素应具有光栅 c 的值；

当 a > 0.5 时，新光栅像素的值应为 1.85；

当 0.2 <= a <= 0.5 时，新光栅像素应具有光栅 b 的值。我正在使用以下代码来实现

# create empty copy
E <- raster(a)

E[(a < 0.2)] <- NA
E <- cover(E, c)
E[(a >= 0.2) & (a <= 0.5)] <- NA
E <- cover(E, b)
E[a > 0.5] <- 1.85

这样做是否正确？或者有更好的替代方法。

Answer 1

我会用一些 ifelse 替换栅格中的值来做到这一点：

从空白光栅开始：

> F = raster(a)

然后在 a<0.2 中取值 c 否则使用 NA:

> F[] = ifelse(a[]<0.2, c[], NA)

然后 a 超过 0.5 时将 F 设置为 1.85，否则使用 F:

中的任何值

> F[] = ifelse(a[]>0.5, 1.85, F[])

如果 a 介于 0.2 和 0.5 之间，则取 b 中的值，否则使用 F:

中的值

> F[] = ifelse(a[]>=0.2 & a[]<=0.5, b[], F[])

这与您的 E:

给出的值相同

> all(F[]==E[])
[1] TRUE

或者您可以通过栅格值向量中的条件替换来实现：

> F=raster(a)
> F[a[]<0.2] <- c[a[]<0.2]
> F[a[]>0.5] <- 1.85
> F[a[]>=0.2 & a[]<=0.5] <- b[a[]>=0.2 & a[]<=0.5]
> all(F[]==E[])
[1] TRUE

我不确定在速度、可读性或灵活性方面哪个更好。无论如何，我会将这一切包装成一个函数：

replacer = function(a, b, c,
  low_thresh=0.2, high_thresh=0.5,
  high_value=1.85){
...
}

如果速度是一个问题，然后对其进行基准测试。否则，请使用您认为看起来最好并且最容易维护的东西。

这是三种方法（我已经优化了第三种以避免重复测试）：

replacer_1 = function(a, b, c,
  low_thresh=0.2, high_thresh=0.5,
  high_value=1.85){
    E <- raster(a)
    E[(a < low_thresh)] <- NA
    E <- cover(E, c)
    E[(a >= low_thresh) & (a <= high_thresh)] <- NA
    E <- cover(E, b)
    E[a > high_thresh] <- high_value
    return(E)
}

replacer_2 = function(a, b, c,
  low_thresh=0.2, high_thresh=0.5,
  high_value=1.85){
    F = raster(a)
    F[] = ifelse(a[]<0.2, c[], NA)
    F[] = ifelse(a[]>0.5, 1.85, F[])    
    F[] = ifelse(a[]>=0.2 & a[]<=0.5, b[], F[])
    return(F)
}

replacer_3 = function(a, b, c,
  low_thresh=0.2, high_thresh=0.5,
  high_value=1.85){
    low = a[]<low_thresh
    F[low] <- c[low]
    F[a[]>high_thresh] <- high_value
    mid =a[]>=low_thresh & a[]<=high_thresh 
    F[mid] <- b[mid]
    return(F)
}

检查所有 return 相同的答案：

> E1 = replacer_1(a,b,c)
> E2 = replacer_2(a,b,c)
> E3 = replacer_3(a,b,c)
> all(E1[]==E2[])
[1] TRUE
> all(E1[]==E3[])
[1] TRUE
> all(E2[]==E3[])
[1] TRUE

然后用microbenchmark包比较：

> microbenchmark(replacer_1(a,b,c), replacer_2(a,b,c), replacer_3(a,b,c))
Unit: microseconds
                expr       min        lq       mean    median         uq
 replacer_1(a, b, c) 50687.567 52486.878 54089.3265 53721.770 55221.1125
 replacer_2(a, b, c)  2359.977  2507.661  2667.1080  2581.568  2642.9790
 replacer_3(a, b, c)   485.086   504.377   543.6963   542.970   565.7025

这表明第三种方法比第一种方法快大约 100 倍。享受我救你的所有时光吧！

Answer 2

简单示例数据

library(raster)
a <- b <- c <- raster(ncol=10, nrow=10)
set.seed(2)    
values(a) = sort(runif(100, min=-1, max=1))
values(b) = 3
values(c) = 5

到达那里的一种方法是

m <- reclassify(a, c(-Inf, 0.2, 1, 0.2, Inf, NA))
x <- mask(c, m)
m <- reclassify(a, c(-Inf, 0.2, NA, 0.2, 0.5, 1, 0.5, Inf, NA))
y <- mask(b, m)
z <- init(a, 1.85)

现在

r <- merge(x, y, z)

或

r <- merge(x, y)
r <- cover(r, z)

在 ifelse 上跟随@Spacedman 的脚步；你也可以使用隐藏的 .ifel 方法

s <- raster:::.ifel(a < 0.2, c, raster:::.ifel(a > 0.5, 1.85, b))

或者使用 terra::ifel --- 应该也更快

library(terra)
aa <- rast(a)
bb <- rast(b)
cc <- rast(c)

ss <- ifel(aa < 0.2, cc, ifel(aa > 0.5, 1.85, bb))

满足特定条件时用另一个栅格替换栅格值

Replacing raster values with another raster when meeting certain condition

r

r-raster