Bare minimum capability and library alignment for using GDAL R and Python together, using mdim or classic raster:
Single GDAL, single PROJ, single GEOS, R packages and Python packages linking the same libraries, reticulate bridging cleanly between them, modern numpy, modern Python, modern R, current GDAL master available or latest release.
https://github.com/hypertidy/gdal-r-ci/
Here’s a demonstration of what the motivation for those images is. Here we run docker with config to allow remote read of NetCDF, which is one of the motivators. When we start the container it gives a brief message.
docker run --rm -ti --security-opt seccomp=unconfined ghcr.io/hypertidy/gdal-r-python:dev
── ghcr.io/hypertidy/gdal-r-python ──
GDAL 3.13.0dev PROJ 9.8.1 GEOS 3.13.1
R 4.6.0
Py 3.12.3 numpy 2.4.4
Docs: https://github.com/hypertidy/gdal-r-ci#gdal-r-pythonVariants ‘:dev’ and ‘:latest’ (default) exist for ‘gdal-system’, ‘gdal-r’, ‘gdal-r-full’, and ‘gdal-r-python’.
dsn <- "/vsicurl/https://github.com/OSGeo/gdal/raw/refs/heads/master/autotest/gdrivers/data/netcdf/MODIS_ARRAY.nc"
library(gdalraster)## GDAL 3.13.0beta2 (released NA), GEOS 3.13.1, PROJ 9.8.1
info <- mdim_info(dsn)## {
## "type": "group",
## "driver": "netCDF",
## "name": "/",
## "dimensions": [
## {
## "name": "y",
## "full_name": "/y",
## "size": 4,
## "indexing_variable": {
## "y": {
## "full_name": "/y",
## "datatype": "Float64",
## "dimensions": [
## "/y"
## ],
## "dimension_size": [
## 4
## ],
## "nodata_value": "NaN"
## }
## }
## },
## {
## "name": "x",
## "full_name": "/x",
## "size": 4,
## "indexing_variable": {
## "x": {
## "full_name": "/x",
## "datatype": "Float64",
## "dimensions": [
## "/x"
## ],
## "dimension_size": [
## 4
## ],
## "nodata_value": "NaN"
## }
## }
## }
## ],
## "arrays": {
## "y": {
## "full_name": "/y",
## "datatype": "Float64",
## "dimensions": [
## "/y"
## ],
## "dimension_size": [
## 4
## ],
## "nodata_value": "NaN"
## },
## "x": {
## "full_name": "/x",
## "datatype": "Float64",
## "dimensions": [
## "/x"
## ],
## "dimension_size": [
## 4
## ],
## "nodata_value": "NaN"
## },
## "__xarray_dataarray_variable__": {
## "full_name": "/__xarray_dataarray_variable__",
## "datatype": "Int16",
## "dimensions": [
## "/y",
## "/x"
## ],
## "dimension_size": [
## 4,
## 4
## ],
## "attributes": {
## "crs": "+a=6371007.181 +b=6371007.181 +lon_0=0 +no_defs +proj=sinu +units=m +x_0=0 +y_0=0",
## "res": [231.656358263958, 231.65635826375001],
## "transform": [231.65635826395601, 0, -7274009.6494862903, 0, -231.65635826375001, 5050108.6101527503, 0, 0, 1]
## },
## "srs": {
## "wkt": "PROJCRS[\"unknown\",BASEGEOGCRS[\"unknown\",DATUM[\"unknown\",ELLIPSOID[\"unknown\",6371007.181,0,LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8901]]],CONVERSION[\"unknown\",METHOD[\"Sinusoidal\"],PARAMETER[\"Longitude of natural origin\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"False easting\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"(E)\",east,ORDER[1],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]],AXIS[\"(N)\",north,ORDER[2],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]]",
## "data_axis_to_srs_axis_mapping": [2, 1]
## }
## }
## },
## "structural_info": {
## "NC_FORMAT": "CLASSIC"
## }
## }
ds <- new(GDALRaster, dsn)
str(gdal_version())## chr [1:4] "GDAL 3.13.0beta2 \"Iowa City\", released 2018/99/99" "3130000" ...
str(proj_version())## List of 4
## $ name : chr "9.8.1"
## $ major: int 9
## $ minor: int 8
## $ patch: int 1
str(geos_version())## List of 4
## $ name : chr "3.13.1"
## $ major: int 3
## $ minor: int 13
## $ patch: int 1
reticulate::import("pyogrio")## Module(pyogrio)
reticulate::import("rasterio")$show_versions()## rasterio info:
## rasterio: 1.5.0
## GDAL: 3.13.0beta2
## PROJ: 9.8.1
## GEOS: 3.13.1
## PROJ DATA: /root/.local/share/proj:/usr/local/share/proj:/usr/local/share/proj
## GDAL DATA: None
##
## System:
## python: 3.12.3 (main, Mar 23 2026, 19:04:32) [GCC 13.3.0]
## executable: /opt/gdal-py/bin/python
## machine: Linux-6.8.0-110-generic-x86_64-with-glibc2.39
##
## Python deps:
## affine: 2.4.0
## attrs: 26.1.0
## certifi: 2026.4.22
## click: 8.3.3
## cligj: 0.7.2
## cython: 3.2.4
## numpy: 2.4.4
## setuptools: 82.0.1
reticulate::import("numpy")[["__version__"]]## [1] "2.4.4"
reticulate::import("osgeo.gdal_array") ## this would error if GDAL was built against numpy>2## Module(osgeo.gdal_array)
gdal <- reticulate::import("osgeo.gdal")
gdal$UseExceptions()
mdim <- gdal$OpenEx(dsn, gdal$OF_MULTIDIM_RASTER)
gdal$Open(dsn)## <osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at 0x7a630e5dca80> >
rg <- mdim$GetRootGroup()
rg$GetMDArrayFullNamesRecursive ()## [1] "/y" "/x"
## [3] "/__xarray_dataarray_variable__"
arr <- rg$OpenMDArrayFromFullname("/__xarray_dataarray_variable__")
str(arr$ReadAsArray())## int [1:4, 1:4] -32767 -32767 -32767 -32767 -32767 -32767 -32767 -32767 -32767 -32767 ...
To build this document I put the source .Rmd in ~/demo-doc/demo.Rmd
and do
docker run --rm -ti -v $HOME/demo-doc:/demo-doc --security-opt seccomp=unconfined \
ghcr.io/hypertidy/gdal-r-python:dev \
Rscript -e 'rmarkdown::render("/demo-doc/demo.Rmd", output_format = "github_document")'
rmarkdown is installed there along with all the other things we needed in this document.
We need the -v mount to persist the document to our host machine, and
the --security-opt because of the remote NetCDF required.
In GDAL itself we see the following capabilities:
library(gdalraster)
## many drivers (formats)
str(gdalraster::gdal_formats())## 'data.frame': 211 obs. of 17 variables:
## $ short_name : chr "DERIVED" "GTI" "SNAP_TIFF" "GTiff" ...
## $ extensions : chr "" "gti.gpkg gti.fgb gti" "" "tif tiff" ...
## $ raster : logi TRUE TRUE TRUE TRUE TRUE TRUE ...
## $ multidim_raster : logi FALSE FALSE FALSE FALSE FALSE FALSE ...
## $ vector : logi FALSE FALSE FALSE FALSE FALSE FALSE ...
## $ geography_network : logi FALSE FALSE FALSE FALSE FALSE FALSE ...
## $ rw_flag : chr "ro" "ro" "ro" "rw+u" ...
## $ virtual_io : logi FALSE TRUE TRUE TRUE TRUE TRUE ...
## $ subdatasets : logi FALSE FALSE FALSE TRUE FALSE FALSE ...
## $ long_name : chr "Derived datasets using VRT pixel functions" "GDAL Raster Tile Index" "Sentinel Application Processing GeoTIFF" "GeoTIFF" ...
## $ sql_dialects : chr "" "" "" "" ...
## $ creation_datatypes : chr "" "" "" "Byte Int8 UInt16 Int16 UInt32 Int32 Float32 Float64 CInt16 CInt32 CFloat32 CFloat64" ...
## $ creation_field_types : chr "" "" "" "" ...
## $ creation_field_subtypes: chr "" "" "" "" ...
## $ multiple_vec_layers : logi FALSE FALSE FALSE FALSE FALSE FALSE ...
## $ read_field_domains : logi FALSE FALSE FALSE FALSE FALSE FALSE ...
## $ creation_fld_dom_types : chr "" "" "" "" ...
## GEOS embedded
gdalraster::has_geos()## [1] TRUE
## Spatialite available
gdalraster::has_spatialite()## [1] TRUE
## CPLHTTP services
gdalraster::http_enabled()## [1] TRUE
## GeoParquet capability details
str(gdalraster::gdal_get_driver_md("Parquet"))## List of 27
## $ DCAP_VECTOR : chr "YES"
## $ DCAP_CREATE_LAYER : chr "YES"
## $ DMD_LONGNAME : chr "(Geo)Parquet"
## $ DMD_EXTENSIONS : chr "parquet"
## $ DMD_EXTENSION : chr "parquet"
## $ DMD_HELPTOPIC : chr "drivers/vector/parquet.html"
## $ DCAP_VIRTUALIO : chr "YES"
## $ DCAP_MEASURED_GEOMETRIES : chr "YES"
## $ DCAP_Z_GEOMETRIES : chr "YES"
## $ DCAP_REOPEN_AFTER_WRITE_REQUIRED: chr "YES"
## $ DCAP_CAN_READ_AFTER_DELETE : chr "YES"
## $ DCAP_CREATE_FIELD : chr "YES"
## $ DMD_CREATIONFIELDDATATYPES : chr "Integer Integer64 Real String Date Time DateTime Binary IntegerList Integer64List RealList StringList"
## $ DMD_CREATIONFIELDDATASUBTYPES : chr "Boolean Int16 Float32 JSON UUID"
## $ DMD_CREATION_FIELD_DEFN_FLAGS : chr "WidthPrecision Nullable Comment AlternativeName Domain"
## $ DMD_SUPPORTED_SQL_DIALECTS : chr "OGRSQL SQLITE"
## $ DMD_OPENOPTIONLIST : chr "<OpenOptionList> <Option name='GEOM_POSSIBLE_NAMES' type='string' description='Comma separated list of possibl"| __truncated__
## $ DCAP_OPEN : chr "YES"
## $ DCAP_CREATE : chr "YES"
## $ DCAP_UPDATE : chr "YES"
## $ DMD_UPDATE_ITEMS : chr "Features"
## $ DCAP_DELETE_FIELD : chr "YES"
## $ DCAP_REORDER_FIELDS : chr "YES"
## $ GDAL_DMD_ALTER_FIELD_DEFN_FLAGS : chr "Name Type WidthPrecision"
## $ ARROW_VERSION : chr "24.0.0"
## $ ARROW_DATASET : chr "YES"
## $ DS_LAYER_CREATIONOPTIONLIST : chr "<LayerCreationOptionList>\n <Option name=\"COMPRESSION\" type=\"string-select\" description=\"Compression meth"| __truncated__
## Zarr capability details
str(gdalraster::gdal_get_driver_md("Zarr"))## List of 22
## $ DCAP_RASTER : chr "YES"
## $ DCAP_MULTIDIM_RASTER : chr "YES"
## $ DMD_LONGNAME : chr "Zarr"
## $ DMD_EXTENSIONS : chr "zarr"
## $ DMD_EXTENSION : chr "zarr"
## $ DMD_CREATIONDATATYPES : chr "Int8 Byte Int16 UInt16 Int32 UInt32 Int64 UInt64 Float16 Float32 Float64 CFLoat16 CFloat32 CFloat64"
## $ DCAP_VIRTUALIO : chr "YES"
## $ DMD_SUBDATASETS : chr "YES"
## $ DCAP_CREATE_SUBDATASETS : chr "YES"
## $ DMD_OPENOPTIONLIST : chr "<OpenOptionList> <Option name='LIST_ALL_ARRAYS' type='boolean' description='Whether to list all arrays, and n"| __truncated__
## $ DMD_MULTIDIM_DATASET_CREATIONOPTIONLIST: chr "<MultiDimDatasetCreationOptionList> <Option name='FORMAT' type='string-select' default='ZARR_V2'> <Value>"| __truncated__
## $ DCAP_OPEN : chr "YES"
## $ DCAP_CREATE : chr "YES"
## $ DCAP_CREATECOPY : chr "YES"
## $ DCAP_CREATE_MULTIDIMENSIONAL : chr "YES"
## $ DCAP_UPDATE : chr "YES"
## $ DMD_UPDATE_ITEMS : chr "GeoTransform SRS NoData RasterValues DatasetMetadata BandMetadata"
## $ COMPRESSORS : chr "blosc,zlib,gzip,lzma,zstd,lz4"
## $ FILTERS : chr "delta"
## $ BLOSC_COMPRESSORS : chr "blosclz,lz4,lz4hc,snappy,zlib,zstd"
## $ DMD_MULTIDIM_ARRAY_CREATIONOPTIONLIST : chr "<MultiDimArrayCreationOptionList>\n <Option name=\"COMPRESS\" type=\"string-select\" description=\"Compression"| __truncated__
## $ DMD_CREATIONOPTIONLIST : chr "<CreationOptionList>\n <Option name=\"COMPRESS\" type=\"string-select\" description=\"Compression method\" def"| __truncated__
# if spatialite isn't compiled in, ST_AsText is unknown.
# We use a built-in shapefile shipped with gdalraster (no temp files).
dsn <- system.file("extdata/poly_multipoly.shp", package = "gdalraster")
result <- ogr_execute_sql(
dsn,
"SELECT ST_AsText(ST_PointFromText('POINT(1 2)')) AS pt",
dialect = "SQLite"
)
#the most direct possible test:
ver <- ogr_execute_sql(
dsn,
"SELECT spatialite_version() AS v",
dialect = "SQLite"
)
result$getNextFeature()## OGR feature (attributes)
## $FID
## integer64
## [1] 0
##
## $pt
## [1] POINT(1 2)
ver$getNextFeature()## OGR feature (attributes)
## $FID
## integer64
## [1] 0
##
## $v
## [1] 5.1.0
