Using Dask to assess landcover¶

Technical

Description & purpose: This notebook introduces dask and how it can be used on the Hub to process some larger jobs. As a demonstrator we will calculate the modal landcover class per pixel over a decade of landcover data from the European Space Agency (ESA). The data collection on EODH that we will use is the ESA Land Cover Climate Change Initiative (CCI): Global Land Cover Maps. Much of the workflow is expressed as lazy Dask graphs and are only computed when .compute() is called. Although this example is simple, it demonstrates how Dask can parallelise across CPU cores or a distributed cluster with zero code changes.

First we need to set up our Jupyter environment.

In [ ]:

# Run this cell if pyeodh is not installed, or needs updating
%pip install --upgrade pyeodh
%pip install rioxarray xarray numpy pandas pystac matplotlib

# Run this cell if pyeodh is not installed, or needs updating %pip install --upgrade pyeodh %pip install rioxarray xarray numpy pandas pystac matplotlib

In [2]:

# Import packages
import os
import warnings

import pyeodh
import pystac
import urllib.request

import numpy as np
import pandas as pd
import xarray as xr
import dask
import dask.array as da
import matplotlib.pyplot as plt
import rioxarray

%matplotlib inline
warnings.filterwarnings("ignore", category=UserWarning)

# Import packages import os import warnings import pyeodh import pystac import urllib.request import numpy as np import pandas as pd import xarray as xr import dask import dask.array as da import matplotlib.pyplot as plt import rioxarray %matplotlib inline warnings.filterwarnings("ignore", category=UserWarning)

This notebook is written and constructed in a way that is modular. Some of the inputs can be configured in the next cell. The main processing steps are then presented as a series of functions - those functions can be easily changed, run per cell or chained together and run in a final cell.

In [3]:

# CONFIGURATION

START_YEAR = 2006
END_YEAR   = 2015   # inclusive 

DATETIME_RANGE = f"{START_YEAR}-01-01T00:00:00Z/{END_YEAR}-12-31T23:59:59Z"

# Bounding box e.g. mainland Britain in EPSG:4326 [west, south, east, north]
# (lon_min, lat_min, lon_max, lat_max) can be defined using https://geojson.io/

# Area around Oxford
BBOX = [-1.4724114697858113, 51.606628996417555, -1.0989795612960052, 51.84678057888078]

COLLECTION_ID = "land_cover"   # ESA Land Cover Climate Change Initiative (CCI): Global Land Cover Maps

# CONFIGURATION START_YEAR = 2006 END_YEAR = 2015 # inclusive DATETIME_RANGE = f"{START_YEAR}-01-01T00:00:00Z/{END_YEAR}-12-31T23:59:59Z" # Bounding box e.g. mainland Britain in EPSG:4326 [west, south, east, north] # (lon_min, lat_min, lon_max, lat_max) can be defined using https://geojson.io/ # Area around Oxford BBOX = [-1.4724114697858113, 51.606628996417555, -1.0989795612960052, 51.84678057888078] COLLECTION_ID = "land_cover" # ESA Land Cover Climate Change Initiative (CCI): Global Land Cover Maps

Find the data¶

First we search EODH using pyeodh for the ESA Land Cover data.

In [5]:

# Find the landcover items on EODH

def search_landcover_items(collection = COLLECTION_ID, bbox = BBOX, datetime_range = DATETIME_RANGE):
    """
    Use pyeodh to search the EODH STAC catalogue and return a list of
    item objects matching the query.
    """

    client = pyeodh.Client(
        base_url="https://eodatahub.org.uk"
    ).get_catalog_service()
    
    print(f" Searching collection '{collection}' …")
    print(f" bbox          : {bbox}")
    print(f" datetime range: {datetime_range}")

    # get the items
    items = client.search(
        collections=[collection],
        bbox=bbox,
        datetime=datetime_range,
    )

    # Count the number of items returned by the search 
    total_items = sum(1 for _ in items)
    print(f"Total items: {total_items}")

    return items

items = search_landcover_items()

# Find the landcover items on EODH def search_landcover_items(collection = COLLECTION_ID, bbox = BBOX, datetime_range = DATETIME_RANGE): """ Use pyeodh to search the EODH STAC catalogue and return a list of item objects matching the query. """ client = pyeodh.Client( base_url="https://eodatahub.org.uk" ).get_catalog_service() print(f" Searching collection '{collection}' …") print(f" bbox : {bbox}") print(f" datetime range: {datetime_range}") # get the items items = client.search( collections=[collection], bbox=bbox, datetime=datetime_range, ) # Count the number of items returned by the search total_items = sum(1 for _ in items) print(f"Total items: {total_items}") return items items = search_landcover_items()

 Searching collection 'land_cover' …
 bbox          : [-1.4724114697858113, 51.606628996417555, -1.0989795612960052, 51.84678057888078]
 datetime range: 2006-01-01T00:00:00Z/2015-12-31T23:59:59Z
Total items: 10

In [6]:

# Get details about the items

def inspect_items(items):
    """
    Print a summary of STAC items with one row per item.
    """
    rows = []
    for item in items:
        dt = item.datetime
        assets = list(item.assets.keys())
        rows.append({
            "id":       item.id,
            "datetime": dt,
            "assets":   ", ".join(assets),
        })

    df = pd.DataFrame(rows).sort_values("datetime").reset_index(drop=True)
    print("STAC item summary:")
    print(df.to_string(index=False))
    print()
    return df

summary = inspect_items(items)

# Determine which asset to load
asset_key = 'GeoTIFF'

# Get details about the items def inspect_items(items): """ Print a summary of STAC items with one row per item. """ rows = [] for item in items: dt = item.datetime assets = list(item.assets.keys()) rows.append({ "id": item.id, "datetime": dt, "assets": ", ".join(assets), }) df = pd.DataFrame(rows).sort_values("datetime").reset_index(drop=True) print("STAC item summary:") print(df.to_string(index=False)) print() return df summary = inspect_items(items) # Determine which asset to load asset_key = 'GeoTIFF'

STAC item summary:
                                        id                  datetime          assets
ESACCI-LC-L4-LCCS-Map-300m-P1Y-2006-v2.0.7 2006-07-02 00:00:00+00:00 NetCDF, GeoTIFF
ESACCI-LC-L4-LCCS-Map-300m-P1Y-2007-v2.0.7 2007-07-02 00:00:00+00:00 GeoTIFF, NetCDF
ESACCI-LC-L4-LCCS-Map-300m-P1Y-2008-v2.0.7 2008-07-01 12:00:00+00:00 GeoTIFF, NetCDF
ESACCI-LC-L4-LCCS-Map-300m-P1Y-2009-v2.0.7 2009-07-02 00:00:00+00:00 GeoTIFF, NetCDF
ESACCI-LC-L4-LCCS-Map-300m-P1Y-2010-v2.0.7 2010-07-02 00:00:00+00:00 GeoTIFF, NetCDF
ESACCI-LC-L4-LCCS-Map-300m-P1Y-2011-v2.0.7 2011-07-02 00:00:00+00:00 GeoTIFF, NetCDF
ESACCI-LC-L4-LCCS-Map-300m-P1Y-2012-v2.0.7 2012-07-01 12:00:00+00:00 NetCDF, GeoTIFF
ESACCI-LC-L4-LCCS-Map-300m-P1Y-2013-v2.0.7 2013-07-02 00:00:00+00:00 GeoTIFF, NetCDF
ESACCI-LC-L4-LCCS-Map-300m-P1Y-2014-v2.0.7 2014-07-02 00:00:00+00:00 NetCDF, GeoTIFF
ESACCI-LC-L4-LCCS-Map-300m-P1Y-2015-v2.0.7 2015-07-02 00:00:00+00:00 NetCDF, GeoTIFF

How-to Highlight

Take the paginated responses and extract the information we want to create a list.

In [7]:

def to_pystac_items(items):
    pystac_items = []
    for item in items:
        raw = item.to_dict() # dict with 'type', 'geometry', 'properties', 'assets', etc.
        pystac_item = pystac.Item.from_dict(raw)
        pystac_items.append(pystac_item)

    # Sort chronologically so the time axis is ordered
    pystac_items.sort(key=lambda i: i.datetime or i.properties.get("datetime", ""))
    return pystac_items
a = to_pystac_items(items)

def to_pystac_items(items): pystac_items = [] for item in items: raw = item.to_dict() # dict with 'type', 'geometry', 'properties', 'assets', etc. pystac_item = pystac.Item.from_dict(raw) pystac_items.append(pystac_item) # Sort chronologically so the time axis is ordered pystac_items.sort(key=lambda i: i.datetime or i.properties.get("datetime", "")) return pystac_items a = to_pystac_items(items)

In [8]:

# check that the first three items have the info we need
for item in a[:3]:
        print(f"\n Item: {item.id}  |  datetime: {item.datetime}")
        print(f"  {'Asset key':<20} {'media_type':<20} {'href[:60]'}")
        print(f"  {'-'*20} {'-'*20} {'-'*60}")
        for key, asset in item.assets.items():
            print(f"  {key:<20} {str(asset.media_type):<20} {asset.href[:60]}")

# check that the first three items have the info we need for item in a[:3]: print(f"\n Item: {item.id} | datetime: {item.datetime}") print(f" {'Asset key':<20} {'media_type':<20} {'href[:60]'}") print(f" {'-'*20} {'-'*20} {'-'*60}") for key, asset in item.assets.items(): print(f" {key:<20} {str(asset.media_type):<20} {asset.href[:60]}")

 Item: ESACCI-LC-L4-LCCS-Map-300m-P1Y-2006-v2.0.7  |  datetime: 2006-07-02 00:00:00+00:00
  Asset key            media_type           href[:60]
  -------------------- -------------------- ------------------------------------------------------------
  NetCDF               None                 https://dap.ceda.ac.uk/neodc/esacci/land_cover/data/land_cov
  GeoTIFF              None                 https://dap.ceda.ac.uk/neodc/esacci/land_cover/data/land_cov

 Item: ESACCI-LC-L4-LCCS-Map-300m-P1Y-2007-v2.0.7  |  datetime: 2007-07-02 00:00:00+00:00
  Asset key            media_type           href[:60]
  -------------------- -------------------- ------------------------------------------------------------
  GeoTIFF              None                 https://dap.ceda.ac.uk/neodc/esacci/land_cover/data/land_cov
  NetCDF               None                 https://dap.ceda.ac.uk/neodc/esacci/land_cover/data/land_cov

 Item: ESACCI-LC-L4-LCCS-Map-300m-P1Y-2008-v2.0.7  |  datetime: 2008-07-01 12:00:00+00:00
  Asset key            media_type           href[:60]
  -------------------- -------------------- ------------------------------------------------------------
  GeoTIFF              None                 https://dap.ceda.ac.uk/neodc/esacci/land_cover/data/land_cov
  NetCDF               None                 https://dap.ceda.ac.uk/neodc/esacci/land_cover/data/land_cov

Build a data cube¶

The following section sets up a data cube of just the area of interest.

In [9]:

# Extract one URL per item (defaults to GeoTIFF)
pystac_items = a

def get_asset_url(item) -> tuple[str, str]:
    """Return (url, format) preferring GeoTIFF over NetCDF."""
    assets = item.assets
    for key in ("GeoTIFF", "geotiff", "data", "classification", "lcm"):
        if key in assets:
            return assets[key].href, "tiff"
    for key in ("NetCDF", "netcdf", "nc"):
        if key in assets:
            return assets[key].href, "netcdf"
    # fallback: first asset
    first = next(iter(assets.values()))
    return first.href, "unknown"

# Extract one URL per item (defaults to GeoTIFF) pystac_items = a def get_asset_url(item) -> tuple[str, str]: """Return (url, format) preferring GeoTIFF over NetCDF.""" assets = item.assets for key in ("GeoTIFF", "geotiff", "data", "classification", "lcm"): if key in assets: return assets[key].href, "tiff" for key in ("NetCDF", "netcdf", "nc"): if key in assets: return assets[key].href, "netcdf" # fallback: first asset first = next(iter(assets.values())) return first.href, "unknown"

In [10]:

# Open a single item. This is used as input for the cube.

def open_item(item, override_path):
    """
    Open one STAC item as a lazy (Dask-backed) 2-D DataArray (y, x).
    """
    if override_path:
        url = override_path
        fmt = "netcdf" if override_path.endswith(".nc") else "tiff"
    else:
        url, fmt = get_asset_url(item)
    
    print(f"  Opening [{fmt}] {item.id}  →  {url[:80]}")

    if fmt == "tiff":
        da_item = xr.open_dataarray(
            url,
            engine="rasterio",
            chunks={"x": 512, "y": 512},   # lazy Dask chunks
            mask_and_scale=False,          # keep raw integer class labels
        )

        if "band" in da_item.dims:
            da_item = da_item.squeeze("band", drop=True)

    else:
        # NetCDF: open the first data variable that looks like landcover
        ds = xr.open_dataset(
            url,
            engine="netcdf4",
            chunks={"lat": 512, "lon": 512},
            mask_and_scale=False,
        )

        lc_var = next(
            (v for v in ds.data_vars
             if any(k in v.lower() for k in ("lccs", "lc_class", "class", "land"))),
            list(ds.data_vars)[0],          # fallback: first variable
        )
        print(f"    Using variable: '{lc_var}'")
        da_item = ds[lc_var]
        # Normalise spatial dim names to y/x
        rename = {}
        if "lat" in da_item.dims:  rename["lat"] = "y"
        if "lon" in da_item.dims:  rename["lon"] = "x"
        if rename:
            da_item = da_item.rename(rename)
        # Drop any singleton time dim already present in the file
        if "time" in da_item.dims:
            da_item = da_item.squeeze("time", drop=True)

    # Tag with this item's datetime so concat builds a proper time axis
    da_item = da_item.assign_coords(time=item.datetime).expand_dims("time")
    return da_item

# Open a single item. This is used as input for the cube. def open_item(item, override_path): """ Open one STAC item as a lazy (Dask-backed) 2-D DataArray (y, x). """ if override_path: url = override_path fmt = "netcdf" if override_path.endswith(".nc") else "tiff" else: url, fmt = get_asset_url(item) print(f" Opening [{fmt}] {item.id} → {url[:80]}") if fmt == "tiff": da_item = xr.open_dataarray( url, engine="rasterio", chunks={"x": 512, "y": 512}, # lazy Dask chunks mask_and_scale=False, # keep raw integer class labels ) if "band" in da_item.dims: da_item = da_item.squeeze("band", drop=True) else: # NetCDF: open the first data variable that looks like landcover ds = xr.open_dataset( url, engine="netcdf4", chunks={"lat": 512, "lon": 512}, mask_and_scale=False, ) lc_var = next( (v for v in ds.data_vars if any(k in v.lower() for k in ("lccs", "lc_class", "class", "land"))), list(ds.data_vars)[0], # fallback: first variable ) print(f" Using variable: '{lc_var}'") da_item = ds[lc_var] # Normalise spatial dim names to y/x rename = {} if "lat" in da_item.dims: rename["lat"] = "y" if "lon" in da_item.dims: rename["lon"] = "x" if rename: da_item = da_item.rename(rename) # Drop any singleton time dim already present in the file if "time" in da_item.dims: da_item = da_item.squeeze("time", drop=True) # Tag with this item's datetime so concat builds a proper time axis da_item = da_item.assign_coords(time=item.datetime).expand_dims("time") return da_item

In [11]:

def download_items(pystac_items, out_dir="./lc_data"):
    '''
        Download dataset locally - only run as a back up and if the dataset is small 
        e.g. we have 10 landcover items
    '''
    os.makedirs(out_dir, exist_ok=True)
    local_paths = []
    for item in pystac_items:
        url, fmt = get_asset_url(item)
        ext = "nc" if fmt == "netcdf" else "tif"
        local_path = os.path.join(out_dir, f"{item.id}.{ext}")
        if not os.path.exists(local_path):
            print(f"Downloading {item.id} …")
            urllib.request.urlretrieve(url, local_path)
        else:
            print(f"  Already cached: {item.id}.{ext}")
        local_paths.append(local_path)
    return local_paths


def build_cube(pystac_items, local_paths):
    '''
        Build the data cube from all of the data
    '''
    print(f"\n Opening {len(pystac_items)} items lazily …")
    
    if local_paths:
        # Use local files — no network I/O during Dask compute
        slices = [open_item(item, override_path=path) 
                  for item, path in zip(pystac_items, local_paths)]
    else:
        # Read direct from remote URLs
        slices = [open_item(item) for item in pystac_items]

    print("\n Concatenating along time axis …")
    cube = xr.concat(slices, dim="time").sortby("time")
    cube = cube.rename("landcover")
    print(f"      Cube shape : {dict(cube.sizes)}")
    print(f"      Time range : {cube.time.values[0]}  →  {cube.time.values[-1]}")
    print(f"      Dask chunks: {cube.chunks}")
    return cube

def download_items(pystac_items, out_dir="./lc_data"): ''' Download dataset locally - only run as a back up and if the dataset is small e.g. we have 10 landcover items ''' os.makedirs(out_dir, exist_ok=True) local_paths = [] for item in pystac_items: url, fmt = get_asset_url(item) ext = "nc" if fmt == "netcdf" else "tif" local_path = os.path.join(out_dir, f"{item.id}.{ext}") if not os.path.exists(local_path): print(f"Downloading {item.id} …") urllib.request.urlretrieve(url, local_path) else: print(f" Already cached: {item.id}.{ext}") local_paths.append(local_path) return local_paths def build_cube(pystac_items, local_paths): ''' Build the data cube from all of the data ''' print(f"\n Opening {len(pystac_items)} items lazily …") if local_paths: # Use local files — no network I/O during Dask compute slices = [open_item(item, override_path=path) for item, path in zip(pystac_items, local_paths)] else: # Read direct from remote URLs slices = [open_item(item) for item in pystac_items] print("\n Concatenating along time axis …") cube = xr.concat(slices, dim="time").sortby("time") cube = cube.rename("landcover") print(f" Cube shape : {dict(cube.sizes)}") print(f" Time range : {cube.time.values[0]} → {cube.time.values[-1]}") print(f" Dask chunks: {cube.chunks}") return cube

In [12]:

# only run on manageable datasets
local_paths = download_items(pystac_items, out_dir="./lc_data")

# only run on manageable datasets local_paths = download_items(pystac_items, out_dir="./lc_data")

  Already cached: ESACCI-LC-L4-LCCS-Map-300m-P1Y-2006-v2.0.7.tif
  Already cached: ESACCI-LC-L4-LCCS-Map-300m-P1Y-2007-v2.0.7.tif
  Already cached: ESACCI-LC-L4-LCCS-Map-300m-P1Y-2008-v2.0.7.tif
  Already cached: ESACCI-LC-L4-LCCS-Map-300m-P1Y-2009-v2.0.7.tif
  Already cached: ESACCI-LC-L4-LCCS-Map-300m-P1Y-2010-v2.0.7.tif
  Already cached: ESACCI-LC-L4-LCCS-Map-300m-P1Y-2011-v2.0.7.tif
  Already cached: ESACCI-LC-L4-LCCS-Map-300m-P1Y-2012-v2.0.7.tif
  Already cached: ESACCI-LC-L4-LCCS-Map-300m-P1Y-2013-v2.0.7.tif
  Already cached: ESACCI-LC-L4-LCCS-Map-300m-P1Y-2014-v2.0.7.tif
  Already cached: ESACCI-LC-L4-LCCS-Map-300m-P1Y-2015-v2.0.7.tif

In [13]:

#Build the actual data cube
cube = build_cube(pystac_items, local_paths=local_paths)

#Build the actual data cube cube = build_cube(pystac_items, local_paths=local_paths)

 Opening 10 items lazily …
  Opening [tiff] ESACCI-LC-L4-LCCS-Map-300m-P1Y-2006-v2.0.7  →  ./lc_data/ESACCI-LC-L4-LCCS-Map-300m-P1Y-2006-v2.0.7.tif
  Opening [tiff] ESACCI-LC-L4-LCCS-Map-300m-P1Y-2007-v2.0.7  →  ./lc_data/ESACCI-LC-L4-LCCS-Map-300m-P1Y-2007-v2.0.7.tif
  Opening [tiff] ESACCI-LC-L4-LCCS-Map-300m-P1Y-2008-v2.0.7  →  ./lc_data/ESACCI-LC-L4-LCCS-Map-300m-P1Y-2008-v2.0.7.tif
  Opening [tiff] ESACCI-LC-L4-LCCS-Map-300m-P1Y-2009-v2.0.7  →  ./lc_data/ESACCI-LC-L4-LCCS-Map-300m-P1Y-2009-v2.0.7.tif
  Opening [tiff] ESACCI-LC-L4-LCCS-Map-300m-P1Y-2010-v2.0.7  →  ./lc_data/ESACCI-LC-L4-LCCS-Map-300m-P1Y-2010-v2.0.7.tif
  Opening [tiff] ESACCI-LC-L4-LCCS-Map-300m-P1Y-2011-v2.0.7  →  ./lc_data/ESACCI-LC-L4-LCCS-Map-300m-P1Y-2011-v2.0.7.tif
  Opening [tiff] ESACCI-LC-L4-LCCS-Map-300m-P1Y-2012-v2.0.7  →  ./lc_data/ESACCI-LC-L4-LCCS-Map-300m-P1Y-2012-v2.0.7.tif
  Opening [tiff] ESACCI-LC-L4-LCCS-Map-300m-P1Y-2013-v2.0.7  →  ./lc_data/ESACCI-LC-L4-LCCS-Map-300m-P1Y-2013-v2.0.7.tif
  Opening [tiff] ESACCI-LC-L4-LCCS-Map-300m-P1Y-2014-v2.0.7  →  ./lc_data/ESACCI-LC-L4-LCCS-Map-300m-P1Y-2014-v2.0.7.tif
  Opening [tiff] ESACCI-LC-L4-LCCS-Map-300m-P1Y-2015-v2.0.7  →  ./lc_data/ESACCI-LC-L4-LCCS-Map-300m-P1Y-2015-v2.0.7.tif

 Concatenating along time axis …
      Cube shape : {'time': 10, 'y': 64800, 'x': 129600}
      Time range : 1151798400000000000  →  1435795200000000000
      Dask chunks: ((1, 1, 1, 1, 1, 1, 1, 1, 1, 1), (512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 288), (512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 512, 64))

In [14]:

cube

cube

Out[14]:

<xarray.DataArray 'landcover' (time: 10, y: 64800, x: 129600)> Size: 84GB
dask.array<getitem, shape=(10, 64800, 129600), dtype=uint8, chunksize=(1, 512, 512), chunktype=numpy.ndarray>
Coordinates:
  * x            (x) float64 1MB -180.0 -180.0 -180.0 ... 180.0 180.0 180.0
  * y            (y) float64 518kB 90.0 90.0 89.99 89.99 ... -89.99 -90.0 -90.0
    spatial_ref  int64 8B 0
  * time         (time) object 80B 1151798400000000000 ... 1435795200000000000
Attributes:
    Copyright:       ESA 2017 - UCLouvain
    Dataset:         Global annual land cover map at 300 m based on the full ...
    Description:     Pixel value corresponds to the label of a land cover cla...
    Scaling Factor:  none
    Version:         2.0.7
    AREA_OR_POINT:   Area
    _FillValue:      0
    scale_factor:    1.0
    add_offset:      0.0

xarray.DataArray

'landcover'

• time: 10
• y: 64800
• x: 129600

• dask.array<chunksize=(1, 512, 512), meta=np.ndarray>

  Array Chunk Bytes 78.21 GiB 256.00 kiB Shape (10, 64800, 129600) (1, 512, 512) Dask graph 322580 chunks in 42 graph layers Data type uint8 numpy.ndarray

• Coordinates: (4)
  • x
    (x)
    float64
    -180.0 -180.0 ... 180.0 180.0

    array([-179.998611, -179.995833, -179.993056, ...,  179.993056,  179.995833,
            179.998611], shape=(129600,))

  • y
    (y)
    float64
    90.0 90.0 89.99 ... -90.0 -90.0

    array([ 89.998611,  89.995833,  89.993056, ..., -89.993056, -89.995833,
           -89.998611], shape=(64800,))

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    grid_mapping_name :

    latitude_longitude

    spatial_ref :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    GeoTransform :

    -180.0 0.002777777777778 0.0 90.0 0.0 -0.002777777777778

    array(0)

  • time
    (time)
    object
    1151798400000000000 ... 14357952...

    array([1151798400000000000, 1183334400000000000, 1214913600000000000,
           1246492800000000000, 1278028800000000000, 1309564800000000000,
           1341144000000000000, 1372723200000000000, 1404259200000000000,
           1435795200000000000], dtype=object)

• Indexes: (3)
  • x
    PandasIndex

    PandasIndex(Index([ -179.9986111111111, -179.99583333333334, -179.99305555555554,
           -179.99027777777778, -179.98749999999998, -179.98472222222222,
           -179.98194444444442, -179.97916666666666,  -179.9763888888889,
            -179.9736111111111,
           ...
             179.9736111111399,  179.97638888891768,  179.97916666669548,
             179.9819444444732,    179.984722222251,   179.9875000000288,
             179.9902777778066,  179.99305555558433,  179.99583333336213,
            179.99861111113992],
          dtype='float64', name='x', length=129600))

  • y
    PandasIndex

    PandasIndex(Index([ 89.99861111111112,  89.99583333333334,  89.99305555555556,
            89.99027777777778,  89.98750000000001,  89.98472222222223,
            89.98194444444445,  89.97916666666667,  89.97638888888889,
            89.97361111111111,
           ...
           -89.97361111112549, -89.97638888890329, -89.97916666668105,
           -89.98194444445882, -89.98472222223661, -89.98750000001438,
           -89.99027777779217, -89.99305555556994, -89.99583333334773,
            -89.9986111111255],
          dtype='float64', name='y', length=64800))

  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2006-07-02 00:00:00+00:00', '2007-07-02 00:00:00+00:00',
                   '2008-07-01 12:00:00+00:00', '2009-07-02 00:00:00+00:00',
                   '2010-07-02 00:00:00+00:00', '2011-07-02 00:00:00+00:00',
                   '2012-07-01 12:00:00+00:00', '2013-07-02 00:00:00+00:00',
                   '2014-07-02 00:00:00+00:00', '2015-07-02 00:00:00+00:00'],
                  dtype='datetime64[ns, tzutc()]', name='time', freq=None))

• Attributes: (9)

  Copyright :

  ESA 2017 - UCLouvain

  Dataset :

  Global annual land cover map at 300 m based on the full archives of 300 m MERIS, 1 km SPOT-VEGETATION, 1 km PROBA-V and 1 km AVHRR surface reflectance 7-day composites

  Description :

  Pixel value corresponds to the label of a land cover class defined using UN-LCCS classifiers

  Scaling Factor :

  none

  Version :

  2.0.7

  AREA_OR_POINT :

  Area

  _FillValue :

  0

  scale_factor :

  1.0

  add_offset :

  0.0

In [15]:

# Clip to the bounding box

xdsc = cube.rio.clip_box(
    minx=BBOX[0],
    miny=BBOX[1],
    maxx=BBOX[2],
    maxy=BBOX[3],
)
xdsc

# Clip to the bounding box xdsc = cube.rio.clip_box( minx=BBOX[0], miny=BBOX[1], maxx=BBOX[2], maxy=BBOX[3], ) xdsc

Out[15]:

<xarray.DataArray 'landcover' (time: 10, y: 87, x: 136)> Size: 118kB
dask.array<getitem, shape=(10, 87, 136), dtype=uint8, chunksize=(1, 87, 136), chunktype=numpy.ndarray>
Coordinates:
  * x            (x) float64 1kB -1.474 -1.471 -1.468 ... -1.104 -1.101 -1.099
  * y            (y) float64 696B 51.85 51.84 51.84 51.84 ... 51.61 51.61 51.61
  * time         (time) object 80B 1151798400000000000 ... 1435795200000000000
    spatial_ref  int64 8B 0
Attributes:
    Copyright:       ESA 2017 - UCLouvain
    Dataset:         Global annual land cover map at 300 m based on the full ...
    Description:     Pixel value corresponds to the label of a land cover cla...
    Scaling Factor:  none
    Version:         2.0.7
    AREA_OR_POINT:   Area
    scale_factor:    1.0
    add_offset:      0.0
    _FillValue:      0

xarray.DataArray

'landcover'

• time: 10
• y: 87
• x: 136

• dask.array<chunksize=(1, 87, 136), meta=np.ndarray>

  Array Chunk Bytes 115.55 kiB 11.55 kiB Shape (10, 87, 136) (1, 87, 136) Dask graph 10 chunks in 43 graph layers Data type uint8 numpy.ndarray

• Coordinates: (4)
  • x
    (x)
    float64
    -1.474 -1.471 ... -1.101 -1.099

    axis :

    X

    long_name :

    longitude

    standard_name :

    longitude

    units :

    degrees_east

    array([-1.473611, -1.470833, -1.468056, -1.465278, -1.4625  , -1.459722,
           -1.456944, -1.454167, -1.451389, -1.448611, -1.445833, -1.443056,
           -1.440278, -1.4375  , -1.434722, -1.431944, -1.429167, -1.426389,
           -1.423611, -1.420833, -1.418056, -1.415278, -1.4125  , -1.409722,
           -1.406944, -1.404167, -1.401389, -1.398611, -1.395833, -1.393056,
           -1.390278, -1.3875  , -1.384722, -1.381944, -1.379167, -1.376389,
           -1.373611, -1.370833, -1.368056, -1.365278, -1.3625  , -1.359722,
           -1.356944, -1.354167, -1.351389, -1.348611, -1.345833, -1.343056,
           -1.340278, -1.3375  , -1.334722, -1.331944, -1.329167, -1.326389,
           -1.323611, -1.320833, -1.318056, -1.315278, -1.3125  , -1.309722,
           -1.306944, -1.304167, -1.301389, -1.298611, -1.295833, -1.293056,
           -1.290278, -1.2875  , -1.284722, -1.281944, -1.279167, -1.276389,
           -1.273611, -1.270833, -1.268056, -1.265278, -1.2625  , -1.259722,
           -1.256944, -1.254167, -1.251389, -1.248611, -1.245833, -1.243056,
           -1.240278, -1.2375  , -1.234722, -1.231944, -1.229167, -1.226389,
           -1.223611, -1.220833, -1.218056, -1.215278, -1.2125  , -1.209722,
           -1.206944, -1.204167, -1.201389, -1.198611, -1.195833, -1.193056,
           -1.190278, -1.1875  , -1.184722, -1.181944, -1.179167, -1.176389,
           -1.173611, -1.170833, -1.168056, -1.165278, -1.1625  , -1.159722,
           -1.156944, -1.154167, -1.151389, -1.148611, -1.145833, -1.143056,
           -1.140278, -1.1375  , -1.134722, -1.131944, -1.129167, -1.126389,
           -1.123611, -1.120833, -1.118056, -1.115278, -1.1125  , -1.109722,
           -1.106944, -1.104167, -1.101389, -1.098611])

  • y
    (y)
    float64
    51.85 51.84 51.84 ... 51.61 51.61

    axis :

    Y

    long_name :

    latitude

    standard_name :

    latitude

    units :

    degrees_north

    array([51.845833, 51.843056, 51.840278, 51.8375  , 51.834722, 51.831944,
           51.829167, 51.826389, 51.823611, 51.820833, 51.818056, 51.815278,
           51.8125  , 51.809722, 51.806944, 51.804167, 51.801389, 51.798611,
           51.795833, 51.793056, 51.790278, 51.7875  , 51.784722, 51.781944,
           51.779167, 51.776389, 51.773611, 51.770833, 51.768056, 51.765278,
           51.7625  , 51.759722, 51.756944, 51.754167, 51.751389, 51.748611,
           51.745833, 51.743056, 51.740278, 51.7375  , 51.734722, 51.731944,
           51.729167, 51.726389, 51.723611, 51.720833, 51.718056, 51.715278,
           51.7125  , 51.709722, 51.706944, 51.704167, 51.701389, 51.698611,
           51.695833, 51.693056, 51.690278, 51.6875  , 51.684722, 51.681944,
           51.679167, 51.676389, 51.673611, 51.670833, 51.668056, 51.665278,
           51.6625  , 51.659722, 51.656944, 51.654167, 51.651389, 51.648611,
           51.645833, 51.643056, 51.640278, 51.6375  , 51.634722, 51.631944,
           51.629167, 51.626389, 51.623611, 51.620833, 51.618056, 51.615278,
           51.6125  , 51.609722, 51.606944])

  • time
    (time)
    object
    1151798400000000000 ... 14357952...

    array([1151798400000000000, 1183334400000000000, 1214913600000000000,
           1246492800000000000, 1278028800000000000, 1309564800000000000,
           1341144000000000000, 1372723200000000000, 1404259200000000000,
           1435795200000000000], dtype=object)

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    grid_mapping_name :

    latitude_longitude

    spatial_ref :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    GeoTransform :

    -1.4749999999857184 0.0027777777777779882 0.0 51.84722222221918 0.0 -0.0027777777777780056

    array(0)

• Indexes: (3)
  • x
    PandasIndex

    PandasIndex(Index([-1.4736111110968295, -1.4708333333190353, -1.4680555555412695,
           -1.4652777777635038, -1.4624999999857096, -1.4597222222079438,
           -1.4569444444301496, -1.4541666666523838, -1.4513888888745896,
           -1.4486111110968238,
           ...
           -1.1236111110968068, -1.1208333333190126, -1.1180555555412468,
           -1.1152777777634526, -1.1124999999856868,  -1.109722222207921,
           -1.1069444444301269,  -1.104166666652361, -1.1013888888745669,
            -1.098611111096801],
          dtype='float64', name='x', length=136))

  • y
    PandasIndex

    PandasIndex(Index([ 51.84583333333029,  51.84305555555251,  51.84027777777473,
            51.83749999999696,  51.83472222221918,   51.8319444444414,
            51.82916666666362, 51.826388888885845, 51.823611111108065,
           51.820833333330285, 51.818055555552505,  51.81527777777473,
            51.81249999999695,  51.80972222221917,   51.8069444444414,
            51.80416666666362,  51.80138888888584,  51.79861111110806,
           51.795833333330286, 51.793055555552506, 51.790277777774726,
            51.78749999999695,  51.78472222221917,  51.78194444444139,
            51.77916666666361,  51.77638888888584,  51.77361111110806,
            51.77083333333028,  51.76805555555251,  51.76527777777473,
            51.76249999999695,  51.75972222221917, 51.756944444441395,
           51.754166666663615, 51.751388888885835, 51.748611111108055,
            51.74583333333028,   51.7430555555525,  51.74027777777472,
            51.73749999999695,  51.73472222221917,  51.73194444444139,
            51.72916666666361, 51.726388888885836, 51.723611111108056,
           51.720833333330276,   51.7180555555525,  51.71527777777472,
            51.71249999999694,  51.70972222221916,  51.70694444444139,
            51.70416666666361,  51.70138888888583,  51.69861111110806,
            51.69583333333028,   51.6930555555525,  51.69027777777472,
           51.687499999996945, 51.684722222219165, 51.681944444441385,
           51.679166666663605,  51.67638888888583,  51.67361111110805,
            51.67083333333027,   51.6680555555525,  51.66527777777472,
            51.66249999999694,  51.65972222221916, 51.656944444441386,
           51.654166666663606, 51.651388888885826,  51.64861111110805,
            51.64583333333027,  51.64305555555249,  51.64027777777471,
            51.63749999999694,  51.63472222221916,  51.63194444444138,
            51.62916666666361,  51.62638888888583,  51.62361111110805,
            51.62083333333027, 51.618055555552495, 51.615277777774715,
           51.612499999996935, 51.609722222219155,  51.60694444444138],
          dtype='float64', name='y'))

  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2006-07-02 00:00:00+00:00', '2007-07-02 00:00:00+00:00',
                   '2008-07-01 12:00:00+00:00', '2009-07-02 00:00:00+00:00',
                   '2010-07-02 00:00:00+00:00', '2011-07-02 00:00:00+00:00',
                   '2012-07-01 12:00:00+00:00', '2013-07-02 00:00:00+00:00',
                   '2014-07-02 00:00:00+00:00', '2015-07-02 00:00:00+00:00'],
                  dtype='datetime64[ns, tzutc()]', name='time', freq=None))

• Attributes: (9)

  Copyright :

  ESA 2017 - UCLouvain

  Dataset :

  Global annual land cover map at 300 m based on the full archives of 300 m MERIS, 1 km SPOT-VEGETATION, 1 km PROBA-V and 1 km AVHRR surface reflectance 7-day composites

  Description :

  Pixel value corresponds to the label of a land cover class defined using UN-LCCS classifiers

  Scaling Factor :

  none

  Version :

  2.0.7

  AREA_OR_POINT :

  Area

  scale_factor :

  1.0

  add_offset :

  0.0

  _FillValue :

  0

Analysis¶

Now that the data is in a format that can be used to calculate the land cover mode, we can use the dask .compute function.

In [16]:

# Compute dominant (modal) class per pixel

def dominant_class(cube):
    """
    For each pixel, find which landcover class appears most often
    across all years (modal value over the time axis).

    Uses Dask to compute class counts in parallel — one pass per class label.
    For ESA CCI the LCCS class values are multi-digit integers (e.g. 10, 30,
    190 …)
    """
    print("\n Computing dominant class …")

    # Sample a spatial tile to discover class labels without loading everything
    sample = cube.isel(
        y=slice(0, 200), x=slice(0, 200)
    ).compute()
    unique_classes = np.unique(sample.values)
    unique_classes = unique_classes[unique_classes > 0]   # drop nodata (0 / fill)
    
    print(f" Detected {len(unique_classes)} unique classes: {unique_classes}")

    # For each class build a (y, x) count array — all lazy Dask ops
    count_arrays = []
    for cls in unique_classes:
        count = da.sum(cube.data == int(cls), axis=0)   # sum over time axis
        count_arrays.append(count)

    # Stack to (n_classes, y, x) then argmax → index of winning class
    counts_stack = da.stack(count_arrays, axis=0)       # lazy
    argmax_idx   = da.argmax(counts_stack, axis=0)      # lazy (y, x) of ints

    # ── Trigger computation ──────────────────────────────────────────────
    # Compute argmax_idx as a plain numpy array first, then use it to look up
    # class labels via standard numpy indexing.

    print(" Triggering Dask compute …")
    
    with dask.config.set(scheduler="threads"):
        argmax_np = argmax_idx.compute()                

    # Map argmax indices to class label values using numpy
    dominant_np = unique_classes[argmax_np]

    result = xr.DataArray(
        dominant_np,
        dims=["y", "x"],
        coords={"y": cube.coords["y"], "x": cube.coords["x"]},
        attrs={"long_name": "Dominant landcover class (mode over time)"},
    )

    result = xr.DataArray(
        dominant_np,
        dims=["y", "x"],
        coords={"y": cube.coords["y"], "x": cube.coords["x"]},
        attrs={"long_name": "Dominant landcover class (mode over time)"},
    )

    print("Complete")
    return result, unique_classes
    
dominant, classes = dominant_class(xdsc)

# Compute dominant (modal) class per pixel def dominant_class(cube): """ For each pixel, find which landcover class appears most often across all years (modal value over the time axis). Uses Dask to compute class counts in parallel — one pass per class label. For ESA CCI the LCCS class values are multi-digit integers (e.g. 10, 30, 190 …) """ print("\n Computing dominant class …") # Sample a spatial tile to discover class labels without loading everything sample = cube.isel( y=slice(0, 200), x=slice(0, 200) ).compute() unique_classes = np.unique(sample.values) unique_classes = unique_classes[unique_classes > 0] # drop nodata (0 / fill) print(f" Detected {len(unique_classes)} unique classes: {unique_classes}") # For each class build a (y, x) count array — all lazy Dask ops count_arrays = [] for cls in unique_classes: count = da.sum(cube.data == int(cls), axis=0) # sum over time axis count_arrays.append(count) # Stack to (n_classes, y, x) then argmax → index of winning class counts_stack = da.stack(count_arrays, axis=0) # lazy argmax_idx = da.argmax(counts_stack, axis=0) # lazy (y, x) of ints # ── Trigger computation ────────────────────────────────────────────── # Compute argmax_idx as a plain numpy array first, then use it to look up # class labels via standard numpy indexing. print(" Triggering Dask compute …") with dask.config.set(scheduler="threads"): argmax_np = argmax_idx.compute() # Map argmax indices to class label values using numpy dominant_np = unique_classes[argmax_np] result = xr.DataArray( dominant_np, dims=["y", "x"], coords={"y": cube.coords["y"], "x": cube.coords["x"]}, attrs={"long_name": "Dominant landcover class (mode over time)"}, ) result = xr.DataArray( dominant_np, dims=["y", "x"], coords={"y": cube.coords["y"], "x": cube.coords["x"]}, attrs={"long_name": "Dominant landcover class (mode over time)"}, ) print("Complete") return result, unique_classes dominant, classes = dominant_class(xdsc)

 Computing dominant class …
 Detected 13 unique classes: [ 10  11  30  40  60  70 100 110 130 180 190 200 210]
 Triggering Dask compute …
Complete

Output summary¶

Finally, we pull out the statistics and also create a map of the dominant class for each pixel and a plot of class share.

In [17]:

# Summarise and plot

def summarise_and_plot(dominant, unique_classes):

    vals = dominant.values.ravel()
    counts = {int(cls): int(np.sum(vals == cls)) for cls in unique_classes}
    total = sum(counts.values())

    df = pd.DataFrame([
        {"class_id": k, "pixel_count": v, "pct": round(100 * v / total, 2)}
        for k, v in sorted(counts.items(), key=lambda x: -x[1])
        if v > 0
    ])

    print("\n Dominant class frequency table:")
    print(df.to_string(index=False))

    # Generate plot
    fig, (ax_map, ax_bar) = plt.subplots(1, 2, figsize=(14, 6))
    fig.suptitle("Dominant Landcover Class (10-year mode)", fontsize=13)

    # Map
    #n_classes = len(df)
    cmap = plt.colormaps.get_cmap("tab20")  
    im = ax_map.imshow(
        dominant.values,
        cmap=cmap,
        interpolation="nearest",
        aspect="auto",
    )
    plt.colorbar(im, ax=ax_map, label="Class ID", shrink=0.7)
    ax_map.set_title("Spatial map of dominant class")
    ax_map.set_xlabel("x pixel")
    ax_map.set_ylabel("y pixel")

    # Bar chart
    ax_bar.barh(
        df["class_id"].astype(str),
        df["pct"],
        color=[cmap(i) for i in range(len(df))],
        edgecolor="grey", linewidth=0.4,
    )
    ax_bar.set_xlabel("% of pixels")
    ax_bar.set_title("Share of dominant class")
    ax_bar.invert_yaxis()

    plt.tight_layout()
    plt.show()

summarise_and_plot(dominant, classes)

# Summarise and plot def summarise_and_plot(dominant, unique_classes): vals = dominant.values.ravel() counts = {int(cls): int(np.sum(vals == cls)) for cls in unique_classes} total = sum(counts.values()) df = pd.DataFrame([ {"class_id": k, "pixel_count": v, "pct": round(100 * v / total, 2)} for k, v in sorted(counts.items(), key=lambda x: -x[1]) if v > 0 ]) print("\n Dominant class frequency table:") print(df.to_string(index=False)) # Generate plot fig, (ax_map, ax_bar) = plt.subplots(1, 2, figsize=(14, 6)) fig.suptitle("Dominant Landcover Class (10-year mode)", fontsize=13) # Map #n_classes = len(df) cmap = plt.colormaps.get_cmap("tab20") im = ax_map.imshow( dominant.values, cmap=cmap, interpolation="nearest", aspect="auto", ) plt.colorbar(im, ax=ax_map, label="Class ID", shrink=0.7) ax_map.set_title("Spatial map of dominant class") ax_map.set_xlabel("x pixel") ax_map.set_ylabel("y pixel") # Bar chart ax_bar.barh( df["class_id"].astype(str), df["pct"], color=[cmap(i) for i in range(len(df))], edgecolor="grey", linewidth=0.4, ) ax_bar.set_xlabel("% of pixels") ax_bar.set_title("Share of dominant class") ax_bar.invert_yaxis() plt.tight_layout() plt.show() summarise_and_plot(dominant, classes)

 Dominant class frequency table:
 class_id  pixel_count   pct
       11         6990 59.08
      130         2211 18.69
      190         1159  9.80
       40          478  4.04
       30          351  2.97
       60          295  2.49
      210          119  1.01
      100           76  0.64
       10           48  0.41
      110           39  0.33
       70           29  0.25
      180           27  0.23
      200           10  0.08

---

Author(s): Alastair Graham

Date created: 2026-03-16

Date last modified: 2026-03-16

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.