Set up the Python environment
If you have not done it yet, create a Python virtual 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 CDS 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
Plot Snow depth water equivalent map
Visualize GLOFAS products
Also for this exercise, we are going to retrieve data from the CDS API, as described in the code blocks below:
Retrieve data
The script shows how to retrieve the control reforecasts product from 1999 to 2018, relative to the date 2019-01-03, for two station coordinates, one on the river network of the Thames and the other one on the Po river.
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
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
YEARS = range(1999,2019)
# read dataset
ds = xr.open_dataset("glofas_reforecast_Po_january_03.grib",engine="cfgrib")
da = ds.isel(latitude=0,longitude=0).dis24
# plotting parameters
n = len(YEARS)
colors = plt.cm.jet(np.linspace(0,1,n))
ls = ["-","--",":"]
# plot
fig,ax = plt.subplots(1,1,figsize=(16,8))
steps = list(range(24,dsa.shape[1]*24+24,24))
for i,year in enumerate(YEARS):
y = da.sel(time=f"{year}-01-03").values
ax.plot(steps,y,label=f"{year}-01-03",color=colors[i],ls=ls[i%3])
plt.legend(bbox_to_anchor=(0.9,-0.1),ncol=7)
ax.set_xlim([0,steps[-1]])
ax.set_xlabel("leadtime hours")
ax.set_ylabel("m**3 s**-1")
plt.title(f"Po river discharge at {float(ds.latitude.values),round(float(ds.longitude.values),2)}")
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(2017,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.isel(latitude=0,longitude=0).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.isel(latitude=0,longitude=0).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(2017,7,19), datetime(2017,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(2017,7,19,3),85),rotation=0)
plt.xticks(rotation=70)
plt.title(f"Forecast reference time: 11-07-2007 - Thames river discharge at {float(ds.latitude.values),round(float(ds.longitude.values-360),2)}")
