我有三个栅格图层,我想根据某些条件更新它们。
当土地覆盖是“人为”时,我希望能够使用输入积雪中的值替换修改后积雪中的值。只要所有三个光栅文件的每个单元格中都有值,我就能够完成这项工作。
c[lc=='anthropogenic']<-b[lc=='anthropogenic']
但是,当栅格包含没有值的像元(位于研究区域的边缘)时,即使范围完全匹配,它也会抛出错误。
c[lc=='anthropogenic']<-b[lc=='anthropogenic']
#Error: [`[<-`] length of cells and values do not match
这里有一些虚拟代码(相当愚蠢),可以重现该问题。
library(terra)
library(geodata)
#import world countries dataset
w <- world(path=".", resolution = 1)
#select three countries on coastline
sel<-w[w$NAME_0%in%c('Ghana', 'Togo', 'Benin'),]
#add dummy variable that is numeric (class would be equivalent)
sel$num<-c(1, 2, 3)
#import soil dataset
soil_full<-soil_af(var="pH", depth=15, path=".")
#crop to extent of three countries
soil<-crop(soil_full, ext(sel))
soil_dummy<-soil
#rasterize country dataset to same extent and resolution as soil dataset
values(soil_dummy)<-NA
sel_rast<-rasterize(x=sel, y=soil_dummy, field="num",
background=NA, update=TRUE, touches=TRUE)
names(sel_rast)<-"num"
par(mfrow=c(1,2))
plot(sel_rast)
plot(soil)
sel_rast[sel_rast==1]<-soil[sel_rast==1]
#Error: [`[<-`] length of cells and values do not match
sel_rast
#class : SpatRaster
#dimensions : 921, 853, 1 (nrow, ncol, nlyr)
#resolution : 0.008333333, 0.008333333 (x, y)
#extent : -3.258333, 3.85, 4.741667, 12.41667 (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 (EPSG:4326)
#source(s) : memory
#name : num
#min value : 1
#max value : 3
soil
#class : SpatRaster
#dimensions : 921, 853, 1 (nrow, ncol, nlyr)
#resolution : 0.008333333, 0.008333333 (x, y)
#extent : -3.258333, 3.85, 4.741667, 12.41667 (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 (EPSG:4326)
#source(s) : memory
#name : pH_5-15cm
#min value : 4.68
#max value : 7.85
#attempt with a mask
test<-mask(soil, sel_rast)
sel_rast[sel_rast==1]<-test[sel_rast==1]#fails
#crop to a subset within the landmass
sel_rast.c<-crop(sel_rast, ext(sel_rast)/10)
soil.c<-crop(soil, ext(sel_rast)/10)
plot(soil.c)
plot(sel_rast.c)
sel_rast.c[sel_rast.c==3]<-soil.c[sel_rast.c==3]#works
plot(sel_rast.c)
相关问题但不涉及具体问题: 有条件更新无漏洞 使用 terra 有条件地更新另一个栅格的栅格值 将栅格与多边形进行比较,但对于我的实际分析来说,它需要全部是栅格 根据另一栅格中的条件从一个栅格中提取值并通过矢量文件中的多边形进行区分
最小示例数据
library(terra)
lc <- rast(ncol=4, nrow=5)
lc <- init(lc, "col")
snow <- init(lc, "row")
modsnow <- init(lc, "cell")
您可以使用
terra::ifel
d <- ifel(lc == 2, modsnow, snow)
在底层,它使用
terra::cover
、terra::mask
和 terra::classify
的组合
您还可以使用
terra::lapp
f <- function(x, y, z) {
i <- which(x==2)
y[i] <- z[i]
y
}
x <- lapp(c(lc, snow, modsnow), f)