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

编程问答 2022-06-03

问题描述

我有三个栅格。

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

这是否是正确的方法？还是有更好的替代方法。

解决方法

我可以通过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,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,high_value=1.85){
    F = raster(a)
    F[] = ifelse(a[]<0.2,NA)
    F[] = ifelse(a[]>0.5,F[])    
    F[] = ifelse(a[]>=0.2 & a[]<=0.5,F[])
    return(F)
}

replacer_3 = function(a,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,c)
> E2 = replacer_2(a,c)
> E3 = replacer_3(a,c)
> all(E1[]==E2[])
[1] TRUE
> all(E1[]==E3[])
[1] TRUE
> all(E2[]==E3[])
[1] TRUE

然后使用microbenchmark包进行比较：

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

这告诉我第三种方法的速度大约是第一种方法的100倍。享受我拯救您的所有时间！

,

请尽量简化示例数据，以便更轻松地跟踪发生的情况

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,Inf,NA))
x <- mask(c,m)
m <- reclassify(a,NA,0.5,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,raster:::.ifel(a > 0.5,b))

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

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

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