# Climate Adaptation Digital Twin: average temperature in Germany and Heating Degree Days in Darmstadt This notebook will provide you guidance on how to access and use Climate Adaptation Digital Twin (Climate DT) [datasets](https://earthdatahub.destine.eu/collections/climate-dt-2/datasets/IFS-NEMO-SSP3-7.0-sfc-hourly-standard) on [Earth Data Hub](https://earthdatahub.destine.eu) (EDH). Access to these datasets is granted via [Destination Earth](https://platform.destine.eu/)’s Standard API keys, but restricted to users with [Upgraded Access](https://platform.destine.eu/access-policy-upgrade). In this tutorial we will learn how to use Xarray to access the [Climate DT standard resolution hourly data on single-level or surface, IFS-NEMO, SSP3-7.0](https://earthdatahub.destine.eu/collections/climate-dt-2/datasets/IFS-NEMO-SSP3-7.0-sfc-hourly-standard) dataset on EDH. However, you can easily replicate the procedure for any other dataset belonging to the same collection. If you are running this tutorial in your own environment, you can download the convenience modules costing.py and display.py from [DESP-UserWorkflowService-Templates/EarthDataHub](https://github.com/SercoSPA/DESP-UserWorkflowService-Templates/tree/main/EarthDataHub) ## Goal of this tutorial The first goal of this tutorial is to plot the average 2 metre temperature in Germany for each year between 2015 and 2049. The second goal of this tutorial is to calculate the Heating Degree Days in Darmstadt in the same years. ## What you will learn: - how to access the dataset - select and reduce the data - plot the results --- ## Working with EDH data Datasets on EDH are typically very large and remotely hosted. Typical use imply a selection of the data followed by one or more reduction steps to be performed in a local or distributed Dask environment. The structure of a workflow that uses EDH data typically looks like this: 1. Data access 2. Data selection 3. Data reduction (optional) 4. Data download 5. Further operations and visualisation [Xarray](https://docs.xarray.dev/en/stable) and [Dask](https://www.dask.org) work together following a lazy principle. This means that when you access and manipulate a Zarr store the data is in not immediately downloaded and loaded in memory. Instead, Dask constructs a task graph that represents the operations to be performed. A smart user will first reduce the amount of data that needs to be downloaded and explicitly call `compute()` on it. Once the `compute()` operation is complete the data is loaded into memory and available for subsequent fast processing. ### 1. Data access ```python # Optional: ignore warnings import warnings warnings.filterwarnings('ignore') ``` To access any data on EDH, you need to obtain an API key. To obtain an API key you need to register to the [DestinE platform](https://platform.destine.eu). Then, you can go to [Earth Data Hub Account Settings - Quota & API keys](https://earthdatahub.destine.eu/account-settings) where you can find out your default API key or create others. Finally, you need to instruct your tools to use your API key when downloading the data. One way to do it, is to incorporate the API key into the dataset’s url. If you wish to follow this method, cut and paste your personal API key below: ⤵ ```python EDH_KEY = "your_EDH_API_key" #e.g. EDH_KEY="edh_key_MHOQO_Y_h4ywpaMQe4cndtydvbpU6ySBNARWprm5tgz" url = f"https://edh:{EDH_KEY}@api.earthdatahub.destine.eu/climate-dt-2/IFS-NEMO-SSP3-7.0-sfc-hourly-standard-v0.zarr" ``` There are more ways in which you can provide your API key to Xarray, especially if you plan to use the system frequently. Please refer to [https://earthdatahub.destine.eu/getting-started](https://earthdatahub.destine.eu/getting-started) for a complete overview. We are now ready to access the data with Xarray: ```python import xarray as xr xr.set_options(keep_attrs=True) ds = xr.open_dataset( url, chunks={}, engine="zarr", zarr_format=3, ) ds ``` ```text Size: 22TB Dimensions: (time: 306816, latitude: 512, longitude: 1025) Coordinates: * time (time) datetime64[ns] 2MB 2015-01-01 ... 2049-12-31T23:00:00 * latitude (latitude) float64 4kB -90.0 -89.65 -89.3 ... 89.3 89.65 90.0 * longitude (longitude) float64 8kB -180.0 -179.6 -179.3 ... 179.6 180.0 Data variables: (12/36) avg_ie (time, latitude, longitude) float32 644GB dask.array avg_iews (time, latitude, longitude) float32 644GB dask.array avg_inss (time, latitude, longitude) float32 644GB dask.array avg_ishf (time, latitude, longitude) float32 644GB dask.array avg_sdlwrf (time, latitude, longitude) float32 644GB dask.array avg_sdswrf (time, latitude, longitude) float32 644GB dask.array ... ... tciw (time, latitude, longitude) float32 644GB dask.array tclw (time, latitude, longitude) float32 644GB dask.array tcw (time, latitude, longitude) float32 644GB dask.array tcwv (time, latitude, longitude) float32 644GB dask.array u10 (time, latitude, longitude) float32 644GB dask.array v10 (time, latitude, longitude) float32 644GB dask.array Attributes: GRIB_edition: 2 GRIB_centre: ecmf GRIB_centreDescription: European Centre for Medium-Range Weather Forecasts GRIB_subCentre: 1003 Conventions: CF-1.7 institution: European Centre for Medium-Range Weather Forecasts ``` ⚠ At this point, no data has been downloaded yet, nor loaded in memory. ### 2. Data selection We first select the 2 metres temperature. ```python t2m = ds.t2m t2m ``` ```text Size: 644GB dask.array Coordinates: * time (time) datetime64[ns] 2MB 2015-01-01 ... 2049-12-31T23:00:00 * latitude (latitude) float64 4kB -90.0 -89.65 -89.3 ... 89.3 89.65 90.0 * longitude (longitude) float64 8kB -180.0 -179.6 -179.3 ... 179.6 180.0 Attributes: (12/20) GRIB_paramId: 167 GRIB_dataType: fc GRIB_numberOfPoints: 196608 GRIB_typeOfLevel: heightAboveGround GRIB_stepUnits: 1 GRIB_stepType: instant ... ... GRIB_shortName: 2t GRIB_units: K long_name: 2 metre temperature units: K standard_name: air_temperature last_ingested_time_index: 306815 ``` As xarray showes, the dimension of the `t2m` DataArray is around 600GiB (in memory). We will now try to narrow down the selection as much as possible. We select only the Germany area: ```python germany = {'latitude': slice(46, 56), 'longitude': slice(5, 16)} t2m_germany = t2m.sel(**germany) t2m_germany ``` ```text Size: 1GB dask.array Coordinates: * time (time) datetime64[ns] 2MB 2015-01-01 ... 2049-12-31T23:00:00 * latitude (latitude) float64 224B 46.32 46.67 47.03 ... 55.13 55.48 55.83 * longitude (longitude) float64 248B 5.273 5.625 5.977 ... 15.12 15.47 15.82 Attributes: (12/20) GRIB_paramId: 167 GRIB_dataType: fc GRIB_numberOfPoints: 196608 GRIB_typeOfLevel: heightAboveGround GRIB_stepUnits: 1 GRIB_stepType: instant ... ... GRIB_shortName: 2t GRIB_units: K long_name: 2 metre temperature units: K standard_name: air_temperature last_ingested_time_index: 306815 ``` Notice that the size of the array dropped down to around 1GiB. However, due to the chunked structure of the DataArray, xarray must download every chunk that includes a portion of the selected data. To estimate the size of the download, we can use the *costing.py* module. This estimate must be done before we apply any reduction operation. ```python import costing costing.estimate_download_size(t2m_germany) ``` ```text Estimated download size: 0.94 GiB (963.0 MiB) ``` ### 3. Data reduction We average the 2 metres temperature by year. This also is a lazy operation. ```python t2m_germany_yearly_mean = t2m_germany.resample(time='YE').mean(dim="time") t2m_germany_yearly_mean ``` ```text Size: 122kB dask.array Coordinates: * time (time) datetime64[ns] 280B 2015-12-31 2016-12-31 ... 2049-12-31 * latitude (latitude) float64 224B 46.32 46.67 47.03 ... 55.13 55.48 55.83 * longitude (longitude) float64 248B 5.273 5.625 5.977 ... 15.12 15.47 15.82 Attributes: (12/20) GRIB_paramId: 167 GRIB_dataType: fc GRIB_numberOfPoints: 196608 GRIB_typeOfLevel: heightAboveGround GRIB_stepUnits: 1 GRIB_stepType: instant ... ... GRIB_shortName: 2t GRIB_units: K long_name: 2 metre temperature units: K standard_name: air_temperature last_ingested_time_index: 306815 ``` ### 4. Data download This is the phase where we explicitly trigger the download of the data. Remember to assign the return of the `compute()` function to a new variable, so that the data is kept in memory. ```python %%time t2m_germany_yearly_mean_computed = t2m_germany_yearly_mean.compute() ``` ```text CPU times: user 25.8 s, sys: 2.16 s, total: 27.9 s Wall time: 1min 6s ``` ### 5. Visualisation We will now create and display an animation of the average winter 2 metres temperature in Germany, for the years 2020-2039. ```python t2m_germany_yearly_mean_computed_degC = t2m_germany_yearly_mean_computed - 273.15 t2m_germany_yearly_mean_computed_degC.attrs["units"] = "°C" ``` ```python import pandas as pd import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature from cartopy import crs from matplotlib.animation import FuncAnimation from IPython.display import HTML from mpl_toolkits.axes_grid1 import make_axes_locatable fig, ax = plt.subplots(subplot_kw={'projection': ccrs.PlateCarree()}) ax.add_feature(cfeature.COASTLINE, linewidth=0.5) ax.add_feature(cfeature.BORDERS, linestyle=':', linewidth=0.5) ax.add_feature(cfeature.OCEAN, facecolor='lightblue', zorder=2) ax.gridlines(draw_labels=True, zorder=3, color="white", alpha=0.5) t2m_germany_yearly_mean_computed_degC.isel(time=0).plot(ax=ax, transform=ccrs.PlateCarree(), cmap='RdBu_r', add_colorbar=True, cbar_kwargs={'orientation': 'vertical', 'shrink':0.9,'pad': 0.15}) # function to update the plot for each frame (each timestep) def update(frame): data = t2m_germany_yearly_mean_computed_degC.isel(time=frame) plot = data.plot( ax=ax, transform=ccrs.PlateCarree(), cmap='RdBu_r', vmin=-25, vmax=25, add_colorbar=False ) ax.set_title(f"Time: {pd.Timestamp(data['time'].values).strftime('%Y-%m-%d')}") return plot anim = FuncAnimation(fig, update, frames=len(t2m_germany_yearly_mean_computed_degC['time']), repeat=True) # Create the animation plt.close() # close the static plot to avoid duplicate display HTML(anim.to_jshtml()) # display the animation in the notebook ``` ```text ``` ## Heating Degree Days (HDD) in Darmstadt Let us now investigate the Heating Degree Days in Darmstadt. We start again from the 2 metres temperature. ```python t2m ``` ```text Size: 644GB dask.array Coordinates: * time (time) datetime64[ns] 2MB 2015-01-01 ... 2049-12-31T23:00:00 * latitude (latitude) float64 4kB -90.0 -89.65 -89.3 ... 89.3 89.65 90.0 * longitude (longitude) float64 8kB -180.0 -179.6 -179.3 ... 179.6 180.0 Attributes: (12/20) GRIB_paramId: 167 GRIB_dataType: fc GRIB_numberOfPoints: 196608 GRIB_typeOfLevel: heightAboveGround GRIB_stepUnits: 1 GRIB_stepType: instant ... ... GRIB_shortName: 2t GRIB_units: K long_name: 2 metre temperature units: K standard_name: air_temperature last_ingested_time_index: 306815 ``` We narrow down the selection to the data that is closer to Darmstadt. ```python darmstadt = {"latitude": 49.88, "longitude": 8.65} base_temperature = 288.15 #[K] t2m_darmstadt = t2m.sel(darmstadt, method="nearest") t2m_darmstadt ``` ```text Size: 1MB dask.array Coordinates: * time (time) datetime64[ns] 2MB 2015-01-01 ... 2049-12-31T23:00:00 latitude float64 8B 49.84 longitude float64 8B 8.789 Attributes: (12/20) GRIB_paramId: 167 GRIB_dataType: fc GRIB_numberOfPoints: 196608 GRIB_typeOfLevel: heightAboveGround GRIB_stepUnits: 1 GRIB_stepType: instant ... ... GRIB_shortName: 2t GRIB_units: K long_name: 2 metre temperature units: K standard_name: air_temperature last_ingested_time_index: 306815 ``` We estimate the cost of the download with the `costing.py` module. ```python costing.estimate_download_size(t2m_darmstadt) ``` ```text Estimated download size: 0.94 GiB (963.0 MiB) ``` We compute the HDD in Darmstadt with a very simplified formula: ```python t2m_darmstadt_daily_mean = t2m_darmstadt.resample(time='1D').mean(dim='time') diff = (base_temperature - t2m_darmstadt_daily_mean) hdd = diff.where(diff > 0).groupby("time.year").sum() ``` We explicitly trigger the download of the data. Remember to assign the return of the `compute()` function to a new variable, so that the data is kept in memory. ```python %%time hdd_computed = hdd.compute() ``` ```text CPU times: user 15.3 s, sys: 1.5 s, total: 16.8 s Wall time: 1min 5s ``` We can finally visualize the HDD in Darmstadt. ```python hdd ``` ```text Size: 140B dask.array Coordinates: * year (year) int64 280B 2015 2016 2017 2018 ... 2046 2047 2048 2049 latitude float64 8B 49.84 longitude float64 8B 8.789 Attributes: (12/20) GRIB_paramId: 167 GRIB_dataType: fc GRIB_numberOfPoints: 196608 GRIB_typeOfLevel: heightAboveGround GRIB_stepUnits: 1 GRIB_stepType: instant ... ... GRIB_shortName: 2t GRIB_units: K long_name: 2 metre temperature units: K standard_name: air_temperature last_ingested_time_index: 306815 ``` ```python plt.style.use("seaborn-v0_8-darkgrid") fig, ax = plt.subplots(figsize=(12, 5)) plt.bar(hdd_computed.year, hdd.values, color='#ff0000', alpha=0.7) plt.xlabel('time') plt.ylabel('HDD [°C]') plt.grid(axis='y', alpha=0.75) plt.title('Heating Degrees Days in Darmstadt') plt.xticks(hdd.year[::2], rotation=45); ``` ![image](https://earthdatahub.destine.eu/docs/_images/notebooks_insula-code-lab_tutorial-climate-dt-heating-degree-days-in-darmstadt_39_0.png)