On-Demand Earth System Data Cubes (ESDCs) in Python

---

GitHub: https://github.com/davemlz/cubo

Documentation: https://cubo.readthedocs.io/

PyPI: https://pypi.org/project/cubo/

Conda-forge: https://anaconda.org/conda-forge/cubo

Tutorials: https://cubo.readthedocs.io/en/latest/tutorials.html

Paper: https://doi.org/10.1109/IGARSS53475.2024.10640742

---

News

Important

⭐ Pinned (2026-02-12): Our cubo paper is out in IEEE! Check it here: Montero, D., Aybar, C., Ji, C., Kraemer, G., Sochting, M., Teber, K., & Mahecha, M.D. (2024). On-Demand Earth System Data Cubes.

Overview

SpatioTemporal Asset Catalogs (STAC) provide a standardized format that describes geospatial information. Multiple platforms are using this standard to provide clients several datasets. Nice platforms such as Planetary Computer use this standard. Additionally, Google Earth Engine (GEE) also provides a gigantic catalogue that users can harness for different tasks in Python.

cubo is a Python package that provides users of STAC and GEE an easy way to create On-Demand Earth System Data Cubes (ESDCs). This is perfectly suitable for Deep Learning (DL) tasks. You can easily create a lot of ESDCs by just knowing a pair of coordinates and the edge size of the cube in pixels!

Check the simple usage of cubo with STAC here:

import cubo
import xarray as xr

da = cubo.create(
    lat=4.31, # Central latitude of the cube
    lon=-76.2, # Central longitude of the cube
    collection="sentinel-2-l2a", # Name of the STAC collection
    bands=["B02","B03","B04"], # Bands to retrieve
    start_date="2021-06-01", # Start date of the cube
    end_date="2021-06-10", # End date of the cube (included)
    edge_size=64, # Edge size of the cube (px)
    resolution=10, # Pixel size of the cube (m)
)

This chunk of code just created an xr.DataArray object given a pair of coordinates, the edge size of the cube (in pixels), and additional information to get the data from STAC (Planetary Computer by default, but you can use another provider!). Note that you can also use the resolution you want (in meters) and the bands that you require.

Now check the simple usage of cubo with GEE here:

import cubo
import xarray as xr

da = cubo.create(
    lat=51.079225, # Central latitude of the cube
    lon=10.452173, # Central longitude of the cube
    collection="COPERNICUS/S2_SR_HARMONIZED", # Id of the GEE collection
    bands=["B2","B3","B4"], # Bands to retrieve
    start_date="2016-06-01", # Start date of the cube
    end_date="2017-07-01", # End date of the cube (included)
    edge_size=128, # Edge size of the cube (px)
    resolution=10, # Pixel size of the cube (m)
    gee=True # Use GEE instead of STAC
)

This chunk of code is very similar to the STAC-based cubo code. Note that the collection is now the ID of the GEE collection to use, and note that the gee argument must be set to True (or stac="ee").

And you can even combine multiple collections from STAC and GEE in one cube:

import cubo

dt = cubo.create(
    lat=51.3397, # Central latitude of the cube
    lon=12.3731, # Central longitude of the cube
    collection=["sentinel-2-l2a", "COPERNICUS/S2_CLOUD_PROBABILITY"], # Collections to use
    bands={"sentinel-2-l2a": ["B02", "B03", "B04"], "COPERNICUS/S2_CLOUD_PROBABILITY": "probability"}, # Bands of each collection
    stac={
        "https://planetarycomputer.microsoft.com/api/stac/v1": "sentinel-2-l2a", # STAC collections
        "ee": "COPERNICUS/S2_CLOUD_PROBABILITY", # GEE collections
    },
    start_date="2022-06-14", # Start date of the cube
    end_date="2022-06-17", # End date of the cube (included)
    edge_size=128, # Edge size of the cube (px)
    resolution=10, # Pixel size of the cube (m)
    datatree=True, # One node per collection
)

All collections share the same pixel grid, so they are co-registered by design!

How does it work?

The thing is super easy and simple.

1. You have the coordinates of a point of interest. The cube will be created around these coordinates (i.e., these coordinates will be approximately the spatial center of the cube).
2. Internally, the coordinates are transformed to the projected UTM coordinates [x,y] in meters (i.e., local UTM CRS). They are rounded to the closest pair of coordinates that are divisible by the resolution you requested.
3. The edge size you provide is used to create a Bounding Box (BBox) for the cube in the local UTM CRS given the exact amount of pixels (Note that the edge size should be a multiple of 2, otherwise it will be rounded, usual edge sizes for ML are 64, 128, 256, 512, etc.).
4. Additional information is used to retrieve the data from the STAC catalogue or from GEE: starts and end dates, name of the collection(s), endpoint of the catalogue(s), etc.
5. Then, by using stackstac and pystac_client the cube is retrieved as a xr.DataArray. In the case of GEE, the cube is retrieved via xee using exactly the same pixel grid. When using multiple collections, the cubes are merged, or returned as a xr.DataTree.
6. Success! That's what cubo is doing for you, and you just need to provide the coordinates, the edge size, and the additional info to get the cube.

Installation

Install the latest version from PyPI:

pip install cubo

Install cubo with the required GEE dependencies from PyPI:

pip install cubo[ee]

Upgrade cubo by running:

pip install -U cubo

Install the latest version from conda-forge:

conda install -c conda-forge cubo

Install the latest dev version from GitHub by running:

pip install git+https://github.com/davemlz/cubo

Features

Main function: create()

cubo is pretty straightforward, everything you need is in the create() function:

da = cubo.create(
    lat=4.31,
    lon=-76.2,
    collection="sentinel-2-l2a",
    bands=["B02","B03","B04"],
    start_date="2021-06-01",
    end_date="2021-06-10",
    edge_size=64,
    resolution=10,
)

Using different units for edge_size

By default, the units of edge_size are pixels. But you can modify this using the units argument:

da = cubo.create(
    lat=4.31,
    lon=-76.2,
    collection="sentinel-2-l2a",
    bands=["B02","B03","B04"],
    start_date="2021-06-01",
    end_date="2021-06-10",
    edge_size=1500,
    units="m",
    resolution=10,
)

Tip

You can use "px" (pixels), "m" (meters), or any unit available in scipy.constants.

da = cubo.create(
    lat=4.31,
    lon=-76.2,
    collection="sentinel-2-l2a",
    bands=["B02","B03","B04"],
    start_date="2021-06-01",
    end_date="2021-06-10",
    edge_size=1.5,
    units="kilo",
    resolution=10,
)

Using another endpoint

By default, cubo uses Planetary Computer. But you can use another STAC provider endpoint if you want:

da = cubo.create(
    lat=4.31,
    lon=-76.2,
    collection="sentinel-2-l2a",
    bands=["rededge1","rededge2","rededge3"],
    start_date="2020-01-01",
    end_date="2020-06-01",
    edge_size=128,
    resolution=20,
    stac="https://earth-search.aws.element84.com/v1"
)

Keywords for searching data

You can pass kwargs to pystac_client.Client.search() if required:

da = cubo.create(
    lat=4.31,
    lon=-76.2,
    collection="sentinel-2-l2a",
    bands=["B02","B03","B04"],
    start_date="2021-01-01",
    end_date="2021-06-10",
    edge_size=64,
    resolution=10,
    query={"eo:cloud_cover": {"lt": 40}} # kwarg to pass
)

Multiple collections

You can use multiple collections by passing a list to collection and a dictionary with the bands of each collection to bands. By default, the cubes are merged along the band dimension (the time dimension is the union of the timestamps of all collections):

da = cubo.create(
    lat=51.3397,
    lon=12.3731,
    collection=["sentinel-2-l2a", "sentinel-1-rtc"],
    bands={"sentinel-2-l2a": ["B04", "B08"], "sentinel-1-rtc": ["vv", "vh"]},
    start_date="2022-06-01",
    end_date="2022-06-15",
    edge_size=64,
    resolution=10,
)

Tip

If the collections come from different STAC Catalogues (or from GEE), use a dictionary in stac with endpoints as keys and the collection(s) of each endpoint as values. Use the "ee" key for GEE collections.

DataTrees and Datasets

Collections rarely share timestamps. Use datatree=True to get a xr.DataTree with one node per collection (each one with its own time and band dimensions, but sharing the spatial coordinates), and dataset=True to get one variable per band:

dt = cubo.create(
    lat=51.3397,
    lon=12.3731,
    collection=["sentinel-2-l2a", "sentinel-1-rtc"],
    bands={"sentinel-2-l2a": ["B04", "B08"], "sentinel-1-rtc": ["vv", "vh"]},
    start_date="2022-06-01",
    end_date="2022-06-15",
    edge_size=64,
    resolution=10,
    datatree=True,
    dataset=True,
)

dt["sentinel-1-rtc"]["vv"]

Tip

Keyword arguments (stackstac_kw, gee_kw, and search keyword arguments such as query) apply to all collections by default. Use collections as keys to set them per collection, e.g., query={"sentinel-2-l2a": {"eo:cloud_cover": {"lt": 10}}}.

Pixel coordinates

By default, the x and y coordinates are the top-left corners of the pixels (for both STAC and GEE). Use xy_coords="center" to get the center of the pixels instead:

da = cubo.create(
    lat=4.31,
    lon=-76.2,
    collection="sentinel-2-l2a",
    bands=["B02","B03","B04"],
    start_date="2021-06-01",
    end_date="2021-06-10",
    edge_size=64,
    resolution=10,
    xy_coords="center",
)

License

The project is licensed under the MIT license.

Citation

If you use this work, please consider citing the following paper:

@inproceedings{montero2024cubo,
  title = {On-Demand Earth System Data Cubes},
  url = {http://dx.doi.org/10.1109/IGARSS53475.2024.10640742},
  DOI = {10.1109/igarss53475.2024.10640742},
  booktitle = {IGARSS 2024 - 2024 IEEE International Geoscience and Remote Sensing Symposium},
  publisher = {IEEE},
  author = {Montero,  David and Aybar,  César and Ji,  Chaonan and Kraemer,  Guido and S\"{o}chting,  Maximilian and Teber,  Khalil and Mahecha,  Miguel D.},
  year = {2024},
  month = jul,
  pages = {7529–7532}
}

Logo Attribution

The logo and images were created using dice icons created by Freepik - Flaticon.