gdalcubes 0.1.0.dev0

gdalcubes.py

Python bindings for gdalcubes - On-demand Earth observation data cubes from GDAL image collections.

Installation

You need development files for GDAL, netCDF (C library), SQLite3, libcurl, with CMake ≥ 3.18 and a C++17 compiler.

On macOS:

brew install gdal netcdf

On Debian/Ubuntu:

apt install libgdal-dev libnetcdf-dev libsqlite3-dev libcurl4-openssl-dev

Then install the package:

pip install .

Pre-built wheels bundle GDAL and friends, so most users will not need any of the system packages above once those are published.

Optional extras

gdalcubes.py provides some integrations with the Python ecosystem. These integrations require extra Python dependencies. The table below lists these extra dependencies and when they are used:

Extra	Installs	For
gdalcubes[xarray]	xarray, dask, netCDF4	reading results with to_xarray() and extract_geom()
gdalcubes[tif]	rioxarray	reading GeoTIFF collections
gdalcubes[stac]	pystac-client, planetary-computer	building collections from STAC

After selecting the extra dependencies you want, you can install them with the following command:

pip install 'gdalcubes[xarray,tif,stac]'

Getting started

This is the same walkthrough as the R package README, using a public Landsat 8 surface reflectance dataset over the Amazon (~2 GB). The full script is in examples/local-files.py. Below we take it one step at a time.

1. Configure gdalcubes

Set how many threads to use and turn on a progress bar for the long-running steps.

import os
import gdalcubes as gc

gc.set_threads(os.cpu_count() or 1)
gc.set_progress(True)

2. Download the example imagery

import urllib.request
import zipfile

urllib.request.urlretrieve(
    "https://hs-bochum.sciebo.de/s/8XcKAmPfPGp2CYh/download", "L8_Amazon.zip"
)

with zipfile.ZipFile("L8_Amazon.zip") as archive:
    archive.extractall("L8_Amazon")

3. Build an image collection

Scanning the images once extracts their spatial extent, acquisition time, and band layout into an image collection, a small SQLite database that references the files (no second copy of the data).

How metadata is read from file names is described by a collection format. Presets for several Sentinel, Landsat, and MODIS products ship with the package. You can list them with gc.list_collection_formats().

import glob

files = glob.glob("L8_Amazon/**/*.tif", recursive=True)
col = gc.create_image_collection(files, fmt="L8_SR", out="L8.db")

4. Define a data cube view

A data cube view defines the geometry of the target cube independently of the imagery: spatial reference system, extent, pixel size, and temporal resolution.

Here we ask for a coarse overview, the full collection extent at 1 km pixels, one cell per year, and let gdalcubes reproject, crop, resample, and aggregate on the fly.

view = gc.cube_view(
    srs="EPSG:3857",
    extent=col,
    dx=1000,
    dy=1000,
    dt="P1Y",
    aggregation="median",
    resampling="bilinear",
)

5. Create the cube

raster_cube() returns immediately: it is a lazy proxy, nothing is computed yet.

cube = gc.raster_cube(col, view)

print(cube)
# <gdalcubes.Cube 10 band(s) [AEROSOL, B01, ..., RADSAT_QA],
#  t: 2013-01-01..2019-12-31 (n=7), y: ... (n=559), x: ... (n=783), srs: EPSG:3857>

6. Compute a median RGB composite

Operations compose as plain functions and stay lazy. First pick the three bands, then reduce the time axis to a single median per band.

rgb_median = gc.reduce_time(
    gc.select_bands(cube, ["B02", "B03", "B04"]),
    ["median(B02)", "median(B03)", "median(B04)"],
)

7. Compute an NDVI composite

ndvi = gc.apply_pixel(
    gc.select_bands(cube, ["B04", "B05"]),
    "(B05-B04)/(B05+B04)",
    names=["NDVI"],
)

8. Export the result

Evaluation starts here. Write a single netCDF file:

gc.write_netcdf(ndvi, "ndvi.nc")

Or one (optionally cloud-optimized) GeoTIFF per time slice, which any GDAL-based tool (rasterio, QGIS, ...) reads directly:

gc.write_tif(ndvi, "ndvi_tifs", prefix="NDVI_", cog=True)

9. Use the result with xarray

Step 8 already evaluated the cube and wrote it to disk, so reading it back costs nothing more. The function to_xarray() reads those files rather than evaluating again, and returns an xarray.Dataset backed by dask:

ds = gc.to_xarray(ndvi)

ds["NDVI"].isel(time=0).mean().compute()

It does not matter which format you wrote, netCDF or the GeoTIFF collection both work.

For a cube that has not been written yet, you can use the to_xarray() giving a path. It evaluates once and reads from there:

ds = gc.to_xarray(rgb_median, "rgb_median.nc")

The resulting xarray.Dataset has one variable per band, dimensions time, y, x, and a crs coordinate, so rioxarray picks up the projection directly.

Load from disk

To load a cube, netCDF or GeoTIFF, as an xarray.Dataset you can use the specialized open functions from gdalcubes.py. For example, to load an existing cube saved as netCDF you can use:

ds = gc.open_xarray("ndvi.nc")

To read a GeoTIFF collection from a directory, you can use:

ds = gc.open_tif_collection("ndvi_tifs", prefix="NDVI_")

Once you load it, you can continue the operations using the xarray / dask infrastructure.

10. Extract time series of feature locations

extract_geom() evaluates a cube and extracts values from feature locations specified in a vector data source (any format GDAL reads, or a geopandas.GeoDataFrame). The function returns a pandas.DataFrame with one row per feature, time step, and pixel. The fid column is the position of the feature in the input, so the rows join back to their features directly:

df = gc.extract_geom(ndvi, "samples.gpkg")

Pass time_column="date" to extract each feature at its own time step only, as recorded in that attribute.

11. Read chunks in memory

A cube is evaluated one chunk at a time, in the chunk size it was created with. The chunks_iter() hands those chunks to Python as the core computes them, so a cube can be processed without writing it to disk first:

for chunk in gc.chunks_iter(ndvi):
    if chunk.status is gc.ChunkStatus.ERROR:
        raise RuntimeError(f"chunk {chunk.info.id} failed to read")

    if chunk.values is None:
        continue

    print(chunk.info.start, chunk.values.shape)

The chunk.values is a C-contiguous float64 array of shape (bands, t, y, x), with NaN for no data and row 0 at the northern edge. It shares the core's own buffer, which it keeps alive, so nothing is copied.

Next, the chunk.info says where the chunk is: start and shape in cube coordinates, left, right, bottom and top in the cube's spatial reference system, and t0 and t1 for its first and last time step. The info.slices gives the same position as slices, which is all it takes to put the chunks back together:

import numpy as np

values = np.full(tuple(ndvi.size), np.nan)

for chunk in gc.chunks_iter(ndvi):
    if chunk.values is not None:
        values[(slice(None), *chunk.info.slices)] = chunk.values

To drive the loop yourself, or to read several chunks at once, use chunks_count() and chunks_read().

Reading releases the GIL while the core works, so Python threads read in parallel, exactly as the core does internally when it evaluates a whole cube:

from concurrent.futures import ThreadPoolExecutor

with ThreadPoolExecutor(max_workers=8) as pool:
    chunks = pool.map(lambda i: gc.chunks_read(ndvi, i), range(gc.chunks_count(ndvi)))

This is not faster than writing the cube, it is just a different way to access the cube data.

Data cubes from STAC

Instead of indexing local files, an image collection can be built directly from STAC search results. No image file is opened, the metadata comes from the STAC items, so indexing thousands of remote COGs is instant. The full example script is in examples/stac-mpc.py.

Install the extra client libraries:

pip install gdalcubes[stac]

1. Search a catalog

Search with any STAC client (here pystac-client). URL signing happens at search time.

import planetary_computer
import pystac_client

# create catalog client
catalog = pystac_client.Client.open(
    "https://planetarycomputer.microsoft.com/api/stac/v1",
    modifier=planetary_computer.sign_inplace,
)

# search
items = catalog.search(
    collections=["sentinel-2-l2a"],
    bbox=[6.1, 50.2, 6.9, 50.8],
    datetime="2023-06/2023-08",
)

items = items.item_collection()

2. Build the collection from the items

from gdalcubes import stac

col = gc.stac_image_collection(
    items, source = "MPC", assets=["B04", "B08", "SCL"]
)

3. Build a cloud-masked cube

The mask argument drops pixels whose SCL class is one of the listed values (cloud, shadow, ...).

view = gc.cube_view(
    srs="EPSG:32632",
    extent=col,
    dx=100,
    dy=100,
    dt="P1M",
    aggregation="median",
    resampling="bilinear",
)

cube = gc.raster_cube(
    col, view, mask={
        "band": "SCL",
        "values": [3, 8, 9, 10]
    }
)

Adapters

Because providers implement STAC in different ways, provider behavior is factored into an explicit adapter, a frozen bundle of pure functions (asset selection and band naming, href → GDAL URI, datetime extraction, item filtering, GDAL configuration) with spec-compliant defaults.

Presets are plain constants (stac.PLANETARY_COMPUTER, stac.EARTH_SEARCH).

You can customize the Default adapter by subclassing or changing attributes:

adapter = stac.DEFAULT.replace(resolve_uri=lambda href: "/vsicurl/" + href)

Note that signed URLs (e.g. Planetary Computer SAS tokens) expire, so STAC-built collections are best treated as ephemeral, the default is a temporary database, rebuilt from a fresh search when needed.

Development

Using uv, from the repository root:

uv pip install -e '.[test,xarray,tif]'   # build and install (re-run after C++ changes)
uv run pytest                            # run the tests
uvx ruff format gdalcubes examples       # format
uvx ruff check gdalcubes                 # lint
uvx mypy gdalcubes                       # type-check