GDAL7

The goal of GDAL7 is to model the GDAL API in R, generated from GDAL's own SWIG interface definitions and bound with cpp11 and S7.

Installation

You need GDAL with its development headers (gdal-config on its PATH, or pkg-config able to find gdal), then:

remotes::install_github("rgdal-dev/GDAL7")

Or from a clone:

git clone https://github.com/rgdal-dev/GDAL7
R CMD INSTALL GDAL7

The cpp11 registration files (src/cpp11.cpp, R/cpp11.R) are committed, so no generation step is needed to install. You only need cpp11::cpp_register() after adding or changing a [[cpp11::register]] function; data-raw/orchestrate.R does that as part of regenerating the bindings.

Example

The examples below run against a small GeoTIFF that ships with the package, so they need no network.

What an object knows about itself is a property, read from GDAL at the moment it is asked for rather than copied when the object was made, and settable by assignment where GDAL lets it be set. What an object does, and anything that needs an argument, is a function.

library(GDAL7)

dsn <- system.file("extdata/test.tif", package = "GDAL7")
ds <- gdal_open(dsn)

#### What every dataset, band and driver carries

basename(ds@description)                # what the object calls itself
#> [1] "test.tif"
ds@metadata_domain_list                 # the metadata domains it has
#> [1] "IMAGE_STRUCTURE"     "DERIVED_SUBDATASETS" ""

get_metadata_list(ds, "")               # metadata as KEY=VALUE strings
#> [1] "AREA_OR_POINT=Area"
get_metadata_dict(ds, "")               # the same, as a named vector
#> AREA_OR_POINT 
#>        "Area"
get_metadata_item(ds, "AREA_OR_POINT")  # one item, NA when not set
#> [1] "Area"

#### What a dataset knows

substr(ds@crs, 1, 40)                   # the CRS as WKT2
#> [1] "GEOGCRS[\"WGS 84\",ENSEMBLE[\"World Geodeti"
substr(ds@projection, 1, 40)            # GDAL's own WKT1 spelling
#> [1] "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROI"
ds@geotransform                         # pixel and line to x and y
#>        origin_x     pixel_width    row_rotation        origin_y column_rotation 
#>            -180              18               0              90               0 
#>    pixel_height 
#>             -18
basename(ds@file_list)                  # the files it is made of
#> [1] "test.tif"
ds@gcp_count
#> [1] 0
ds@layer_count                          # vector layers, none here
#> [1] 0

c(ds@raster_xsize, ds@raster_ysize, ds@raster_count)
#> [1] 20 10  2

#### Bands

# A band is taken by number, so that is a function; what it then knows about
# itself is a property.
band <- get_raster_band(ds, 1L)
band
#> <GDALRasterBand>
#>   Band:       1
#>   Size:       20 x 10
#>   Type:       Int16
#>   Block size: 20 x 10
#>   NoData:     -32768
#>   Color:      Gray

band@data_type_name
#> [1] "Int16"
band@block_size
#>  x  y 
#> 20 10
band@nodata_value
#> [1] -32768

A band belongs to its dataset, so closing the dataset retires the band with it. Reaching for one afterwards is an error rather than a crash:

gdal_close(ds)
band@xsize
#> Error:
#> ! This GDALRasterBand cannot be used: the GDALDataset it belongs to has been closed

Reading pixels

A read names a window and, separately, the size to return it at. GDAL picks an overview level that can serve the output size, so a large window at a small output size costs only the bytes of that level. The fixture here is a small COG with two overview levels.

ds <- gdal_open(system.file("extdata/overviews.tif", package = "GDAL7"))
band <- get_raster_band(ds, 1)

band@overview_sizes
#>   xsize ysize
#> 1   256   128
#> 2   128    64

# 512x256 down to 4x2, averaged.
read_raster(band, out_size = c(4, 2), resample = "average")
#> [1]  318  446  574  702  830  958 1086 1214

# An 8x8 window down to 2x2. The window may be fractional.
read_raster(band, window = c(0, 0, 8, 8), out_size = c(2, 2),
            resample = "average")
#> [1]  8 12 24 28

gdal_close(ds)

gdal_info() gathers everything gdalinfo reports that does not need a pass over the pixels, in a single call rather than one per property:

info <- gdal_info(system.file("extdata/overviews.tif", package = "GDAL7"))
info$size
#> xsize ysize 
#>   512   256
info$geotransform
#> [1] -180.000000    0.703125    0.000000   90.000000    0.000000   -0.703125
info$band_info
#>   band  type block_x block_y nodata scale offset color unit overviews
#> 1    1 Int16     128     128 -32768    NA     NA  Gray              2

Remote data

Any DSN GDAL understands works, which is the point of the /vsicurl/ and service drivers: no download step, and only the bytes actually needed are read.

gebco <- "/vsicurl/https://data.source.coop/alexgleith/gebco-2024/GEBCO_2024.tif"

ds <- gdal_open(gebco)
c(ds@raster_xsize, ds@raster_ysize)
get_raster_band(ds, 1L)@block_size
gdal_close(ds)

wmts <- paste0(
  "WMTS:https://services.arcgisonline.com/arcgis/rest/services/",
  "World_Imagery/MapServer/WMTS/1.0.0/WMTSCapabilities.xml,layer=World_Imagery"
)

ds <- gdal_open(wmts)
ds@metadata_domain_list
gdal_close(ds)

Multidimensional data

NetCDF, Zarr and HDF5 hold arrays of any number of dimensions rather than bands of pixels. GDAL7 navigates them and reads them.

ds <- gdal_open(system.file("extdata/multidim.zarr", package = "GDAL7"),
                multidim = TRUE)
grp <- get_root_group(ds)
grp
#> <GDALGroup>
#>   Name: /
#>   Arrays (4): lon, temperature, lat, time

arr <- open_mdarray(grp, "temperature")
arr
#> <GDALMDArray>
#>   Name: temperature
#>   Type: Float64
#>   Dimensions: time=3, lat=4, lon=5
#>   Unit: degC
#>   NoData: -999

arr@dimensions
#>   name size         type direction indexed
#> 1 time    3     TEMPORAL              TRUE
#> 2  lat    4 HORIZONTAL_Y     NORTH    TRUE
#> 3  lon    5 HORIZONTAL_X      EAST    TRUE

A read takes an origin, a count along each dimension and a step. The dim of what comes back is the reverse of the array's own dimension order, and carries the dimension names, so there is never a question of which axis is which.

values <- read_mdarray(arr)
dim(values)
#>  lon  lat time 
#>    5    4    3

# The first time step, every second longitude.
read_mdarray(arr, start = c(1, 1, 1), count = c(1, 4, 3), step = c(1, 1, 2))
#> , , 1
#> 
#>      [,1] [,2] [,3] [,4]
#> [1,]    1    6   11   16
#> [2,]    3   NA   13   18
#> [3,]    5   10   15   20

Where the values sit, and what the format says about them:

arr@dimension_values
#> $time
#> [1] 0 1 2
#> 
#> $lat
#> [1] -40 -41 -42 -43
#> 
#> $lon
#> [1] 140 141 142 143 144

arr@attributes
#> $coordinates
#> [1] "lat lon"
#> 
#> $long_name
#> [1] "air temperature"
#> 
#> $valid_range
#> [1] -50  50

A view slices lazily, in GDAL's own syntax, and a two-dimensional array can be handed to the raster side of the package as an ordinary dataset.

first <- get_view(arr, "[0,:,:]")
first@dimensions
#>   name size         type direction indexed
#> 1  lat    4 HORIZONTAL_Y     NORTH    TRUE
#> 2  lon    5 HORIZONTAL_X      EAST    TRUE

raster <- as_classic_dataset(first)
c(raster@raster_xsize, raster@raster_ysize)
#> [1] 5 4

gdal_close(ds)

Lazy arrays

With the altarr package installed, as_altarr() turns an array, a path and an array name, or a band into an ordinary R array whose values stay in the source until something reads them. Base R indexing and reductions read it a batch of chunks at a time, and each batch is one GDAL read, announced first with advise_read() so that Zarr decodes the batch's chunks on GDAL's threads at once. The chunks are the source's own, which an array reports as block_size:

zarr <- system.file("extdata/multidim.zarr", package = "GDAL7")

ds <- gdal_open(zarr, multidim = TRUE)
open_mdarray(get_root_group(ds), "temperature")@block_size
#> time  lat  lon 
#>    1    4    5
gdal_close(ds)

x <- as_altarr(zarr, array = "temperature")
dim(x)
#> [1] 5 4 3

x[2, 3, ]               # one cell, every time step
#> [1] 12 32 52
sum(x, na.rm = TRUE)
#> [1] 1822
altarr::altarr_stats(x)[["fetch_calls"]]
#> [1] 3

Built from a path, the array saves with saveRDS() as a recipe of about a kilobyte, and reopens its source when it is first read in a new session. A band works the same way, and an overview is a band, so as_altarr(get_overview(band, 1)) is a lazy view of a reduced level.

Vector data

A whole layer arrives as a data frame in one call, through the column-oriented Arrow path GDAL added in 3.6, so there is no per-feature work in R at all. Geometry comes back as WKB, in a list column of raw vectors:

gpkg <- system.file("extdata/test.gpkg", package = "GDAL7")

ds <- gdal_open(gpkg)
ds@layers
#>     name geometry_type feature_count fast_arrow
#> 1 places         Point             5       TRUE

places <- read_vector(ds)
places[c("name", "population")]
#>        name population
#> 1    Hobart     247086
#> 2 Melbourne    5031195
#> 3    Sydney    5312163
#> 4     Perth    2141834
#> 5    Darwin     147255

gdal_close(ds)

Filters are set on the layer and apply to everything read afterwards. A spatial filter is what makes a large layer cheap, because a driver with a spatial index uses it rather than reading every feature:

ds <- gdal_open(gpkg)
layer <- get_layer(ds, "places")

set_filter(layer, where = "population > 1e6", bbox = c(140, -45, 155, -30))
read_vector(layer)$name
#> [1] "Melbourne" "Sydney"

read_vector(execute_sql(ds, "SELECT name FROM places ORDER BY name LIMIT 2"))
#>   OGC_FID   name
#> 1       1 Darwin
#> 2       2 Hobart
gdal_close(ds)

The same path runs the other way, so a data frame with a WKB column writes back out as a layer:

path <- tempfile(fileext = ".gpkg")
write_vector(places, path, layer = "places", crs = "EPSG:4326",
             geometry_type = "Point")

identical(read_vector(path), places)
#> [1] TRUE
unlink(path)

Creating and writing

A dataset is made with a size, a band count and a type, written into, and closed. Closing is what finishes the file.

path <- tempfile(fileext = ".tif")

ds <- gdal_create(path, 64, 48, bands = 1, type = "Float32")
ds@crs <- "EPSG:4326"
ds@geotransform <- c(-180, 360 / 64, 0, 90, 0, -180 / 48)

write_raster(ds, list(as.double(seq_len(64 * 48))))
gdal_close(ds)

info <- gdal_info(path)
info$size
#> xsize ysize 
#>    64    48
info$band_info[c("band", "type", "block_x", "block_y")]
#>   band    type block_x block_y
#> 1    1 Float32      64      32
substr(info$projection, 1, 30)
#> [1] "GEOGCS[\"WGS 84\",DATUM[\"WGS_198"

Creation options are a table, read out of the driver's own metadata, rather than a string of XML to parse or a list of names to remember:

options <- driver_options("GTiff")
options[options$name %in% c("COMPRESS", "TILED", "BLOCKXSIZE"),
        c("name", "type", "default")]
#>          name          type default
#> 1    COMPRESS string-select    <NA>
#> 20      TILED       boolean      NO
#> 24 BLOCKXSIZE           int     256

options$choices[[which(options$name == "COMPRESS")]][1:6]
#> [1] "NONE"      "LZW"       "PACKBITS"  "JPEG"      "CCITTRLE"  "CCITTFAX3"

They are asked about before anything is made, so a name the driver does not have, or a value it will not take, is an error rather than a file left behind:

validate_creation_options("GTiff", c(COMPRESS = "DEFLATE"))
#> [1] TRUE
validate_creation_options("GTiff", c(COMPRES = "DEFLATE"))
#> [1] FALSE

gdal_create(tempfile(fileext = ".tif"), 4, 4,
            options = c(COMPRESS = "NOT_A_CODEC"))
#> Error:
#> ! The GTiff driver does not take these creation options: 'NOT_A_CODEC' is an unexpected value for COMPRESS creation option of type string-select.

gdal_create_copy() puts an existing dataset through a driver, which is how a COG gets made: the driver builds the overviews and lays the file out itself.

cog_path <- tempfile(fileext = ".tif")

cog <- gdal_create_copy(path, cog_path, driver = "COG",
                        options = c(COMPRESS = "DEFLATE", BLOCKSIZE = "128"),
                        progress = FALSE)
gdal_close(cog)

reopened <- gdal_open(cog_path)
get_metadata_dict(reopened, "IMAGE_STRUCTURE")
#>      LAYOUT COMPRESSION  INTERLEAVE 
#>       "COG"   "DEFLATE"      "BAND"
get_raster_band(reopened, 1)@block_size
#>   x   y 
#> 128 128
gdal_close(reopened)

gdal_drivers() says which drivers can do which of those things, in one pass over the driver manager:

# Raster drivers that can make a dataset from nothing, rather than only copy
# an existing one.
drivers <- gdal_drivers(c("DCAP_RASTER", "DCAP_CREATE"))
nrow(drivers)
#> [1] 36
head(drivers[c("short_name", "copy", "vsi", "extensions")], 5)
#>   short_name  copy   vsi extensions
#> 1      GTiff  TRUE  TRUE   tif tiff
#> 2        VRT  TRUE  TRUE        vrt
#> 3       NITF  TRUE  TRUE        ntf
#> 4        HFA  TRUE  TRUE        img
#> 5        MEM FALSE FALSE       <NA>

Virtual file systems

Every GDAL path works, not only files on disk, and they compose. A dataset built at a /vsimem/ path never reaches the disk, and its bytes come back as raw:

ds <- gdal_create("/vsimem/small.tif", 4, 4, bands = 1, type = "Byte")
write_raster(ds, list(as.double(seq_len(16))))
gdal_close(ds)

vfs_stat("/vsimem/small.tif")$size
#> [1] 162
bytes <- vfs_read_file("/vsimem/small.tif")
rawToChar(bytes[1:2])
#> [1] "II"

vfs_unlink("/vsimem/small.tif")

GDAL reads its configuration options at the moment it needs them, so setting one changes how the next call behaves. with_gdal_config() sets them for one expression and puts back exactly what was there:

with_gdal_config(c(GDAL_CACHEMAX = "16"), gdal_config("GDAL_CACHEMAX"))
#> [1] "16"
gdal_config("GDAL_CACHEMAX")
#> NULL

Algorithms

GDAL 3.11 added a registry of the algorithms its own command line is built from. GDAL7 binds the registry rather than each utility, so what gdal can do arrives with GDAL rather than with a GDAL7 release.

gdal_has_algorithms() says whether this build has any. Having the API is not the same as having the algorithms: a GDAL can be built with them turned off, and before 3.12 they registered themselves as a side effect of being loaded, which a static link leaves out. The Windows build here links GDAL statically, and answers FALSE today; Linux and macOS have the algorithms.

gdal_algorithms()
#> [1] "convert"  "dataset"  "info"     "mdim"     "pipeline" "raster"   "vector"  
#> [8] "vsi"

head(gdal_algorithms("raster"), 12)
#>  [1] "as-features"  "aspect"       "blend"        "calc"         "clean-collar"
#>  [6] "clip"         "color-map"    "compare"      "contour"      "convert"     
#> [11] "create"       "edit"

Every algorithm describes itself, so there is no table of arguments in this package to fall out of date.

info <- gdal_algorithm_info("raster reproject")
info$description
#> [1] "Reproject a raster dataset."

head(info$arguments[c("name", "type", "required")], 8)
#>           name         type required
#> 1         help      boolean    FALSE
#> 2     help-doc      boolean    FALSE
#> 3   json-usage      boolean    FALSE
#> 4       config  string_list    FALSE
#> 5 input-format  string_list    FALSE
#> 6  open-option  string_list    FALSE
#> 7        input dataset_list     TRUE
#> 8        quiet      boolean    FALSE

Arguments go in by name. An algorithm told to write to memory hands back a dataset and touches no disk; one given a file name writes it, closes it and hands back the path.

ds <- gdal_open(system.file("extdata/test.tif", package = "GDAL7"))

reprojected <- gdal_run("raster reproject", list(
  input = ds,
  output = "",
  "output-format" = "MEM",
  "dst-crs" = "EPSG:3857",
  bbox = c(-180, -85, 180, 85),
  "bbox-crs" = "EPSG:4326"
), progress = FALSE)

c(reprojected@raster_xsize, reprojected@raster_ysize)
#> [1] 19 19
substr(reprojected@projection, 1, 40)
#> [1] "PROJCS[\"WGS 84 / Pseudo-Mercator\",GEOGCS"

gdal_close(reprojected)
gdal_close(ds)

Pipelines work the same way, with the steps written as GDAL writes them:

piped <- gdal_run("raster pipeline", list(
  pipeline = paste(
    "read", system.file("extdata/overviews.tif", package = "GDAL7"),
    "! reproject --dst-crs EPSG:3857 --resampling average",
    "! write --output-format MEM streamed"
  )
), progress = FALSE)
#> Warning: GDAL: Clamping output bounds to (-20037508.342789,-20037508.342789) ->
#> (20037508.342789, 20037508.342789)

c(piped@raster_xsize, piped@raster_ysize)
#> [1] 497 497
gdal_close(piped)

A long algorithm draws a progress bar and stops on Ctrl-C.

Constants and capabilities

The enumerators and metadata keys that GDAL declares come through as two named vectors, rather than two hundred exported names:

gdal_constants("GDT_")[1:6]
#> GDT_Unknown    GDT_Byte    GDT_Int8  GDT_UInt16   GDT_Int16  GDT_UInt32 
#>           0           1          14           2           3           4
gdal_string_constants("DCAP_")[1:3]
#>                      DCAP_OPEN                    DCAP_CREATE 
#>                    "DCAP_OPEN"                  "DCAP_CREATE" 
#>   DCAP_CREATE_MULTIDIMENSIONAL 
#> "DCAP_CREATE_MULTIDIMENSIONAL"

A few bindings call GDAL functions newer than the minimum GDAL7 requires. They always exist and always dispatch; calling one that the GDAL in use is too old for raises an error naming the release it needs. gdal7_capabilities() is how to ask first:

gdal_release()[["release"]]
#> [1] "3.12.4"
gdal7_capabilities()
#>                              binding   gdal available
#> 1     dataset_mark_suppress_on_close 3.12.0      TRUE
#> 2 dataset_get_close_reports_progress 3.13.0     FALSE
#> 3                 dataset_as_mdarray 3.12.0      TRUE
#> 4             dataset_is_thread_safe 3.10.0      TRUE
#> 5    dataset_get_thread_safe_dataset 3.10.0      TRUE

Code of Conduct

Please note that the GDAL7 project is released with a Contributor Code of Conduct. By contributing to this project, you agree to abide by its terms.