Create a data cube

Fixed for every cube

The call above is the default cube's; the summary each run offers to download writes the call of its own cube.

Scaling

Sentinel-2 L2A: reflectance = DN × 0.0001 − 0.1 from processing baseline 04.00, DN × 0.0001 before; DN 0 and 65535 are no data. Landsat 8/9 L2: reflectance = DN × scale + offset, both read from each asset's raster:bands; DN 0 is no data. Negative reflectance is kept.

Clouds and shadows removed

Sentinel-2 L2A: SCL classes 0 no data, 1 saturated or defective, 2 cast shadow, 3 cloud shadow, 8 cloud, medium probability, 9 cloud, high probability, 10 thin cirrus. Landsat 8/9 L2: qa_pixel bits 0 fill, 1 dilated cloud, 2 cirrus, 3 cloud, 4 cloud shadow, 9 high cloud confidence, 11 high cloud shadow confidence, 15 high cirrus confidence; a cell goes when any is set.

Resampling

bilinear for reflectance, near for the mask.

Time

Monthly cubes have bins of one month, calendar months from the start's month to the end's; the others bins of 16 days from the start day, the last of which may pass the end day. A bin that reaches beyond the dates asked is marked partial; only acquisitions on the days asked go in.

Grid

The UTM zone of the area's centre, in whole cells whose edges are multiples of the pixel size (EPSG:32723 for the default area).

File

A classic netCDF file, one variable per band of dimensions (time, y, x), y from north to south, written by write_netcdf() of gdalcubes.py.

How much it loads

The first Regularize loads the engine: Pyodide 314.0.7, 13.5 MB; its packages (numpy, orjson, tqdm, micropip), 3.3 MB; the gdalcubes 0.1.0.dev1 wheel, 11.1 MB: 27.9 MB in all, uncompressed; the browser's cache may keep them for later runs. The imagery of the area comes on top, its size set by the area, the pixel size and the dates.

Opening the file in Python

"""Open a cube from Create a data cube in Python: python open_example.py cube.nc

Needs numpy and netCDF4 (pip install numpy netCDF4). The file is classic netCDF with one
variable per band, of dimensions (time, y, x), y from north to south. The time variable holds
whole numbers of the cube's step unit since its first day: "months since 2024-06-01T00:00:00"
for monthly cubes, "days since ..." for 16-day ones. netCDF4 and xarray do not turn months or
years into dates, so the dates are counted here from the calendar.
"""

from __future__ import annotations

import calendar
import sys
from datetime import datetime, timedelta

import netCDF4
import numpy as np

# units of a fixed length; months and years are counted on the calendar
FIXED_UNITS = {"seconds": 1, "minutes": 60, "hours": 3600, "days": 86400}
MONTHS_PER_UNIT = {"months": 1, "years": 12}


def _add_months(origin: datetime, months: int) -> datetime:
    """The same day and time, `months` calendar months later; refuses a day the month lacks."""
    year, month = divmod(origin.month - 1 + months, 12)
    year += origin.year
    if origin.day > calendar.monthrange(year, month + 1)[1]:
        raise ValueError(f"{origin:%Y-%m-%d} plus {months} months has no such day.")
    return origin.replace(year=year, month=month + 1)


def interval_starts(units: str, offsets: list[float]) -> list[datetime]:
    """The start of each interval, from the time variable's units and values."""
    unit, since = units.split(" since ")
    origin = datetime.fromisoformat(since)
    starts = []
    for offset in offsets:
        if offset != int(offset):
            raise ValueError(f"The time value {offset} is not a whole number of {unit}.")
        if unit in MONTHS_PER_UNIT:
            starts.append(_add_months(origin, int(offset) * MONTHS_PER_UNIT[unit]))
        elif unit in FIXED_UNITS:
            starts.append(origin + timedelta(seconds=int(offset) * FIXED_UNITS[unit]))
        else:
            raise ValueError(f"Time in {unit!r} is not handled here.")
    return starts


def open_cube(path: str) -> dict:
    """The cube's interval starts, coordinates, CRS and bands, missing values as NaN."""
    with netCDF4.Dataset(path) as dataset:
        time = dataset["time"]
        bands = {
            name: np.ma.filled(variable[:].astype("float64"), np.nan)
            for name, variable in dataset.variables.items()
            if variable.dimensions == ("time", "y", "x")
        }
        return {
            "dates": interval_starts(time.units, [float(value) for value in time[:]]),
            "x": np.asarray(dataset["x"][:]),
            "y": np.asarray(dataset["y"][:]),
            "crs": dataset["crs"].spatial_ref,
            "bands": bands,
        }


if __name__ == "__main__":
    cube = open_cube(sys.argv[1] if len(sys.argv) > 1 else "cube.nc")
    print("bands:", list(cube["bands"]), "x:", cube["x"].size, "y:", cube["y"].size)
    # the area's mean of each band at each interval: a time series
    for index, start in enumerate(cube["dates"]):
        means = {}
        for name, values in cube["bands"].items():
            valid = values[index][np.isfinite(values[index])]
            means[name] = round(float(valid.mean()), 4) if valid.size else None
        print(start.date().isoformat(), means)

A run's record

Each run's record names the search and when it ran, the acquisitions and the windows read from each file, the runtime and the file's SHA-256, so the cube can be traced. The same settings later may give another file: the catalogue and its imagery change.

Credits

• Sentinel-2 L2A: Copernicus Sentinel data, via Microsoft Planetary Computer.
• Landsat 8/9 L2: Landsat data courtesy of the U.S. Geological Survey, via Microsoft Planetary Computer.
• Microsoft Planetary Computer: the STAC catalogue and the imagery.
• gdalcubes.py: 0.1.0.dev1, its WebAssembly wheel.
• Pyodide: 314.0.7, Python in the browser.
• Basemap: OpenFreeMap © OpenMapTiles Data from OpenStreetMap (licence).
• Basemap: EOxCloudless https://cloudless.eox.at by EOX IT Services GmbH (Contains modified Copernicus Sentinel data 2016) (licence).