GDAL API Bindings for R gdalraster logo

https://firelab.github.io/gdalraster/

 

Chris Toney (@ctoney)
USDA Forest Service, Rocky Mountain Research Station

Michael D. Sumner
Australian Antarctic Division

 

Logo designed by Sylvain Beorchia

R package gdalraster gdalraster logo

osgeo-like bindings + custom tools

  • GDAL Raster and Vector Data Models
  • Modularized, low-level access to GDAL objects
  • Thin layer of abstraction over Driver, Dataset, RasterBand, OGRLayer, GDALAlgorithm, …
    → Rcpp modules w/ C++ classes exposed in R
  • Safe object lifetimes
  • ℹ️ Extensive documentation built-in and online

Comprehensive API bindings gdalraster logo

→ Raster and Vector APIs, plus:

  • GDAL’s modernized CLI framework!
  • GDAL traditional utilities
  • Spatial reference systems and transformation (GDAL PROJ)
  • Geometry API operating on WKB/WKT
  • Virtual Systems Interface (VSI)
  • Partial support for multidimensional raster
  • Custom tools 🔧

Installation gdalraster logo

CRAN (incl. Windows and macOS binaries)

install.packages("gdalraster")

conda-forge (release builds)

conda install --channel conda-forge r-gdalraster

GitHub

remotes::install_github("firelab/gdalraster")

R-universe (dev builds)

install.packages("gdalraster", repos = c("https://firelab.r-universe.dev", "https://cran.r-project.org"))

Short overview/highlights of each

Raster

Vector

Spatial ref and transformation

Geometry API

Virtual file systems

Bindings to the GDAL CLI framework

Raster Data Model

Class GDALRaster

Holds a pointer to a raster Dataset

Constructor w/ optional args:

ds <- new(GDALRaster, filename, read_only = TRUE, open_options = NULL, shared = TRUE)

Exposes methods of Dataset and RasterBand, e.g. …

ds$info()
ds$getDriverShortName()
ds$getDriverLongName()
ds$getRasterXSize()
ds$getRasterYSize()
ds$getRasterCount()
ds$addBand(dataType, options)
ds$getGeoTransform()
ds$setGeoTransform(transform)
ds$getProjection()
ds$setProjection(projection)
ds$getDescription(band)
ds$setDescription(band, desc)
ds$getBlockSize(band)
ds$getDataTypeName(band)
ds$getNoDataValue(band)
ds$setNoDataValue(band, value)
ds$deleteNoDataValue(band)
ds$getMetadata(band, dom)
ds$setMetadata(band, metadata, dom)
ds$getMaskFlags(band)
# many more...

GDALRaster pixel I/O methods

  • Wraps GDALRasterIO()

  • Helpers to iterate on block boundaries, or user-defined multi-block chunks

ds$get_block_indexing(band)
ds$make_chunk_index(band, max_pixels)

ds$read(band, xoff, yoff, xsize, ysize, out_xsize, out_ysize)
ds$readBlock(band, xblockoff, yblockoff)
ds$readChunk(band, chunk_def)

ds$write(band, xoff, yoff, xsize, ysize, rasterData)
ds$writeBlock(band, xblockoff, yblockoff, rasterData)
ds$writeChunk(band, chunk_def, rasterData)

Raster creation techniques

  • If the format driver supports: gdal_formats()
  • Return a live dataset object
    → enables GDAL in-memory (MEM) format in R
  • Programmatic: create(), createCopy()
  • Utilities: translate(), warp(), …
  • GDAL CLI algorithms
    • API usage of the CLI framework operates on dataset objects, not just files

No-copy MEM raster from R data

rvector_to_MEM()

  • Expose pixel data in an existing R vector to GDAL
  • Pointer to R’s underlying C array (writable)
  • Safe object lifetimes
v <- sample(0:255, 20, replace = TRUE)
(ds_mem <- rvector_to_MEM(v, xsize = 5, ysize = 4))
#> C++ object of class <GDALRaster>
#>   • Driver: In Memory Raster (MEM)
#>   • DSN: <data pointer>
#>   • Dimensions: 5, 4, 1
#>   • CRS: not set
#>   • Pixel resolution: 1.000000, 1.000000
#>   • Bbox: 0.000000, 0.000000, 5.000000, 4.000000

Raster tools

  • autoCreateWarpedVRT()
  • buildVRT()
  • calc()
  • combine()
  • dem_proc()
  • fillNodata()
  • footprint()
  • is_los_visible()
  • make_chunk_index()
  • polygonize()
  • rasterize()
  • sieveFilter()
  • translate()
  • warp()
  • buildRAT()
  • mdim_as_classic()
  • pixel_extract()
  • bandCopyWholeRaster()
  • copyDatasetFiles()
  • deleteDataset()
  • inspectDataset()
  • renameDataset()

Raster display

read_ds(rast, bands = c(7, 5, 4)) |>
  plot_raster(
    minmax_pct_cut = c(2, 98),
    main = "Mount Hood July 2022 Landsat 7-5-4 (SWIR)")

Raster Attribute Tables

Read/display/build/write

Vector Data Model

Feature-complete Vector API

  • Dataset creation
  • Test dataset/layer capabilities
  • Complete schema management:
    • Layers and fields
    • Create/delete/rename
    • Read/write field domains
    • Execute SQL
  • GDAL utilities: ogrinfo(), ogr2ogr()
  • GDAL algorithms
    • Vector geoprocessing
  • Vector display
  • Persistent dataset connection
  • Efficient I/O
  • Attribute and spatial filters
  • Write:
    • New features in a layer
    • Upsert
    • Edit/overwrite
    • Delete
  • Transactions
  • Arrow Array Stream read/write
  • Create/append/inspect SOZip

Class GDALVector

Pointers to one OGRLayer and its Dataset

Constructor w/ optional args:

lyr <- new(GDALVector, dsn, layer, read_only, open_options, spatial_filter, dialect)

Exposes methods on OGRLayer, e.g. …

lyr$getName()
lyr$getFieldNames()
lyr$testCapability()
lyr$getFIDColumn()
lyr$getGeomType()
lyr$getGeometryColumn()
lyr$getSpatialRef()
lyr$bbox()
lyr$getLayerDefn()
lyr$getFieldDomain(domain_name)
lyr$setAttributeFilter(query)
lyr$getAttributeFilter()
lyr$setIgnoredFields(fields)
lyr$setSelectedFields(fields)
lyr$getIgnoredFields()
lyr$setSpatialFilter(wkt)
lyr$setSpatialFilterRect(bbox)
lyr$getSpatialFilter()
lyr$getFeatureCount()
# ...

GDALVector I/O methods

  • Traditional row iteration
  • GDAL’s Arrow-based column-oriented API
  • R-side ArrowArrayStream delivered by nanoarrow
lyr$getNextFeature()
lyr$setNextByIndex(i)
lyr$getFeature(fid)
lyr$resetReading()

lyr$fetch(n)  # analog of DBI::dbFetch()

lyr$setFeature(feature)
lyr$createFeature(feature)
lyr$batchCreateFeature(feature_set)
lyr$upsertFeature(feature)
lyr$deleteFeature(fid)
lyr$startTransaction()
lyr$commitTransaction()
lyr$rollbackTransaction()

lyr$getArrowStream()
lyr$releaseArrowStream()
lyr$writeArrowBatch(df)

Vector read benchmarks gdalraster logo

  • NZ Building Outlines, 1.5 GB GeoPackage
  • Load 3.3 million polygon features with 13 attributes each
Method Time (s)
gdalraster row-oriented → base R data frame 11.2
gdalraster row-oriented + convert to sf 25.0
sf row-oriented → sf data frame 77.2
gdalraster Arrow stream → base R data frame 2.9
sf Arrow stream → sf data frame 11.2

Spatial Reference System API

Spatial reference and transformation

PROJ via GDAL Spatial Reference System API

 
  • Set of srs_*()functions
    • Convert among SRS formats
    • Query attributes of an SRS definition
    • Query the PROJ database
  • Transform coordinates, geometries, bounding boxes

Geometry API

Geometry API operating on WKB/WKT

GEOS via GDAL Geometry API

  GEOS
  • Flexibly operate on either WKB or WKT
    • Lists of WKB raw vectors or WKT strings
    • g_wk2wk(): auto convert one to the other
  • g_*() functions: geom factory, inspect/test, binary predicates, binary ops, unary ops, measurements, utilities

Geom API is convenient for interactive use and suitable for processing

Plotting uses the wk package

# https://postgis.net/docs/ST_PointOnSurface.html
g <- "POLYGON ((130 120, 120 190, 30 140, 50 20, 190 20,
                170 100, 90 60, 90 130, 130 120))"
pt <- g_point_on_surface(g)
plot_geom(g, main = "https://postgis.net/docs/ST_PointOnSurface.html")
plot_geom(pt, pch = 16, add = TRUE)
centroid <- g_create("POINT", g_centroid(g))
plot_geom(centroid, add = TRUE)

GDAL Virtual File Systems

Bindings to GDAL’s Virtual Systems Interface (VSI)

File system operations abstracted

  • URLs, cloud storage services, Zip/GZip/7z, in-memory files, … as well as regular local file systems
    → Single interface works the same for any storage backend
  • 23 low-level vsi_*() functions
  • Class VSIFile wrapping VSIVirtualHandle

Virtual file systems vsi_*()

Includes analogs of POSIX file system functions and other utilities, e.g. …

vsi_copy_file()
vsi_glob()
vsi_mkdir()
vsi_read_dir()
vsi_rename()
vsi_rmdir()
vsi_stat()
vsi_sync()
vsi_unlink()
vsi_curl_clear_cache()
vsi_get_disk_free_space()
vsi_get_file_metadata()
vsi_get_fs_options()
vsi_is_local()
vsi_set_path_option()
vsi_supports_rnd_write()
vsi_supports_seq_write()
vsi_uri_to_vsi_path()

Class VSIFile: Standard C file I/O

Wraps VSIVirtualHandle: Binary file I/O through VSI

Constructor w/ optional args:

vf <- new(VSIFile, filename, access, options)

Methods of VSIFile

vf$seek(offset, origin)
vf$tell()
vf$rewind()
vf$read(nbytes)
vf$write(object)
vf$eof()
vf$truncate(new_size)
vf$flush()
vf$ingest(max_size)
vf$close()
vf$open()
vf$get_filename()
vf$get_access()
vf$set_access(access)

GDAL CLI Bindings

Unified CLI (GDAL > 3.11)

GDAL’s new framework for algorithms and utilities 🚀

  • Commands and sub-commands with consistent interface
  • GDALAlgorithm - API usage
    • Automatic discovery
    • Invoke programmatically
    • Operate on dataset objects
  • New GDAL capabilities → CLI

Class GDALAlg low-level API

Wraps the GDALAlgorithm (CLI) C API

Constructor w/ optional args:

alg <- new(GDALAlg, cmd, args)

Methods of GDALAlg

alg$info()
alg$argInfo(arg_name)
alg$usage()
alg$usageAsJSON()

alg$setArg(arg_name, arg_value)
alg$parseCommandLineArgs()
alg$getExplicitlySetArgs()
alg$run()
alg$output()
alg$outputs()

alg$close()
alg$release()

Function interface to GDAL CLI

  • Algorithm discovery and usage info
  • Run an algorithm or CLI pipeline and obtain its output
gdal_commands(contains = NULL, recurse = TRUE, cout = TRUE)

gdal_usage(cmd = NULL)

gdal_run(
  cmd,
  args,
  close = FALSE,
  quiet = FALSE,
  setVectorArgsFromObject = TRUE
)

Algorithm discovery

gdal_commands() returns a data frame invisibly

vec <- gdal_commands("vector")
#> → "vector"
#> ℹ Vector commands.
#> ℹ <https://gdal.org/programs/gdal_vector.html>

#> → "vector buffer"
#> ℹ Compute a buffer around geometries of a vector dataset.
#> ℹ <https://gdal.org/programs/gdal_vector_buffer.html>

#> → "vector check-coverage"
#> ℹ Check a polygon coverage for validity
#> ℹ <https://gdal.org/programs/gdal_vector_check_coverage.html>

#> → "vector check-geometry"
#> ℹ Check a dataset for invalid geometries
#> ℹ <https://gdal.org/programs/gdal_vector_check_geometry.html>

#> → "vector clean-coverage"
#> ℹ Alter polygon boundaries to make shared edges identical, removing gaps and overlaps
#> ℹ <https://gdal.org/programs/gdal_vector_clean_coverage.html>

#> → "vector clip"
#> ℹ Clip a vector dataset.
#> ℹ <https://gdal.org/programs/gdal_vector_clip.html>

#> → "vector combine"
#> ℹ Combine features into collections
#> ℹ <https://gdal.org/programs/gdal_vector_combine.html>

#> → "vector compare"
#> ℹ Compare two vector datasets.
#> ℹ <https://gdal.org/programs/gdal_vector_compare.html>

#> → "vector concat"
#> ℹ Concatenate vector datasets.
#> ℹ <https://gdal.org/programs/gdal_vector_concat.html>

#> → "vector concave-hull"
#> ℹ Compute the concave hull of geometries of a vector dataset.
#> ℹ <https://gdal.org/programs/gdal_vector_concave_hull.html>

#> → "vector convert"
#> ℹ Convert a vector dataset.
#> ℹ <https://gdal.org/programs/gdal_vector_convert.html>

#> → "vector convex-hull"
#> ℹ Compute the convex hull of geometries of a vector dataset.
#> ℹ <https://gdal.org/programs/gdal_vector_convex_hull.html>

#> → "vector create"
#> ℹ Create a vector dataset.
#> ℹ <https://gdal.org/programs/gdal_vector_create.html>

#> → "vector dissolve"
#> ℹ Dissolves multipart features
#> ℹ <https://gdal.org/programs/gdal_vector_dissolve.html>

#> → "vector edit"
#> ℹ Edit metadata of a vector dataset.
#> ℹ <https://gdal.org/programs/gdal_vector_edit.html>

#> → "vector explode"
#> ℹ Explode fields or geometries of a vector dataset.
#> ℹ <https://gdal.org/programs/gdal_vector_explode.html>

#> → "vector explode-collections"
#> ℹ Explode geometries of type collection of a vector dataset.
#> ℹ <https://gdal.org/programs/gdal_vector_explode_collections.html>

#> → "vector export-schema"
#> ℹ Export the OGR_SCHEMA from a vector dataset.
#> ℹ <https://gdal.org/programs/gdal_vector_export_schema.html>

#> → "vector filter"
#> ℹ Filter a vector dataset.
#> ℹ <https://gdal.org/programs/gdal_vector_filter.html>

#> → "vector grid"
#> ℹ Create a regular grid from scattered points.
#> ℹ <https://gdal.org/programs/gdal_vector_grid.html>

#> → "vector grid average"
#> ℹ Create a regular grid from scattered points using moving average interpolation.
#> ℹ <https://gdal.org/programs/gdal_vector_grid.html>

#> → "vector grid average-distance"
#> ℹ Create a regular grid from scattered points using the average distance between the grid node (center of the search ellipse) and all of the data points in the search ellipse.
#> ℹ <https://gdal.org/programs/gdal_vector_grid.html>

#> → "vector grid average-distance-points"
#> ℹ Create a regular grid from scattered points using the average distance between the data points in the search ellipse.
#> ℹ <https://gdal.org/programs/gdal_vector_grid.html>

#> → "vector grid count"
#> ℹ Create a regular grid from scattered points using the number of points in the search ellipse.
#> ℹ <https://gdal.org/programs/gdal_vector_grid.html>

#> → "vector grid invdist"
#> ℹ Create a regular grid from scattered points using weighted inverse distance interpolation.
#> ℹ <https://gdal.org/programs/gdal_vector_grid.html>

#> → "vector grid invdistnn"
#> ℹ Create a regular grid from scattered points using weighted inverse distance interpolation nearest neighbour.
#> ℹ <https://gdal.org/programs/gdal_vector_grid.html>

#> → "vector grid linear"
#> ℹ Create a regular grid from scattered points using linear/barycentric interpolation.
#> ℹ <https://gdal.org/programs/gdal_vector_grid.html>

#> → "vector grid maximum"
#> ℹ Create a regular grid from scattered points using the maximum value in the search ellipse.
#> ℹ <https://gdal.org/programs/gdal_vector_grid.html>

#> → "vector grid minimum"
#> ℹ Create a regular grid from scattered points using the minimum value in the search ellipse.
#> ℹ <https://gdal.org/programs/gdal_vector_grid.html>

#> → "vector grid nearest"
#> ℹ Create a regular grid from scattered points using nearest neighbor interpolation.
#> ℹ <https://gdal.org/programs/gdal_vector_grid.html>

#> → "vector grid range"
#> ℹ Create a regular grid from scattered points using the difference between the minimum and maximum values in the search ellipse.
#> ℹ <https://gdal.org/programs/gdal_vector_grid.html>

#> → "vector index"
#> ℹ Create a vector index of vector datasets.
#> ℹ <https://gdal.org/programs/gdal_vector_index.html>

#> → "vector info"
#> ℹ Return information on a vector dataset.
#> ℹ <https://gdal.org/programs/gdal_vector_info.html>

#> → "vector layer-algebra"
#> ℹ Perform algebraic operation between 2 layers.
#> ℹ <https://gdal.org/programs/gdal_vector_layer_algebra.html>

#> → "vector make-point"
#> ℹ Create point geometries from attribute fields
#> ℹ <https://gdal.org/programs/gdal_vector_make_point.html>

#> → "vector make-valid"
#> ℹ Fix validity of geometries of a vector dataset.
#> ℹ <https://gdal.org/programs/gdal_vector_make_valid.html>

#> → "vector partition"
#> ℹ Partition a vector dataset into multiple files.
#> ℹ <https://gdal.org/programs/gdal_vector_partition.html>

#> → "vector pipeline"
#> ℹ Process a vector dataset applying several steps.
#> ℹ <https://gdal.org/programs/gdal_vector_pipeline.html>

#> → "vector rasterize"
#> ℹ Burns vector geometries into a raster.
#> ℹ <https://gdal.org/programs/gdal_vector_rasterize.html>

#> → "vector rename-layer"
#> ℹ Rename layer(s) of a vector dataset.
#> ℹ <https://gdal.org/programs/gdal_vector_rename_layer.html>

#> → "vector reproject"
#> ℹ Reproject a vector dataset.
#> ℹ <https://gdal.org/programs/gdal_vector_reproject.html>

#> → "vector segmentize"
#> ℹ Segmentize geometries of a vector dataset.
#> ℹ <https://gdal.org/programs/gdal_vector_segmentize.html>

#> → "vector select"
#> ℹ Select a subset of fields from a vector dataset.
#> ℹ <https://gdal.org/programs/gdal_vector_select.html>

#> → "vector set-field-type"
#> ℹ Modify the type of a field of a vector dataset.
#> ℹ <https://gdal.org/programs/gdal_vector_set_field_type.html>

#> → "vector set-geom-type"
#> ℹ Modify the geometry type of a vector dataset.
#> ℹ <https://gdal.org/programs/gdal_vector_set_geom_type.html>

#> → "vector simplify"
#> ℹ Simplify geometries of a vector dataset.
#> ℹ <https://gdal.org/programs/gdal_vector_simplify.html>

#> → "vector simplify-coverage"
#> ℹ Simplify shared boundaries of a polygonal vector dataset.
#> ℹ <https://gdal.org/programs/gdal_vector_simplify_coverage.html>

#> → "vector sort"
#> ℹ Spatially order the features in a layer
#> ℹ <https://gdal.org/programs/gdal_vector_sort.html>

#> → "vector sql"
#> ℹ Apply SQL statement(s) to a dataset.
#> ℹ <https://gdal.org/programs/gdal_vector_sql.html>

#> → "vector swap-xy"
#> ℹ Swap X and Y coordinates of geometries of a vector dataset.
#> ℹ <https://gdal.org/programs/gdal_vector_swap_xy.html>

#> → "vector update"
#> ℹ Update an existing vector dataset with an input vector dataset.
#> ℹ <https://gdal.org/programs/gdal_vector_update.html>

head(vec)
#>                 command
#> 1                vector
#> 2         vector buffer
#> 3 vector check-coverage
#> 4 vector check-geometry
#> 5 vector clean-coverage
#> 6           vector clip
#>                                                                           description
#> 1                                                                    Vector commands.
#> 2                             Compute a buffer around geometries of a vector dataset.
#> 3                                               Check a polygon coverage for validity
#> 4                                              Check a dataset for invalid geometries
#> 5 Alter polygon boundaries to make shared edges identical, removing gaps and overlaps
#> 6                                                              Clip a vector dataset.
#>                                                         URL
#> 1                https://gdal.org/programs/gdal_vector.html
#> 2         https://gdal.org/programs/gdal_vector_buffer.html
#> 3 https://gdal.org/programs/gdal_vector_check_coverage.html
#> 4 https://gdal.org/programs/gdal_vector_check_geometry.html
#> 5 https://gdal.org/programs/gdal_vector_clean_coverage.html
#> 6           https://gdal.org/programs/gdal_vector_clip.html

Usage info

gdal_usage() prints to console (built on the low-level API)

gdal_usage("vector rasterize")
#> 
#> Usage: vector rasterize [OPTIONS] <INPUT> <OUTPUT>
#> 
#> Burns vector geometries into a raster. 
#> 
#> Positional arguments:
#>   -i, --input <INPUT>
#>     Input vector datasets
#>     [required]
#>   -o, --output <OUTPUT>
#>     Output raster dataset
#>     [required]
#> 
#> Common options:
#>   -q, --quiet
#>     Quiet mode (no progress bar or warning message)
#> 
#> Options:
#>   -f, --of, --format, --output-format <OUTPUT-FORMAT>
#>     Output format
#>   --co, --creation-option <KEY>=<VALUE>
#>     Creation option
#>     [0 or more values]
#>     [packed values not allowed, repeated arg allowed]
#>   --overwrite
#>     Whether overwriting existing output dataset is allowed
#>     [default: FALSE]
#>   -b, --band <BAND>
#>     The band(s) to burn values into (1-based index)
#>     [0 or more values]
#>     [packed values allowed, repeated arg allowed]
#>   --invert
#>     Invert the rasterization
#>     [default: FALSE]
#>   --all-touched
#>     Enables the ALL_TOUCHED rasterization option
#>   --burn <BURN>
#>     Burn value
#>     [0 or more values]
#>     [packed values allowed, repeated arg allowed]
#>   -a, --attribute-name <ATTRIBUTE-NAME>
#>     Attribute name
#>   --3d
#>     Indicates that a burn value should be extracted from the Z values of the feature
#>   -l, --input-layer <INPUT-LAYER>
#>     Input layer name
#>     [mutually exclusive with --sql]
#>   --where <WHERE>
#>     SQL where clause
#>   --sql <SQL>
#>     SQL select statement
#>     [mutually exclusive with --input-layer]
#>   --dialect <DIALECT>
#>     SQL dialect
#>   --nodata <NODATA>
#>     Assign a specified nodata value to output bands
#>   --init <INIT>
#>     Pre-initialize output bands with specified value
#>     [0 or more values]
#>     [packed values allowed, repeated arg allowed]
#>   --like <DATASET>
#>     Dataset to use as a template for bounds, CRS and resolution
#>   --crs <CRS>
#>     Override the projection for the output file
#>   --transformer-option <NAME>=<VALUE>
#>     Set a transformer option suitable to pass to GDALCreateGenImgProjTransformer2
#>     [0 or more values]
#>     [packed values allowed, repeated arg allowed]
#>   --extent <xmin>,<ymin>,<xmax>,<ymax>
#>     Set the target georeferenced extent
#>     [4 values]
#>     [packed values allowed, repeated arg not allowed]
#>   --resolution <xres>,<yres>
#>     Set the target resolution
#>     [2 values]
#>     [packed values allowed, repeated arg not allowed]
#>     [mutually exclusive with --size]
#>   --tap, --target-aligned-pixels
#>     (target aligned pixels) Align the coordinates of the extent of the output file to the values of the resolution
#>   --size <xsize>,<ysize>
#>     Set the target size in pixels and lines
#>     [2 values]
#>     [packed values allowed, repeated arg not allowed]
#>     [mutually exclusive with --resolution]
#>   --ot, --datatype, --output-data-type <OUTPUT-DATA-TYPE>
#>     Output data type
#>     [UInt8|Int8|UInt16|Int16|UInt32|Int32|UInt64|Int64|CInt16|CInt32|Float16|Float32|Float64|CFloat32|CFloat64]
#>   --optimization <OPTIMIZATION>
#>     Force the algorithm used (results are identical)
#>     [AUTO|RASTER|VECTOR]
#>     [default: AUTO]
#>   --add
#>     Add to existing raster
#>     [default: FALSE]
#>   --update
#>     Whether to open existing dataset in update mode
#>     [default: FALSE]
#> 
#> Advanced options:
#>   --oo, --open-option <KEY>=<VALUE>
#>     Open options
#>     [0 or more values]
#>     [packed values not allowed, repeated arg allowed]
#>   --if, --input-format <INPUT-FORMAT>
#>     Input formats
#>     [0 or more values]
#>     [packed values allowed, repeated arg allowed]
#> 
#> For more details: <https://gdal.org/programs/gdal_vector_rasterize.html>

gdal_run() → object of class GDALAlg

Can assign the return value if algorithm object is needed:

## convert storml_elev.tif to GeoPackage raster
f_tif <- system.file("extdata/storml_elev.tif", package = "gdalraster")
f_gpkg <- file.path(tempdir(), "storml_elev.gpkg")

args <- c("--overwrite", f_tif, f_gpkg)
alg <- gdal_run("raster convert", args)
#> 0...10...20...30...40...50...60...70...80...90...100 - done.

# obtain the output as a GDALRaster object
(ds <- alg$output())
#> C++ object of class GDALRaster
#>  Driver : GeoPackage (GPKG)
#>  DSN    : /tmp/Rtmp1MQDgw/storml_elev.gpkg
#>  Dim    : 143, 107, 1
#>  CRS    : NAD83 / UTM zone 12N (EPSG:26912)
#>  Res    : 30.000000, 30.000000
#>  Bbox   : 323476.071971, 5101871.983031, 327766.071971, 5105081.983031

gdal_run_piped() → output object

Piped variation - pass the output dataset to plot_raster():

# Visualize terrain ruggedness index (TRI)
system.file("extdata/storml_elev.tif", package="gdalraster") |>
  gdal_run_piped("raster tri", "", "MEM") |>
  plot_raster(legend = TRUE,
              minmax_pct_cut = c(0, 99),
              col_map_fn = rev(ltc::ltc("heatmap3")),
              main = "Storm Lake AOI: terrain ruggedness index")

GDAL CLI in an R native pipeline |>

Using in-memory GDAL datasets (MEM):

f_elev |>
  gdal_run_piped("raster color-map", "", "MEM",
    list(color_map = f_pal)
  ) |>
  gdal_run_piped("raster blend", "", "MEM",
    list(operator = "hsv-value",
         overlay =
           gdal_run_piped(f_elev, "raster hillshade", "", "MEM",
             list(zfactor = 1.5)))
  ) |>
  plot_raster(main = "Storm Lake AOI colorized shaded relief")

Wrap Up

Stable and well tested

On CRAN since May 2023

Focused project

Tracking upstream development

Has stickers 😊

💙 Made by hand w/ love & care 💙

https://github.com/firelab/gdalraster

Huge thanks       GEOS