Climate DT Parameter - Data Access using DEDL HDA

pip install --user --quiet --upgrade destinelab

import destinelab as deauth
import json
import importlib.metadata

import requests
from requests.adapters import HTTPAdapter
from urllib3.util.retry import Retry
import os
from getpass import getpass
from tqdm import tqdm
import time
from datetime import datetime
from urllib.parse import unquote
from IPython.display import JSON
import ipywidgets as w

DESP_USERNAME = input("Please input your DESP username: ")
DESP_PASSWORD = getpass("Please input your DESP password: ")

auth = deauth.AuthHandler(DESP_USERNAME, DESP_PASSWORD)
access_token = auth.get_token()
if access_token is not None:
    print("DEDL/DESP Access Token Obtained Successfully")
else:
    print("Failed to Obtain DEDL/DESP Access Token")

auth_headers = {"Authorization": f"Bearer {access_token}"}

auth.is_DTaccess_allowed(access_token)

COLLECTION_ID="EO.ECMWF.DAT.D1.DT_CLIMATE.G1.SCENARIOMIP_SSP3-7.0_IFS-NEMO.R1"

HDA_STAC_ENDPOINT="https://hda.data.destination-earth.eu/stac/v2"
print("STAC endpoint: ", HDA_STAC_ENDPOINT)

HDA_DISCOVERY_ENDPOINT = HDA_STAC_ENDPOINT+'/collections'
print("HDA discovery endpoint: ", HDA_DISCOVERY_ENDPOINT)

discovery_json = requests.get(HDA_DISCOVERY_ENDPOINT,params = 
                             {"q": '"Climate Change Adaptation Digital Twin","Future Projection","IFS-NEMO"'}).json()

print("\nThe discovery result give us:\nthe collection ID : ", json.dumps(discovery_json["collections"][0].get("id")))
print("\nIts time extension : ", json.dumps(discovery_json["collections"][0].get("extent").get("temporal").get("interval")))

print("\nThe available parameters: ")
print("(Type in the text box to narrow the list. When only one item remains, its details will appear)")

parameters=discovery_json["collections"][0]["cube:variables"]
keys = sorted(parameters)
txt = w.Text(description="")
out = w.Output()
def run(_):
    q = txt.value.lower()
    with out:
        out.clear_output(); 
        matches = [k for k in keys if q in k.lower()]
        for k in matches:
            print(k)
            if len(matches)==1:            
                var = discovery_json["collections"][0]["cube:variables"][k]
                print(var["description"], "\n")
                print(json.dumps(var["attrs"], indent=2))
                #print(json.dumps(discovery_json["collections"][0]["cube:variables"][k], indent=2))

                
txt.observe(run, names="value"); display(txt, out); run(None)

#print("\nThe information related to each parameter, e.g. the 2 metre temperature: ")
#JSON(discovery_json["collections"][0].get("cube:variables").get("2_metre_temperature(sfc)"))

components_of_wind=[]
for var_name, var_info in parameters.items():
    if var_name.endswith("component_of_wind(pl)"):
        var_type = var_info.get("type")
        var_unit = var_info.get("unit")
        var_url = var_info.get("attrs").get("url")

        components_of_wind.append ( {
            "type": var_type,
            "unit": var_unit,
            "long_name": var_info.get("attrs").get("long_name"),
            "shortName": var_info.get("attrs").get("shortName"),
            "standard_name": var_info.get("attrs").get("standard_name"),
            "url": var_url,
            "parameter_ID": var_info.get("attrs").get("parameter_ID"),
            "product_type": var_info.get("attrs").get("product_type"),
            "levtype": var_info.get("attrs").get("levtype"),
            "levelist": var_info.get("attrs").get("levelist"),
            "time": var_info.get("attrs").get("time")
        })
components_of_wind

response = requests.post(HDA_STAC_ENDPOINT+"/search", headers=auth_headers, json={
 "collections": [COLLECTION_ID],
    "datetime": "2028-06-10T00:00:00Z",
    "query":  {
    "ecmwf:resolution":{"eq": "high"},
    "ecmwf:levtype":{"eq": components_of_wind[0]["levtype"]},         
    "ecmwf:levelist": {"eq":  ["1000"]},      # 1000 hPa ~100 m (near sea level)],
    "ecmwf:time":{"eq": ["1200"]},
    "ecmwf:param":{"eq": [components_of_wind[0]["parameter_ID"],components_of_wind[1]["parameter_ID"]]}
    }
})

if(response.status_code!= 200):
    (print(response.text))
response.raise_for_status()

product = response.json()["features"][0]
JSON(product)

link = next((l for l in product.get('links', []) if l.get("rel") == "retrieve"), None)

if link:
    href = link.get("href")
    body = link.get("body")   # optional: depends on extension
    print("order endpoint:", href)
    print("order body, same as the polytope format:")
    print(json.dumps(body, indent=4))
else:
    print(f"No link with rel='{target_rel}' found")

#Sometimes requests to polytope get timeouts, it is then convenient define a retry strategy
retry_strategy = Retry(
    total=10,  # Total number of retries
    status_forcelist=[500, 502, 503, 504],  # List of 5xx status codes to retry on
    allowed_methods=["GET",'POST'],  # Methods to retry
    backoff_factor=1  # Wait time between retries (exponential backoff)
)

# Create an adapter with the retry strategy
adapter = HTTPAdapter(max_retries=retry_strategy)

# Create a session and mount the adapter
session = requests.Session()
session.mount("https://", adapter)


response = session.post(href, json=body, headers=auth_headers)

if response.status_code != 200:
    print(response.content)
response.raise_for_status()
    
ordered_item = response.json()
    
product_id = ordered_item["id"]
storage_tier = ordered_item["properties"].get("storage:tier", "online")
order_status = ordered_item["properties"].get("order:status", "unknown")
federation_backend = ordered_item["properties"].get("federation:backends", [None])[0]
    
print(f"Product ordered: {product_id}")
print(f"Provider: {federation_backend}")
print(f"Storage tier: {storage_tier} (product must have storage tier \"online\" to be downloadable)")
print(f"Order status: {order_status}")    

self_url = f"{HDA_STAC_ENDPOINT}/collections/{COLLECTION_ID}/items/{product_id}"
item = {}

#timeout and step for polling (sec)
TIMEOUT = 300
STEP = 1
ONLINE_STATUS = "online"

for i in range(0, TIMEOUT, STEP):
    print(f"Polling {i + 1}/{TIMEOUT // STEP}")
    
    response = session.get(self_url, headers=auth_headers)
    if response.status_code != 200:
        print(response.content)
    response.raise_for_status()
    item = response.json()
    
    storage_tier = item["properties"].get("storage:tier", ONLINE_STATUS)
    
    if storage_tier == ONLINE_STATUS:
        download_url = item["assets"]["downloadLink"]["href"]
        print("Product is ready to be downloaded.")
        print(f"Asset URL: {download_url}")
        break
    
    time.sleep(STEP)
else:
    order_status = item["properties"].get("order:status", "unknown")
    print(f"We could not download the product after {TIMEOUT // STEP} tries. Current order status is {order_status}")

response = session.get(download_url, stream=True, headers=auth_headers)
response.raise_for_status()
    
content_disposition = response.headers.get('Content-Disposition')
total_size = int(response.headers.get("content-length", 0))
if content_disposition:
    filename = content_disposition.split('filename=')[1].split('"')[1]
else:
    filename = os.path.basename(url)
        
# Open a local file in binary write mode and write the content
print(f"downloading {filename}")
    
with tqdm(total=total_size, unit="B", unit_scale=True) as progress_bar:
    with open(filename, 'wb') as f:
        for data in response.iter_content(1024):
            progress_bar.update(len(data))
            f.write(data)

import earthkit.data as ekd
ds = ekd.from_source("file", filename)
ds.ls()

import earthkit.plots
chart = earthkit.plots.Map(domain=[-10, 15, 35, 45])
import numpy as np
from earthkit.regrid import interpolate

# the target grid is a global 1x1 degree regular latitude-longitude grid
out_grid = {"grid": [1,1]}

# perform interpolation for each field and add results to a new fieldlist stored in memory
r = interpolate(ds, out_grid=out_grid, method="linear")
chart.land()
chart.coastlines()
chart.gridlines()
chart.title("Wind speed and direction ")

style = earthkit.plots.styles.Style(
    colors="plasma_r",
    levels=range(0, 22, 2),
    units="m s-1",
)

# Convert GRIB → xarray DataArray
u = r.sel(shortName="u").to_xarray()["u"]
v = r.sel(shortName="v").to_xarray()["v"]

# Compute wind speed
speed = np.hypot(u, v)
speed.name = "wind_speed"

# ---- coloured background ----
chart.quickplot(speed, style=style)

chart.quiver(r)
chart.legend(label="wind speed ({units})", location="right")
chart.show()

chart = earthkit.plots.Map(domain=[-10, 15, 35, 45])
chart.streamplot(r)

chart.land()
chart.coastlines()
chart.gridlines()

chart.title("{variable_name} at {level} hPa - {time:%-d %B %Y %H UTC}")

chart.show()