Search code examples
rconditional-statementsrasterterra

Assign values to raster from another raster based on condition with holes

I have three raster layers that I would to update based on certain conditions.

  • lc : landcover
  • b : input snowcover
  • c : modified snowcover

When the landcover is "anthropogenic", I would like to able to replace the values from the modified snowcover using the values from the input snowcover. I have been able to make this work as long as there are values in every single cell of all three raster files.

c[lc=='anthropogenic']<-b[lc=='anthropogenic']

However, when the raster includes cells without values (at the edges of the study area), it throws an error even though the extents match perfectly.

c[lc=='anthropogenic']<-b[lc=='anthropogenic']
#Error: [`[<-`] length of cells and values do not match

Here is some dummy code (rather silly) which reproduces the issue.

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

enter image description here

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

Related questions but do not deal with the specific issue: Conditional update without holes Conditionally update rast values from another raster using terra Compares a raster with a polygon but it needs to be all rasters for my actual analysis Extracting values from one raster based on condition in another raster and distinguished by polygons in a vector file

Solution

  • Minimal example data

    library(terra)
lc  <- rast(ncol=4, nrow=5)
lc <- init(lc, "col")
snow <- init(lc, "row")
modsnow <- init(lc, "cell")

    You can use terra::ifel

    d <- ifel(lc == 2, modsnow, snow)

    Under the hood, this uses a combination of terra::cover, terra::mask and terra::classify

    You can also use terra::lapp

    f <- function(x, y, z) {
   i <- which(x==2)
   y[i] <- z[i]
   y
}

x <- lapp(c(lc, snow, modsnow), f)