predict function writes only to a single file in case of multi-layer output

[rvalavi]

Thanks again for creating such a great package.

When I use terra::predict function to write multiple layers, all layers are written to a single file, even when multiple output names are provided. In contrast, terra::writeRaster works for multiple layers and filenames. Is this behavior intentional?

I'm using the predict function for a custom C++ model on very large files, and I want to avoid double writing by directly writing the output layers separately to disk.

Here is an example:

library(terra)
#> terra 1.7.78

r <- terra::rast(nrow = 10, ncol = 10)
values(r) <- rnorm(100)

fake_model <- list()

preds <- terra::predict(
    object = r,
    model = fake_model,
    fun = function(obj, dat, ...) {
        return(
            data.frame(
                l1 = rep(0, nrow(dat)),
                l2 = rep(1, nrow(dat))
            )
        )
    },
    filename = c("outputs/prlyr_1.tif", "outputs/prlyr_2.tif"),
    wopt = list(names = c("LYR1", "LYR2"))
)
#> Warning: PROJ: proj_create_from_name: Cannot find proj.db (GDAL error 1)

print(preds)
#> class       : SpatRaster 
#> dimensions  : 10, 10, 2  (nrow, ncol, nlyr)
#> resolution  : 36, 18  (x, y)
#> extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (CRS84) 
#> source      : prlyr_1.tif 
#> names       : LYR1, LYR2 
#> min values  :    0,    1 
#> max values  :    0,    1

dir("outputs/")
#> [1] "prlyr_1.tif"

# now with write raster
terra::writeRaster(
    x = preds, 
    filename = c("outputs/wlyr_1.tif", "outputs/wlyr_2.tif")
)
#> Warning: PROJ: proj_create_from_name: Cannot find proj.db (GDAL error 1)
#> Warning: PROJ: proj_create_from_name: Cannot find proj.db (GDAL error 1)

dir("outputs/")
#> [1] "prlyr_1.tif" "wlyr_1.tif"  "wlyr_2.tif"

sessionInfo()
#> R version 4.4.0 (2024-04-24 ucrt)
#> Platform: x86_64-w64-mingw32/x64
#> Running under: Windows Server 2016 x64 (build 14393)
#> 
#> Matrix products: default
#> 
#> 
#> locale:
#> [1] LC_COLLATE=English_Australia.1252  LC_CTYPE=English_Australia.1252   
#> [3] LC_MONETARY=English_Australia.1252 LC_NUMERIC=C                      
#> [5] LC_TIME=English_Australia.1252    
#> 
#> time zone: Australia/Sydney
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] terra_1.7-78
#> 
#> loaded via a namespace (and not attached):
#>  [1] digest_0.6.35     codetools_0.2-20  fastmap_1.2.0     xfun_0.44        
#>  [5] glue_1.7.0        knitr_1.46        htmltools_0.5.8.1 rmarkdown_2.27   
#>  [9] lifecycle_1.0.4   cli_3.6.2         reprex_2.1.0      withr_3.0.0      
#> [13] compiler_4.4.0    rstudioapi_0.16.0 tools_4.4.0       evaluate_0.23    
#> [17] Rcpp_1.0.12       yaml_2.3.8        rlang_1.1.3       fs_1.6.4

^{Created on 2024-07-23 with reprex v2.1.0}

[rhijmans]

Hi Roozbeh! The way it currently works is that you only can get one file from any terra method, except for writeRaster.
That could be changed, but it would take a considerable amount of work, and could be quite disruptive. Given the long list of issues right now, I do not think it is going to happen anytime soon. I also think that in many cases multilayer files are easier to keep track of, but to each their own!

[rvalavi]

Hi Robert, thanks for your reply. That makes sense. I just wanted to confirm if this is expected behavior or if I might be missing something. Honestly, terra is so fast that passing the predict output to writeRaster is hardly noticeable.

I'll leave it to you to decide whether to close this issue.