terra targets integration

[hansvancalster]

Thanks for creating targets 👍

The reprex below shows how terra::wrap() and terra::rast() conversion works to avoid the "non-exportable object" problem (terra::rast() is an external pointer to a raster file on disk). The same approach also works for terra::vect().

This approach works, but seems wasteful if the raster tif file is very large (often several Gb) because a copy is made and stored. Do alternative approaches exist that do not have this drawback?

library(targets)
#> Warning: package 'targets' was built under R version 4.1.3
tar_script({
    library(terra)
    example_tif_file <- function() {
        system.file("ex/elev.tif", package = "terra")
    }
    create_raster <- function(rasterfile) {
        rast(rasterfile)
    }
    create_packed_raster <- function(rr) {
        wrap(rr)
    }
    list(tarchetypes::tar_file(example_tif, example_tif_file()), 
        tar_target(raster, create_raster(rasterfile = example_tif), 
            format = "qs"), tar_target(packed_raster, create_packed_raster(raster), 
            format = "qs"))
})
library(terra)
#> Warning: package 'terra' was built under R version 4.1.2
#> terra 1.5.21
tar_visnetwork()
#> terra 1.5.21
#> Warning messages:
#> 1: package 'targets' was built under R version 4.1.3 
#> 2: package 'terra' was built under R version 4.1.2 

tar_make()
#> terra 1.5.21
#> * start target example_tif
#> * built target example_tif
#> * start target raster
#> * built target raster
#> * start target packed_raster
#> * built target packed_raster
#> * end pipeline: 1.27 seconds
#> Warning messages:
#> 1: package 'targets' was built under R version 4.1.3 
#> 2: package 'terra' was built under R version 4.1.2 
tar_read(raster)
#> class       : SpatRaster
#> Error in x@ptr$nrow(): external pointer is not valid
packed_raster <- tar_read(packed_raster)
rast(packed_raster)
#> class       : SpatRaster 
#> dimensions  : 90, 95, 1  (nrow, ncol, nlyr)
#> resolution  : 0.008333333, 0.008333333  (x, y)
#> extent      : 5.741667, 6.533333, 49.44167, 50.19167  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : memory 
#> name        : elevation 
#> min value   :       141 
#> max value   :       547

^{Created on 2022-03-16 by the reprex package (v2.0.1.9000)}

[wlandau]

If the large tif file is is an output image, I would suggest writing it once to disk and tracking it with format = "file": https://books.ropensci.org/targets/data.html#external-files. If you use it as input data, you can still treat it as an input file (still format = "file", but only load the subset you need in the target that does the analysis. Maybe this involves writing a function that accepts the file path, creates a terra pointer to it, only reads in the parts of the file it needs, and returns a small object. Otherwise, it's hard to get around copying the whole file, but can be made easier with a custom tar_format() format. The discussion at #805 is similar, the OP has a gist that does a similar thing with parquet (the file is copied, but the object is a pointer where possible).

[mschmidty]

@brownag Your example worked for me initially when I tried the reprex. But then when I tried it with a .tif that I read in from a file, cropped and then tried to write with the same format, the pipeline runs but then when I try to load it I get an error:

r$> tar_load_everything()
Error: [rast] file does not exist: C:/Users/MichaelSchmidt2/AppData/Local/Temp/RtmpquvflJ/spat_5d64ae247b2_23908.tif
In addition: Warning message:
C:/Users/MichaelSchmidt2/AppData/Local/Temp/RtmpquvflJ/spat_5d64ae247b2_23908.tif: No such file or directory (GDAL error 4) 

[brownag]

@mschmidty I think this is because in first example the raster exists/fits in memory, whereas in this second case you are wrapping a SpatRaster backed by a temporary .tif file. This is a related issue on terra repo: rspatial/terra#1150

One option would be to force the loading of the file into memory if possible (with set.values(r))--but that may not scale well.

What you may need instead is terra::wrapCache() which allows for specifying a custom folder to cache the source files for a SpatRaster before wrapping it.

[brownag]

@mschmidty

The following feels really hacky, but there might be some kernel of value in it.

Conforming with the tar_format() write argument requirement that: The function need not return a value, but the file written to path must be a single file, and it cannot be a directory. I figured that perhaps you could bundle the .tif and .tif.aux.xml in a ZIP file (or some other similar single file type).

The following seems to work and perhaps overhead due to compressing/decompressing larger files could be avoided if no compression was used (i.e. utils::zip(..., flags="-r0X"))

# ...

# replace my_files with directory for file-based targets
dir.create("my_files", showWarnings = FALSE, recursive = TRUE)

format_cat_geotiff <- tar_format(
  read = function(path) {
    # replace my_files with directory for file-based targets
    td <- file.path("my_files", basename(path))
    utils::unzip(path, exdir = td)
    terra::rast(list.files(td, pattern = paste0("my_", basename(path), "_file.tif$"),
                           full.names = TRUE, recursive = TRUE))
  },
  write = function(object, path) {
    # bundle .tif and .tif.aux.xml into single file
    utils::zip(path, files = paste0(terra::sources(
      terra::writeRaster(
        x = object,
        # replace my_files with directory for file-based targets
        filename = paste0("my_files/my_", basename(path), "_file.tif"),
        filetype = "GTiff",
        overwrite = TRUE
      )
    ), c("", ".aux.xml")))
    file.rename(paste0(path, ".zip"), path)
  },
  marshal = function(object) terra::wrap(object),
  unmarshal = function(object) terra::unwrap(object)
)

Gives:

> tar_make()
▶ dispatched target r
 (deflated 43%)pace/terra-targets-test/my_files/my_r_file.tif
 (deflated 73%)pace/terra-targets-test/my_files/my_r_file.tif.aux.xml
● completed target r [2.06 seconds]
▶ completed pipeline [3.16 seconds]
 
> tar_load_everything()
> r
class       : SpatRaster 
dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
resolution  : 36, 18  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : my_r_file.tif 
categories  : cover, letters, value 
name        :  cover 
min value   : forest 
max value   :  urban 

[mdsumner]

ah nice, you might consider a format that does support category data without sidecars, I haven't tried it for the use-case but it would be something like

writeRaster(r, "file.img", filetype = "HFA")

that could produce a sidecare xml for the statistics, but the RAT is internal and the .xml is not needed, for completeness could do

setGDALconfig("GDAL_PAM_ENABLED", "NO")
writeRaster(r, "file2.img", filetype = "HFA")

but you should get and reset the original config value (and HFA may not be a suitable format for other reasons), but I'll leave that as an option

[brownag]

+1 to just using a different file format if possible.

FWIW, I realized that with the ZIP file approach the read step I proposed above could be simplified if you just read from the target with GDAL virtual file system /vsizip/ path. So no unzipping or file renaming needed.

#...

format_cat_geotiff <- tar_format(
  read = function(path) {
    terra::rast(paste0("/vsizip/", file.path(path, basename(path))))
  },
  write = function(object, path) {
    # bundle .tif and .tif.aux.xml into single file
    src <- paste0(terra::sources(
      terra::writeRaster(
        x = object,
        filename = basename(path),
        filetype = "GTiff",
        overwrite = TRUE
      )
    ), c("", ".aux.xml"))
    utils::zip(path, files = basename(src))
    unlink(src)
  },
  marshal = function(object) terra::wrap(object),
  unmarshal = function(object) terra::unwrap(object)
)

# Replace the target list below with your own:
list(
  tar_target(
    name = r.zip,
    raster_test(),
    format = format_cat_geotiff 
  )
)

The above cleans up the .tif file initially written by writeRaster() after zipping, so you don't have two copies lying around. The trick is that you need to have your target name end in ".zip" otherwise GDAL will not properly recognize the path. (i.e. /vsizip/_targets/objects/r/r' does not exist in the file system, and is not recognized as a supported dataset name. (GDAL error 4))