vrtility

The goal of vrtility is to make the best use of GDAL’s VRT capabilities for efficient processing of large raster datasets. This package’s primary focus is on the use of GDAL VRT pixel functions using numba. for very efficient raster processing. At present the only function provided is a simple median. There is definitely scope to expand this functionality! Comments and contributions are most welcome!

Features

• No intermediate downloads - the use of nested VRTs enables the download and processing of only the required data in a single gdalwarp call.

• use of numba in python pixel function(s)

TO DO:

• Add efficient masking.
• Add more pixel functions (geometric median in particular).
• clean things up a lot!
• time series functions…

…

Installation

You can install the development version of vrtility from GitHub with:

# install.packages("pak")
pak::pkg_install("Permian-Global-Research/vrtility")
# next set up the required Python environment
vrtility::build_vrtility_python()

Possible extra steps.

For now, {vrtility} requires the use of a virtual python environment - this is a lighter and safer way to manage python envs than with conda or the system python, respectively. It may be necessary to install the virtualenv package for python3.

This can be done on Debian/Ubuntu with:

sudo apt-get install python3-venv

or if on MacOS or Windows, you can install it with pip:

python3 -m pip install virtualenv

Example

Here is a simple example where we: define a bounding box, search a STAC catalog for Sentinel-2 data, create a VRT, and then create a composite image using a median pixel function.

library(vrtility)
library(tictoc)

bbox <- gdalraster::bbox_from_wkt(
  wkt = wk::wkt("POINT (144.3 -7.6)"),
  extend_x = 0.17,
  extend_y = 0.125
)

trs <- to_generic_projected(bbox)

te <- gdalraster::bbox_transform(
  bbox,
  gdalraster::srs_to_wkt("EPSG:4326"),
  trs
)

s2_stac <- sentinel2_stac_query(
  bbox = bbox,
  start_date = "2023-01-01",
  end_date = "2023-12-31",
  max_cloud_cover = 30,
  assets = c(
    "B02",
    "B03",
    "B04",
    "SCL"
  )
)
tic()
median_composite <- vrt_collect(s2_stac) |>
  vrt_set_maskfun(
    "SCL",
    valid_bits = c(4, 5, 6, 7, 11)
  ) |>
  vrt_stack() |>
  vrt_set_pixelfun() |>
  vrt_warp(
    outfile = fs::file_temp(ext = "tif"),
    t_srs = trs,
    te = te
  )
#> 0...10...20...30...40...50...60...70...80...90...100 - done.

toc()
#> 45.749 sec elapsed

plot_raster_src(
  median_composite,
  c(3, 2, 1),
  minmax_pct_cut = c(1, 88)
)