Simple load collection example: ESA WorldCover

This notebook is a simple example of how to load the ESA WorldCover dataset using the load_collection function from the openeo Python package. Here we will load the ESA WorldCover dataset for the year 2021 apply simple UDF and visualize the result.

import openeo
import json
import folium
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors


# connect to the federated backend
connection = openeo.connect("openeo.dataspace.copernicus.eu").authenticate_oidc()

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[1], line 1
----> 1 import openeo
      2 import json
      3 import folium

ModuleNotFoundError: No module named 'openeo'

List the available collections from the backend and fetch its description.

connection.list_collection_ids()

['SENTINEL3_OLCI_L1B',
 'SENTINEL3_SLSTR',
 'SENTINEL_5P_L2',
 'COPERNICUS_VEGETATION_PHENOLOGY_PRODUCTIVITY_10M_SEASON1',
 'COPERNICUS_VEGETATION_PHENOLOGY_PRODUCTIVITY_10M_SEASON2',
 'COPERNICUS_PLANT_PHENOLOGY_INDEX',
 'ESA_WORLDCOVER_10M_2020_V1',
 'ESA_WORLDCOVER_10M_2021_V2',
 'COPERNICUS_VEGETATION_INDICES',
 'SENTINEL2_L1C',
 'SENTINEL2_L2A',
 'SENTINEL1_GRD',
 'COPERNICUS_30',
 'LANDSAT8_L2',
 'SENTINEL3_SYN_L2_SYN',
 'SENTINEL3_SLSTR_L2_LST',
 'SENTINEL1_GLOBAL_MOSAICS']

connection.describe_collection("ESA_WORLDCOVER_10M_2021_V2")

# define aoi and temporal extent
def read_json(filename: str) -> dict:
    with open(filename) as input:
        field = json.load(input)
    return field


dates = ["2021-01-01", "2021-12-31"]
aoi = read_json("aoi/Antwerp.geojson")

m = folium.Map([51.18, 4.5], zoom_start=11)
folium.GeoJson(aoi).add_to(m)
m

Make this Notebook Trusted to load map: File -> Trust Notebook

# load data from the CDSE backend
cube = connection.load_collection(
    "ESA_WORLDCOVER_10M_2021_V2", temporal_extent=dates, spatial_extent=aoi
)

# let us apply a simple rescaling udf

udf = openeo.UDF(
    """
import xarray

def apply_datacube(cube: xarray.DataArray, context: dict) -> xarray.DataArray:
    cube.values = 0.8 * cube.values
    return cube
"""
)

Users can also load the UDF from a python file directly, instead of defining it as a string. For more information, please refer to the openEO UDF documentation.

rescaled_cube = cube.apply(process=udf)

rescaled_cube

# job configuration

job_options = {
    "executor-memory": "2G",
    "executor-memoryOverhead": "3G",
    "python-memory": "3G",
}
# running a batch job

rescaled_cube.execute_batch(
    title="ESA_WORLDCOVER_10M_2021_V2",
    outputfile="ESAWorldCover.nc",
    job_options=job_options,
)

0:00:00 Job 'j-2502101653094c4694b9c4f438184002': send 'start'
0:00:13 Job 'j-2502101653094c4694b9c4f438184002': created (progress 0%)
0:00:18 Job 'j-2502101653094c4694b9c4f438184002': created (progress 0%)
0:00:25 Job 'j-2502101653094c4694b9c4f438184002': created (progress 0%)
0:00:33 Job 'j-2502101653094c4694b9c4f438184002': created (progress 0%)
0:00:43 Job 'j-2502101653094c4694b9c4f438184002': created (progress 0%)
0:00:55 Job 'j-2502101653094c4694b9c4f438184002': created (progress 0%)
0:01:10 Job 'j-2502101653094c4694b9c4f438184002': running (progress N/A)
0:01:29 Job 'j-2502101653094c4694b9c4f438184002': running (progress N/A)
0:01:53 Job 'j-2502101653094c4694b9c4f438184002': running (progress N/A)
0:02:23 Job 'j-2502101653094c4694b9c4f438184002': finished (progress 100%)

# visualize the result

ds = xr.open_dataset("ESAWorldCover.nc")
data = ds[["MAP"]].to_array(dim="bands").squeeze()
# Replace zeros with NaN
arr = data.values
arr[arr == 0] = np.nan

# Generate lat/lon grid
lon, lat = np.meshgrid(
    data.x.values.astype(np.float64), data.y.values.astype(np.float64)
)

# Normalize data using a thematic color scale
cmap = plt.get_cmap("hsv")  # Choose a thematic colormap
norm = mcolors.Normalize(vmin=np.nanmin(arr), vmax=np.nanmax(arr))  # Scale values
colored_data = cmap(norm(arr))  # Convert data to RGBA colors

# Ensure correct shape (H, W, 4)
colored_data = (colored_data[:, :, :4] * 255).astype(
    np.uint8
)  # Convert to uint8 format

# Create Folium map
m = folium.Map(location=[lat.mean(), lon.mean()], zoom_start=12)

folium.raster_layers.ImageOverlay(
    image=colored_data,
    bounds=[[lat.min(), lon.min()], [lat.max(), lon.max()]],
    mercator_project=True,
    opacity=0.7,
).add_to(m)

m

Make this Notebook Trusted to load map: File -> Trust Notebook