Data cubes¶

Technical   Satellite

Description & purpose: This Notebook is designed to showcase different tools and functions that are of scientific interest/use and implement the Earth Observation Data Hub. It is assumed that this notebook will be run within the Notebook Service on the Hub.

Create data cube using sentinel2_ard STAC records¶

Users may have a need to generate a data cube directly from the CEDA Sentinel 2 Analysis Ready Data (ARD) (sentinel2_ard). As the STAC catalogue for that dataset has been created in a specific way, users will need to be aware of this and make alterations from the default. The following code was supplied by Pete Gadomski of Development Seed and enables the generation of a data cube from the ARD catalogue. The STAC catalogue may be reprocessed in future to allow the default settings for odc-stac to be implemented.

In [ ]:

# Note that this uses pystac and pystac_client, not pyeodh

%pip install odc-stac pystac_client

# Note that this uses pystac and pystac_client, not pyeodh %pip install odc-stac pystac_client

In [ ]:

# import os

import dask.distributed
import odc.stac
import rasterio
from pystac.extensions.raster import RasterBand
from pystac_client import Client

# import os import dask.distributed import odc.stac import rasterio from pystac.extensions.raster import RasterBand from pystac_client import Client

In [2]:

client = dask.distributed.Client()

client = dask.distributed.Client()

In [3]:

client = Client.open("https://api.stac.ceda.ac.uk/")
item_collection = client.search(
    collections=["sentinel2_ard"],
    intersects={"type": "Point", "coordinates": [-1.3144, 51.5755]},
    sortby="-properties.datetime",
    max_items=10,
).item_collection()

# Commented out for clarity when submitted to training materials website.
# Uncomment the next line to view the Feature Collection

#item_collection

client = Client.open("https://api.stac.ceda.ac.uk/") item_collection = client.search( collections=["sentinel2_ard"], intersects={"type": "Point", "coordinates": [-1.3144, 51.5755]}, sortby="-properties.datetime", max_items=10, ).item_collection() # Commented out for clarity when submitted to training materials website. # Uncomment the next line to view the Feature Collection #item_collection

odc-stac may refuse to compute the "geobox" and hence the subsequent components needed to generate and view the data cube. The code in the following cell adds the required extensions to enforce this functionality.

How-to Highlight

In [ ]:

for item in item_collection.items:
    asset = item.assets["cog"]
    cog = rasterio.open(asset.href)
    epsg = cog.crs.to_epsg()
    dtypes = cog.dtypes
    shape = cog.shape
    transform = list(cog.transform)

    item.ext.add("proj")
    item.ext.add("raster")

    item.ext.proj.epsg = epsg

    cog = item.assets["cog"]
    cog.ext.raster.bands = [RasterBand.create(data_type=dtype) for dtype in dtypes]
    cog.ext.proj.shape = shape
    cog.ext.proj.transform = transform

for item in item_collection.items: asset = item.assets["cog"] cog = rasterio.open(asset.href) epsg = cog.crs.to_epsg() dtypes = cog.dtypes shape = cog.shape transform = list(cog.transform) item.ext.add("proj") item.ext.add("raster") item.ext.proj.epsg = epsg cog = item.assets["cog"] cog.ext.raster.bands = [RasterBand.create(data_type=dtype) for dtype in dtypes] cog.ext.proj.shape = shape cog.ext.proj.transform = transform

In [8]:

bbox = [-1.2, 51.6, -1.1, 51.7]
dataset = odc.stac.load(item_collection, bands=("B03", "B04", "B08"), chunks={}, bbox=bbox)
dataset

bbox = [-1.2, 51.6, -1.1, 51.7] dataset = odc.stac.load(item_collection, bands=("B03", "B04", "B08"), chunks={}, bbox=bbox) dataset

Out[8]:

<xarray.Dataset> Size: 48MB
Dimensions:      (y: 1121, x: 706, time: 10)
Coordinates:
  * y            (y) float64 9kB 2.005e+05 2.005e+05 ... 1.893e+05 1.893e+05
  * x            (x) float64 6kB 4.554e+05 4.554e+05 ... 4.624e+05 4.624e+05
    spatial_ref  int32 4B 27700
  * time         (time) datetime64[ns] 80B 2023-10-03T11:08:09 ... 2023-11-20...
Data variables:
    B03          (time, y, x) uint16 16MB dask.array<chunksize=(1, 1121, 706), meta=np.ndarray>
    B04          (time, y, x) uint16 16MB dask.array<chunksize=(1, 1121, 706), meta=np.ndarray>
    B08          (time, y, x) uint16 16MB dask.array<chunksize=(1, 1121, 706), meta=np.ndarray>

xarray.Dataset

• Dimensions:
  • y: 1121
  • x: 706
  • time: 10
• Coordinates: (4)
  • y
    (y)
    float64
    2.005e+05 2.005e+05 ... 1.893e+05

    units :

    metre

    resolution :

    -10.0

    crs :

    EPSG:27700

    array([200525., 200515., 200505., ..., 189345., 189335., 189325.])

  • x
    (x)
    float64
    4.554e+05 4.554e+05 ... 4.624e+05

    units :

    metre

    resolution :

    10.0

    crs :

    EPSG:27700

    array([455385., 455395., 455405., ..., 462415., 462425., 462435.])

  • spatial_ref
    ()
    int32
    27700

    spatial_ref :

    PROJCRS["OSGB36 / British National Grid",BASEGEOGCRS["OSGB36",DATUM["Ordnance Survey of Great Britain 1936",ELLIPSOID["Airy 1830",6377563.396,299.3249646,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4277]],CONVERSION["British National Grid",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",49,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",-2,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9996012717,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",400000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",-100000,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Engineering survey, topographic mapping."],AREA["United Kingdom (UK) - offshore to boundary of UKCS within 49°45'N to 61°N and 9°W to 2°E; onshore Great Britain (England, Wales and Scotland). Isle of Man onshore."],BBOX[49.75,-9.01,61.01,2.01]],ID["EPSG",27700]]

    crs_wkt :

    PROJCRS["OSGB36 / British National Grid",BASEGEOGCRS["OSGB36",DATUM["Ordnance Survey of Great Britain 1936",ELLIPSOID["Airy 1830",6377563.396,299.3249646,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4277]],CONVERSION["British National Grid",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",49,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",-2,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9996012717,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",400000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",-100000,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Engineering survey, topographic mapping."],AREA["United Kingdom (UK) - offshore to boundary of UKCS within 49°45'N to 61°N and 9°W to 2°E; onshore Great Britain (England, Wales and Scotland). Isle of Man onshore."],BBOX[49.75,-9.01,61.01,2.01]],ID["EPSG",27700]]

    semi_major_axis :

    6377563.396

    semi_minor_axis :

    6356256.909237285

    inverse_flattening :

    299.3249646

    reference_ellipsoid_name :

    Airy 1830

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    OSGB36

    horizontal_datum_name :

    Ordnance Survey of Great Britain 1936

    projected_crs_name :

    OSGB36 / British National Grid

    grid_mapping_name :

    transverse_mercator

    latitude_of_projection_origin :

    49.0

    longitude_of_central_meridian :

    -2.0

    false_easting :

    400000.0

    false_northing :

    -100000.0

    scale_factor_at_central_meridian :

    0.9996012717

    GeoTransform :

    455380 10 0 200530 0 -10

    array(27700, dtype=int32)

  • time
    (time)
    datetime64[ns]
    2023-10-03T11:08:09 ... 2023-11-...

    array(['2023-10-03T11:08:09.000000000', '2023-10-06T11:21:19.000000000',
           '2023-10-08T11:09:41.000000000', '2023-10-13T11:09:19.000000000',
           '2023-10-26T11:21:19.000000000', '2023-10-28T11:11:51.000000000',
           '2023-11-10T11:23:11.000000000', '2023-11-15T11:22:39.000000000',
           '2023-11-17T11:13:31.000000000', '2023-11-20T11:23:51.000000000'],
          dtype='datetime64[ns]')

• Data variables: (3)
  • B03
    (time, y, x)
    uint16
    dask.array<chunksize=(1, 1121, 706), meta=np.ndarray>

    Array Chunk Bytes 15.10 MiB 1.51 MiB Shape (10, 1121, 706) (1, 1121, 706) Dask graph 10 chunks in 3 graph layers Data type uint16 numpy.ndarray

  • B04
    (time, y, x)
    uint16
    dask.array<chunksize=(1, 1121, 706), meta=np.ndarray>

    Array Chunk Bytes 15.10 MiB 1.51 MiB Shape (10, 1121, 706) (1, 1121, 706) Dask graph 10 chunks in 3 graph layers Data type uint16 numpy.ndarray

  • B08
    (time, y, x)
    uint16
    dask.array<chunksize=(1, 1121, 706), meta=np.ndarray>

    Array Chunk Bytes 15.10 MiB 1.51 MiB Shape (10, 1121, 706) (1, 1121, 706) Dask graph 10 chunks in 3 graph layers Data type uint16 numpy.ndarray

• Indexes: (3)
  • y
    PandasIndex

    PandasIndex(Index([200525.0, 200515.0, 200505.0, 200495.0, 200485.0, 200475.0, 200465.0,
           200455.0, 200445.0, 200435.0,
           ...
           189415.0, 189405.0, 189395.0, 189385.0, 189375.0, 189365.0, 189355.0,
           189345.0, 189335.0, 189325.0],
          dtype='float64', name='y', length=1121))

  • x
    PandasIndex

    PandasIndex(Index([455385.0, 455395.0, 455405.0, 455415.0, 455425.0, 455435.0, 455445.0,
           455455.0, 455465.0, 455475.0,
           ...
           462345.0, 462355.0, 462365.0, 462375.0, 462385.0, 462395.0, 462405.0,
           462415.0, 462425.0, 462435.0],
          dtype='float64', name='x', length=706))

  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2023-10-03 11:08:09', '2023-10-06 11:21:19',
                   '2023-10-08 11:09:41', '2023-10-13 11:09:19',
                   '2023-10-26 11:21:19', '2023-10-28 11:11:51',
                   '2023-11-10 11:23:11', '2023-11-15 11:22:39',
                   '2023-11-17 11:13:31', '2023-11-20 11:23:51'],
                  dtype='datetime64[ns]', name='time', freq=None))

• Attributes: (0)

In [10]:

dataset.odc.geobox
dataset = dataset.compute()

dataset.odc.geobox dataset = dataset.compute()

The final thing to do is run a final test that all is working as expected by plotting an image.

In [11]:

_ = dataset.isel(time=0).to_array("band").plot.imshow(vmin=0, vmax=256)

_ = dataset.isel(time=0).to_array("band").plot.imshow(vmin=0, vmax=256)

---

Author(s): Alastair Graham

Date created: 2024-11-08

Date last modified: 2026-03-20

Licence: This notebook is licensed under Creative Commons Attribution-ShareAlike 4.0 International. The code is released using the BSD-2-Clause license.

Copyright © - All rights reserved.

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.