Climate DT Parameter Series Plot- Data Access using DEDL HDA

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

import destinelab as deauth
import json

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)

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("The discovery result give us:\nthe collection ID : ", discovery_json["collections"][0].get("id"))
print("\nIts time extension : ", discovery_json["collections"][0].get("extent").get("temporal").get("interval"))
#print("\nThe available parameters: ", discovery_json["collections"][0].get("cube:variables").keys())
parameters=discovery_json["collections"][0]["cube:variables"]
print("\nThe available parameters: ")
keys = sorted(parameters)
print(json.dumps(keys, indent=2))
COLLECTION_ID = discovery_json["collections"][0].get("id")
TARGET_COLLECTION_ID ="EO.ECMWF.DAT.D1.DT_CLIMATE.G1.SCENARIOMIP_SSP3-7.0_IFS-NEMO.R1"

if COLLECTION_ID != TARGET_COLLECTION_ID:
    print("Target collection not found in discovery response.")

for var_name, var_info in parameters.items():
    if var_name=="2_metre_temperature(sfc)":
        var_type = var_info.get("type")
        var_unit = var_info.get("unit")
        var_url = var_info.get("attrs").get("url")

        two_metre_temperature = {
            "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")
        }
two_metre_temperature

response = requests.post(HDA_STAC_ENDPOINT+"/search", headers=auth_headers, json={
 "collections": [COLLECTION_ID],
    "datetime": "2020-01-01T00:00:00Z",
    "query":  {
    "ecmwf:resolution":{"eq": "high"},
    "ecmwf:levtype":{"eq": two_metre_temperature["levtype"]},
    "ecmwf:time":{"eq": ["1200"]},
    "ecmwf:param":{"eq": [two_metre_temperature["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=5,  # 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)


# Initialize a list to store filenames
filenames = []

# Define start and end years
start_year = 2020
end_year = 2039

dates=[]
# Loop 
for year in range(start_year, end_year+1,5):
    # Create a datetime object 
    obsdate = datetime(year, 7, 31)
    datechoice = obsdate.strftime("%Y%m%d")
    dates.append(datechoice)

dates

from concurrent.futures import ThreadPoolExecutor, as_completed
#timeout and step for polling (sec)
TIMEOUT = 300
STEP = 1
ONLINE_STATUS = "online"

def download_date(date):

    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 = {}
    
    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)
        
    filenames.append(filename)
    # 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)

    result=filename+' NOT'
    if (os.path.exists(filename)):
        result=filename
    return result

# Run downloads in parallel
with ThreadPoolExecutor(max_workers=4) as executor:
    futures = [executor.submit(download_date, d) for d in dates]

    for future in as_completed(futures):
        result = future.result()
        print(result," downloaded")

import earthkit.data
import earthkit.plots

# Iterate over filenames
data=[]
for filename in filenames:
    data.append( earthkit.data.from_source("file", filename))

STYLE = earthkit.plots.styles.Style(
    colors="Spectral_r",
    levels=range(-20, 45),
    units="celsius",
    # Extend the colorbar at both ends
    extend="both",
)

figure = earthkit.plots.Figure(rows=2, columns=2, size=(10, 10))
for i, year in enumerate(range(start_year, end_year,5)):
    subplot = figure.add_map()
    subplot.contourf(data[i], style=STYLE)
    subplot.title(f"{year} - Climate DT - IFS NEMO {{short_name!u}} {{time:%d %B %Y}}")


figure.title("{variable_name} in {time:%B %d} of {time:%Y}", fontsize=14)

figure.legend(label="{variable_name!l} ({units})")

figure.show()

for filename in filenames:
    os.remove(filename)