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:
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] -32768A 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 closedReading 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 2Remote 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 TRUEA 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 20Where 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 50A 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] 3Built 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")
#> NULLAlgorithms
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 FALSEArguments 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 TRUECode 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.