Set up the Python environment
Activate the conda environment and install the additional Phython package
# activate the local environment. conda activate myenv # install the required packages conda install cartopy matplotlib # and the cdsapi pip install cdsapi
Visualize EFAS products
Retrieve data
These exercises require EFAS model output parameters and auxiliary data. You can retrieve them using the EWDS API, as described in the code blocks below:
Model outputs
Auxiliary data
One can download the auxiliary data by simplying requesting the data through a CDS API request.
Note that area sub-selection is not currently possible for the auxiliary (or time-invariant) data. Therefore any 'area' included in the CDS API request will be ignored and the data downloaded will cover the full EFAS domain.
And then unzip the auxiliary.zip
$ unzip auxiliary.zip Archive: auxiliary.zip extracting: sd_2_4.0.nc extracting: thmax_2_4.0.nc extracting: thmin_2_4.0.nc
Plot map discharge
Plot Soil Wetness Index Map
- The Soil Wetness Index is provided by EFAS starting from 30/05/2024, so we won't need to calculate it after that date. However, for data before that date, we will need to calculate the Soil Wetness Index using the following script. Please note that this process may require a significant amount of memory and time.
- Plot Soil Wetness Index provided by EFAS
Plot Snow depth water equivalent map
Visualize GloFAS products
Also for this exercise, we are going to retrieve data from the EWDS API, as described in the code blocks below:
Retrieve data
The script shows how to retrieve the control reforecasts product from 2003 to 2022, relative to the reference year 2023, for two station coordinates, one on the river network of the Thames and the other one on the Po river.
It's not possible to download only one cell of data , we have to provide and area .In the following script we are adding 0.05 degree to the station coordinates to get the surrounding area of the station . And for a correct interpretation of the data it's better to compare the uparea of the station with the uparea file provided .
This second script shows how to retrieve a point time series reforecast on the river Thames for a single forecast reference time, specifically the 11th of July 2007.
Plot ~20 years of reforecast time series from point coordinates
Plot reforecast ensemble members' time series for a historic flood event
In July 2007, a series of flooding events hit the UK, in particular in some areas of the upper Thames catchment. Up to 120 mm of rain fell between the 19th and 20th of July.
import xarray as xr
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from datetime import datetime,timedelta
forecast_reference_time = datetime(2007,7,11)
# read perturbed ensembles (dataType:pf)
ds = xr.open_dataset("glofas_reforecast_Thames_july_11.grib",engine="cfgrib",backend_kwargs={'filter_by_keys':{'dataType': 'pf'}})
da = ds.sel(latitude=51.3969,longitude=-0.4017, method='nearest').dis24
# read control (dataType:cf)
cds = xr.open_dataset("glofas_reforecast_Thames_july_11.grib",engine="cfgrib",backend_kwargs={'filter_by_keys':{'dataType': 'cf'}})
cda = cds.sel(latitude=51.3969,longitude=-0.4017, method='nearest').dis24
start = forecast_reference_time
end = start + timedelta(45)
time = pd.date_range(start,end)
# plotting parameters
ls = ["-","--",":"]
locator = mdates.AutoDateLocator(minticks=2, maxticks=46)
formatter = mdates.DateFormatter("%b-%d")
# plot
fig,ax = plt.subplots(1,1,figsize=(16,8))
ax.plot(time,cda.values,label=f"control",color="red")
for i,number in enumerate(range(0,10)):
y = da.isel(number=number).values
ax.plot(time,y,label=f"ensemble {number}",color=colors[i],ls=ls[i%3])
plt.legend(bbox_to_anchor=(1,1),ncol=2)
ax.axvspan(datetime(2007,7,19), datetime(2007,7,21), alpha=0.3, color='blue') # highlight event
ax.set_xlabel("date")
ax.set_ylabel("m**3 s**-1")
ax.xaxis.set_major_locator(locator)
ax.xaxis.set_major_formatter(formatter)
ax.set_xlim([time[0],time[-1]])
plt.annotate("peak\nrainfall\nevent",(datetime(2007,7,19,3),85),rotation=0)
plt.xticks(rotation=70)
plt.title(f"Forecast reference time: 11-07-2007 - Thames river discharge at {float(da.latitude.values),round(float(da.longitude.values-360),2)}")
