Code Snippets

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 provides a snapshot of the Hub, the pyeodh API client and the various datasets as of November 2024.

Author(s): Alastair Graham

Date created: 2024-11-08

Date last modified: 2024-11-08

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

Copyright (c) , 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.

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 (as structured in November 2024). The STAC catalogue may be reprocessed in future to allow the default settings for odc-stac to be implemented.

# Note that this uses pystac and pystac_client

import os

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

client = dask.distributed.Client()

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

As of November 2024, the STAC items in sentinel2_ard are missing information because the projection and raster extensions are not installed. This means that odc-stac will 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 those extensions.

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

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

<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)

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

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