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

【问题标题】：Replacing raster values with another raster when meeting certain condition满足特定条件时用另一个栅格替换栅格值
【发布时间】：2022-01-05 18:23:32
【问题描述】：

我有三个栅格。

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

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

当 0.2

# 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

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

【问题讨论】：

你的代码能做到你想要的吗？您是在问这是否有效，或者这是否是最好的方法？如果您使栅格更小（小得多），则可以更容易地判断您的代码是否正常工作，因此您可以通过查看输出栅格是否符合预期来判断。也许从 4x4 栅格开始？
在您的文本中您说“当 a > 0.5 时，新的光栅像素应该具有 1.85 的值”但您的代码看起来像 E[a > 0.5] <- 0.99 - 0.99 应该是 1.85 吗？
在你的文本中你说“当 0.2 E <- cover(E, b) 那么您的意思是“光栅 b 中的值”吗？
抱歉那些错误，我已经纠正了这个问题。有没有更好的方法来做到这一点？

标签： r r-raster

【解决方案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)
}

检查他们都返回相同的答案：

> 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 倍。享受我为你拯救的所有时间！

【讨论】：

速度有点令人担忧，因为我正在为大型光栅做这件事。能否请您提供基准代码？

【解决方案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)

跟随@Spacedman 在ifelse 上的领先；你也可以使用隐藏的 .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))

【讨论】：