Metadata-Version: 2.2
Name: gdalcubes
Version: 0.1.0.dev0
Summary: On-demand Earth observation data cubes from GDAL image collections.
License: MIT License
         
         Copyright (c) 2026 Felipe Carlos
         
         Permission is hereby granted, free of charge, to any person obtaining a copy
         of this software and associated documentation files (the "Software"), to deal
         in the Software without restriction, including without limitation the rights
         to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
         copies of the Software, and to permit persons to whom the Software is
         furnished to do so, subject to the following conditions:
         
         The above copyright notice and this permission notice shall be included in all
         copies or substantial portions of the Software.
         
         THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
         IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
         FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
         AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
         LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
         OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
         SOFTWARE.
         
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: MIT License
Classifier: Programming Language :: C++
Classifier: Programming Language :: Python :: 3
Classifier: Topic :: Scientific/Engineering :: GIS
Project-URL: Homepage, https://github.com/m3nin0-labs/gdalcubes.py
Requires-Python: >=3.10
Requires-Dist: numpy>=1.24
Requires-Dist: orjson>=3.9
Requires-Dist: tqdm>=4.64
Provides-Extra: dev
Requires-Dist: ruff; extra == "dev"
Provides-Extra: test
Requires-Dist: pytest; extra == "test"
Provides-Extra: stac
Requires-Dist: pystac-client; extra == "stac"
Requires-Dist: planetary-computer; extra == "stac"
Provides-Extra: xarray
Requires-Dist: xarray; extra == "xarray"
Requires-Dist: dask; extra == "xarray"
Requires-Dist: netCDF4; extra == "xarray"
Provides-Extra: tif
Requires-Dist: rioxarray; extra == "tif"
Description-Content-Type: text/markdown

## gdalcubes.py

Python bindings for [gdalcubes](https://github.com/appelmar/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:

```sh
brew install gdal netcdf
```

On Debian/Ubuntu:

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

Then install the package:

```sh
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:

```sh
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.

```python
import os
import gdalcubes as gc

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

### 2. Download the example imagery

```python
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()`.

```python
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.

```python
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.

```python
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.

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

### 7. Compute an NDVI composite

```python
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:

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

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

```python
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:

```python
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:

```python
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](https://corteva.github.io/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:

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

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

```python
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:

```python
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:

```python
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:

```python
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:

```python
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:

```sh
pip install gdalcubes[stac]
```

### 1. Search a catalog

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

```python
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

```python
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, ...).

```python
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:

```python
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

The project uses [uv](https://docs.astral.sh/uv/) for environment management. Each step below is a single command you can run from the repository root.

### 1. Build and install in editable mode

This compiles the C++ extension and installs the package with its test and optional dependencies:

```sh
uv pip install -e '.[test,xarray,tif]'
```

Re-run this command whenever you change C++ sources under `src/bindings/`.

### 2. Run the tests

```sh
uv run pytest
```

### 3. Format and lint the Python code

```sh
uvx ruff format gdalcubes examples
uvx ruff check gdalcubes
```

### 4. Type-check the Python code

```sh
uvx mypy gdalcubes
```

### 5. Format the C++ bindings

The binding sources under `src/bindings/` follow the `.clang-format` in that directory:

```sh
uvx clang-format -i src/bindings/*.cpp src/bindings/*.h
```

### 6. Re-vendor the C++ core

The core under `src/gdalcubes/` is a copy of the gdalcubes C++ sources, refreshed from a local checkout of the R package repository:

```sh
tools/sync_core.sh /path/to/gdalcubes
```

Local fixes that upstream does not have yet live in `tools/patches/` and are re-applied by the script after copying. If one no longer applies, upstream changed that file: see `tools/patches/README.md`.

Binary wheels (with GDAL and friends bundled, no system dependencies) are built by `.github/workflows/wheels.yml` via cibuildwheel.
