# Create a data cube

**Task.** How do images taken at irregular times become a regular cube?

A regular data cube made in your browser by [gdalcubes.py](https://github.com/spatial-tools/gdalcubes.py) compiled to WebAssembly. A quick flavour of how images taken on irregular dates become one image per regular interval, from live Sentinel-2 L2A or Landsat 8/9 L2 imagery. Select an area on the map, choose the settings and press Regularize; nothing is fetched before.

Journey: [Building data cubes from image collections, step 4: What does sits_regularize() do?](/regular-cubes#regularize).

**Data.** Computed in your browser by gdalcubes.py from live Planetary Computer imagery.

## The gdalcubes.py call

```python
collection = gc.stac_image_collection(
    items, assets=["B04", "B03", "B02", "B08", "SCL"], out="collection.db"
)
view = gc.cube_view(
    srs="EPSG:32723",
    left=186080, right=189360, bottom=8247660, top=8251060,
    dx=20, dy=20,
    t0="2024-06-01", t1="2024-07-31",
    dt="P1M",
    aggregation="median",
    resampling="bilinear",
)
mask = {"band": "SCL", "values": [0, 1, 2, 3, 8, 9, 10]}
cube = gc.raster_cube(collection, view, chunks=(1, 256, 256), mask=mask, strict=True)
cube = gc.select_bands(cube, ["B04", "B03", "B02", "B08"])
gc.write_netcdf(cube, "cube.nc")
```

Read in gdalcubes 0.1.0.dev1, gdalcubes/cube.py:112.

The page searches Planetary Computer, reads each image's window and calibrates it to reflectance in JavaScript; gdalcubes.py then makes the regular cube from those images.

## What the page shows

No cube has run on page entry. The result reads "Regularize to see your cube." until the reader presses Regularize, which fetches the imagery from Planetary Computer and computes the cube in the browser. The area is chosen on the map with Select area; without a map the default area below is the one that runs.

## The default request

| Setting | Value |
|---|---|
| Area (west, south, east, north, degrees) | -47.93, -15.83, -47.9, -15.8 |
| Collection | Sentinel-2 L2A |
| Bands | B04 red, B03 green, B02 blue, B08 near infrared · mask SCL |
| Dates | 2024-06-01 to 2024-07-31 |
| Interval | P1M |
| Pixel size | 20 m |
| Combine | median |
| Clouds and shadows removed | yes |

The default grid: EPSG:32723, 164 × 170 cells of 20 m (27,880 cells), from 186080 to 189360 m east and 8247660 to 8251060 m north.

| Start | End | Beyond the dates asked |
|---|---|---|
| 2024-06-01 | 2024-06-30 | no |
| 2024-07-01 | 2024-07-31 | no |

## What to try

1. Press Regularize on the default request: 164 × 170 cells of 20 m · 2 intervals of one month, each interval drawn on the map.
2. Switch the collection to Landsat 8/9 L2: the pixel size becomes 30 m and the plan 110 × 113 cells of 30 m · 2 intervals of one month.
3. Choose 16 days (P16D): the same dates make 164 × 170 cells of 20 m · 4 intervals of 16 days.
4. After a run, switch True colour to Near infrared: the picture changes without another run.
5. After a run, change a setting: the cube stays, marked "Settings changed since this cube". An area wider than 0.3° is refused before any request.

## Limits of a job

| Limit | At most |
|---|---|
| Acquisitions read | 50 |
| Cells per interval | 250,000 |
| Intervals | 32 |
| Size of the cube's file | 64.0 MB |
| Size of the images read, staged | 72.0 MB |
| Estimated peak memory | 384.0 MB |
| Span of the area | 0.3° |
| Pages of the search | 5 |

## Fixed for every 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

```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](https://planetarycomputer.microsoft.com/): the STAC catalogue and the imagery.
- [gdalcubes.py](https://pypi.org/project/gdalcubes/0.1.0.dev1/): 0.1.0.dev1, its WebAssembly wheel.
- [Pyodide](https://pyodide.org/): 314.0.7, Python in the browser.
- Basemap: OpenFreeMap © OpenMapTiles Data from OpenStreetMap ([licence](https://www.openstreetmap.org/copyright)).
- Basemap: EOxCloudless https://cloudless.eox.at by EOX IT Services GmbH (Contains modified Copernicus Sentinel data 2016) ([licence](https://creativecommons.org/licenses/by/4.0/)).

Also in the journeys: [Building data cubes from image collections, step 3: What must a regular cube agree on?](/regular-cubes#grid).

In the sits book: [Building regular data cubes](https://e-sensing.github.io/sitsbook/dc_regularize.html).