Dynamic branching with raster tiles

library(geotargets)
library(targets)
library(terra)
library(dplyr)

Computationally intensive raster operations that work in pixel-wise manner may be handled well with dynamic branching over tiled subsets of the raster. tar_terra_tiles() is a target factory that enables creating these dynamic branches so downstream targets can iterate over them. This is useful when, for example, loading an entire raster into memory and doing computations on it results in out of memory errors.

In order to use tar_terra_tiles(), we need to break a raster into smaller pieces. We can do that by providing extents used by the raster. The concept of extent is important, so let’s unpack that a bit more.

What is an extent?

The extent describes the four points that cover the area of a raster. The extent of a raster, r, is printed in the summary:

#example SpatRaster
f <- system.file("ex/elev.tif", package="terra")
r <- rast(f)
r

But we can get the extent with ext (extent):

r_ext <- ext(r)
r_ext

Which maps onto the four corners of the raster here:

# some plot helpers
rect_extent <- function(x, ...) {
  rect(x[1], x[3], x[2], x[4], ...)
}
plot_extents <- function(x, ...) {
  invisible(lapply(x, rect_extent, border = "hotpink", lwd = 2))
}

extend(r, 5) |> plot()
lines(r_ext, col = "hotpink", lty = 2)
points(r_ext, col = "hotpink", pch = 16)

Some geo-computational operations can be done independently of one another—we want to take advantage of that, and we can facilitate this by breaking the raster into smaller pieces, by creating new extents that describe new subsets of the raster.

We can use this extent information downstream in the analysis to describe how to break up a raster. This is similar to how we might want to chunk up a data frame into groups to distribute to different CPU cores. To help with this, we’ve got some helper functions.

Helper functions to create multiple extents of a raster

geotargets provides three helper functions that take a SpatRaster and output the extents for tiles:

tile_n(),
tile_grid(), and
tile_blocksize()

We will demonstrate these now.

`tile_n()`

We can use tile_n(), which is the simplest of the three. It produces about n tiles in a grid.

r_tile_4 <- tile_n(r, 4)
r_tile_4

plot(r)
plot_extents(r_tile_4)
plot(r)
tile_n(r, 6) |> plot_extents()

`tile_grid()`

For more control, use tile_grid(), which allows specification of the number of rows and columns to split the raster into. Here we are specify that we want three columns and 1 row:

r_grid_3x1 <- tile_grid(r, ncol = 3, nrow = 1)
r_grid_3x1
plot(r)
plot_extents(r_grid_3x1)

plot(r)
tile_grid(r, ncol = 2, nrow = 3) |> plot_extents()

`tile_blocksize()`

The third included helper is tile_blocksize(), which tiles by file block size. The block size is a property of raster files, and is the number of pixels (in the x and y direction) that is read into memory at a time. Tiling by multiples of block size may therefore be more efficient because only one block should need to be loaded to create each tile target. You can find the blocksize with fileBlocksize:

fileBlocksize(r)

This tells us that it reads in the raster in 43x95 pixel sizes.

The tile_blocksize function is similar to tile_grid, except instead of saying how many rows and columns, we specify in units of blocksize.

If we just run tile_blocksize() on r we get the extents of the specified blocksize:

tile_blocksize(r)

Which is the same as specifying blocksize for row and column at unit 1:

r_block_size_1x1 <- tile_blocksize(r, n_blocks_row = 1, n_blocks_col = 1)
r_block_size_1x1
plot(r)
plot_extents(r_block_size_1x1)

Here the block size is the same size for the first two blocks, and then a much more narrow block. This is different to the two other tile methods.

Here the column block size is the full width of the raster.

So we could instead have the blocksize extent be written out to 2 blocks in a row, and 1 block size for the columns:

r_block_size_2x1 <- tile_blocksize(r, n_blocks_row = 2, n_blocks_col = 1)
r_block_size_2x1
plot(r)
plot_extents(r_block_size_2x1)

This only works when the SpatRaster points to a file—in-memory rasters have no inherent block size.

sources(r)
#force into memory
r2 <- r + 0
sources(r2)
#this now errors
tile_blocksize(r2)

Example targets pipeline

When developing a targets pipeline using tar_terra_tiles() with tile_blocksize(), it’s a good idea to figure out how many tiles tile_blocksize() will create before implementing tar_terra_tiles(). We’ll start by making a bigger raster to experiment with using terra::disagg(), (which makes a higher resolution raster by breaking the pixels into smaller pixels), and making multiple layers.

# contents of _targets.R
library(targets)
library(geotargets)
library(terra)
geotargets_option_set(gdal_raster_driver = "COG")
list(
    tar_target(
        raster_file,
        system.file("ex/elev.tif", package="terra"),
        format = "file"
    ),
    tar_terra_rast(
        r,
        disagg(rast(raster_file), fact = 10)
    ),
    # add more layers
    tar_terra_rast(
        r_big,
        c(r, r + 100, r * 10, r / 2),
        memory = "transient"
    )
)

tar_make()
#> terra 1.7.83
#> ▶ dispatched target raster_file
#> ● completed target raster_file [0 seconds, 7.994 kilobytes]
#> ▶ dispatched target r
#> ● completed target r [0.01 seconds, 822.781 kilobytes]
#> ▶ dispatched target r_big
#> ● completed target r_big [0.024 seconds, 4.055 megabytes]
#> ▶ ended pipeline [0.415 seconds]
tar_load(r_big)
tile_blocksize(r_big)
#> [[1]]
#>      xmin      xmax      ymin      ymax 
#>  5.741667  6.168333 49.765000 50.191667 
#> 
#> [[2]]
#>      xmin      xmax      ymin      ymax 
#>  6.168333  6.533333 49.765000 50.191667 
#> 
#> [[3]]
#>      xmin      xmax      ymin      ymax 
#>  5.741667  6.168333 49.441667 49.765000 
#> 
#> [[4]]
#>      xmin      xmax      ymin      ymax 
#>  6.168333  6.533333 49.441667 49.765000

Four tiles is reasonable, so we’ll go with that. Note that we have to ensure the r_big target is not in-memory for tar_terra_tiles(), so we set the targets option memory = "transient". See the targets documentation on memory for details.

The process that happens from here can be thought of as split-apply-combine.

Split the raster into pieces using the tar_terra_tiles() target factory
- This returns tiles whose extents are created by one of the tile functions described above (tile_n(), tile_grid(), or tile_blocksize()), supplying this to tile_fun.
Apply a function to the rasters.
- This can be any function that would work on a raster, in the case below we use the app function from terra, which applies some function to the cells of a raster.
- To do this we use tar_terra_rast() and then supply the pattern = map(tiles), where tiles is the name of the target created with tar_terra_tiles(). You can think of pattern = map(tiles) as saying: “Do the task for each of the tiles we have specified and return them as a list”
Combine the list of rasters together.
- In this case we use tar_terra_rast() and use merge() on the tiles.

# contents of _targets.R
library(targets)
library(geotargets)
library(terra)
geotargets_option_set(gdal_raster_driver = "COG")
tar_option_set(memory = "transient")
list(
    tar_target(
        raster_file,
        system.file("ex/elev.tif", package="terra"),
        format = "file"
    ),
    tar_terra_rast(
        r,
        disagg(rast(raster_file), fact = 10)
    ),
    tar_terra_rast(
        r_big,
        c(r, r + 100, r * 10, r / 2),
        memory = "transient"
    ),
    tar_terra_tiles(
        tiles,
        raster = r_big,
        tile_fun = tile_blocksize,
        description = "split raster into tiles"
    ),
    tar_terra_rast(
        tiles_mean,
        app(tiles, \(x) mean(x, na.rm = TRUE)),
        pattern = map(tiles),
        description = "some computationaly intensive task performed on each tile"
    ),
    tar_terra_rast(
        merged_mean,
        merge(sprc(tiles_mean)),
        description = "merge tiles into a single SpatRaster"
    )
)

tar_make()
#> terra 1.7.83
#> ✔ skipped target raster_file
#> ✔ skipped target r
#> ✔ skipped target r_big
#> ▶ dispatched target tiles_exts
#> ● completed target tiles_exts [0.005 seconds, 153 bytes]
#> ▶ dispatched branch tiles_11882e184aa27102
#> ● completed branch tiles_11882e184aa27102 [0.003 seconds, 982.55 kilobytes]
#> ▶ dispatched branch tiles_e39b1bcab4d45ba0
#> ● completed branch tiles_e39b1bcab4d45ba0 [0.042 seconds, 249.778 kilobytes]
#> ▶ dispatched branch tiles_b14b762418a51f63
#> ● completed branch tiles_b14b762418a51f63 [0.002 seconds, 599.615 kilobytes]
#> ▶ dispatched branch tiles_bfb55de5c880456f
#> ● completed branch tiles_bfb55de5c880456f [0.003 seconds, 483.702 kilobytes]
#> ● completed pattern tiles 
#> ▶ dispatched branch tiles_mean_1bcefdb80d0879e0
#> ● completed branch tiles_mean_1bcefdb80d0879e0 [2.127 seconds, 159.194 kilobytes]
#> ▶ dispatched branch tiles_mean_4ff32c943e05dc36
#> ● completed branch tiles_mean_4ff32c943e05dc36 [1.854 seconds, 49.028 kilobytes]
#> ▶ dispatched branch tiles_mean_17723520c935f474
#> ● completed branch tiles_mean_17723520c935f474 [1.606 seconds, 100.948 kilobytes]
#> ▶ dispatched branch tiles_mean_984998bf0804c2cf
#> ● completed branch tiles_mean_984998bf0804c2cf [1.406 seconds, 87.272 kilobytes]
#> ● completed pattern tiles_mean 
#> ▶ dispatched target merged_mean
#> ● completed target merged_mean [0.016 seconds, 879.912 kilobytes]
#> ▶ ended pipeline [7.543 seconds]

We can see from tar_make() output above and the plots below that tiles and tiles_mean are both patterns with four branches each.

library(terra)
#> terra 1.7.83
tar_load(tiles_mean)
op <- par(mfrow = c(2,2))
for (i in seq_along(tiles_mean)) {
    plot(tiles_mean[[i]])
}

par(op)

And combined, they make the full plot again.

plot(tar_read(merged_mean))