Is there a raster value modification function faster than terra::subst?

Question

I have 280 large (400MB) raster files in TIFF format. I am currently using terra::subst to modify the values of each raster file.

For example:

raster_matrix<-matrix(rep(101:200,10000*1600),nrow=1000*40,ncol=1000*40)
raster<-terra::rast(raster_matrix)
raster_subst <- terra::subst(raster, c(151, 152, 153), NA))
raster_subst <- terra::subst(raster_subst, c(160, 161, 162), 1))
raster_subst <- terra::subst(raster_subst, c(170, 171, 172), 2))

However, this process takes too long when applied to all 280 files. Is there a faster and more efficient way to achieve this?

Answer 1

This takes 1/3 of the time on the example dataset below. I could not recreate a SpatRaster as big as your example on my machine so am unsure how it will scale. Hopefully someone more knowledgeable will know how to reduce the time even further.

library(terra)
library(microbenchmark)

raster_matrix<-matrix(rep(101:200,5000*1600),nrow=500*40,ncol=500*40)
raster <- terra::rast(raster_matrix)

results <- microbenchmark(
  separate <- {raster_subst <- terra::subst(raster, c(151, 152, 153), NA);
  raster_subst <- terra::subst(raster_subst, c(160, 161, 162), 1);
  raster_subst <- terra::subst(raster_subst, c(170, 171, 172), 2)},
  together <- subst(raster, 
                        c(151, 152, 153, 160, 161, 162,170, 171, 172), 
                        rep(c(NA, 1, 2), each = 3)),
  times = 10)

results
Unit: seconds
              min       lq     mean   median       uq      max neval
separate 41.73862 41.91065 42.43262 42.24467 43.16119 43.41679    10
together 13.56677 13.82499 13.89107 13.90967 14.02498 14.18583    10

Answer 2

If you want to change the values for ranges (e.g. 150 to 152), terra::classify may be faster. Note argument right=NA to have the interval open at the left and right.

Illustrated by building on l-tyrone's benchmarking.

library(terra)
library(microbenchmark)

raster_matrix <- matrix(101:200,nrow=500,ncol=500)
raster <- terra::rast(raster_matrix)

m1 <- cbind(c(151, 152, 153, 160, 161, 162,170, 171, 172), 
                        rep(c(NA, 1, 2), each = 3))

m2 <- matrix(c(151, 153, NA, 160, 162, 1, 170, 172, 2), ncol=3, byrow=TRUE)
 

microbenchmark(
  separate = {terra::subst(raster, c(151, 152, 153), NA);
        raster_subst <- terra::subst(raster_subst, c(160, 161, 162), 1);
        raster_subst <- terra::subst(raster_subst, c(170, 171, 172), 2)},
  together = subst(raster, m1[,1], m1[,2]),
  reclass1 = classify(raster, m1),  
  reclass2 = classify(raster, m2, right=NA),  
  times = 10
)

#Unit: milliseconds
#     expr     min      lq     mean   median      uq     max neval cld
# separate 14.9612 15.0116 15.66967 15.22865 15.5153 19.4597    10 a  
# together  5.8914  6.0275  6.26044  6.20995  6.4754  6.7361    10  b 
# reclass1  5.9294  6.1203  6.24878  6.21910  6.4330  6.5321    10  b 
# reclass2  4.9246  4.9876  5.16236  5.03045  5.3546  5.6490    10   c