...
| Code Block | ||||||
|---|---|---|---|---|---|---|
| ||||||
# create a local virtual environment, you can call it as you wish, here 'myenv' is used. conda create -n myenv python=3.8 # add repository channel conda config --add channels conda-forge # activate the local environment. conda activate myenv # install the required packages conda install -c conda-forge/label/main xarray cfgrib eccodes conda install cartopy netcdf4 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
| Code Block | ||||||
|---|---|---|---|---|---|---|
| ||||||
import cdsapi
c = cdsapi.Client()
c.retrieve('efas-historical', {
"format": "netcdf",
"hday": ["15","16","17","18"],
"hmonth": "november",
"hyear": "2020",
"model_levels": "surface_level",
"system_version": "version_4_0",
"time": ["00:00","06:00","18:00"],
"variable": "river_discharge_in_the_last_6_hours"
},
'clim_2020111500.nc') |
...
| Code Block | ||||||
|---|---|---|---|---|---|---|
| ||||||
import cdsapi
c = cdsapi.Client()
# forecast snow depth water equivalent
c.retrieve(
'efas-forecast',
{
'format': 'netcdf',
'originating_centre': 'ecmwf',
'variable': 'snow_depth_water_equivalent',
'product_type': 'control_forecast',
'model_levels': 'surface_level',
'year': '2021',
'month': '01',
'day': '30',
'time': '00:00',
'leadtime_hour': [
'6', '12', '18',
'24', '30', '36',
'42', '48', '54',
'60', '66', '72',
'78', '84', '90',
'96', '102',
],
},
'esd_2021013000.nc') |
Auxiliary data
From the CDS Static files link, download the following NetCDF:
thmin1.nc, thmin2.nc, thmin3.nc, thmax1.nc, thmax2.nc, thmax3.nc
| Info | ||||||||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
And then unzip the auxiliary.zip
|
Plot map discharge
| Code Block | ||||||
|---|---|---|---|---|---|---|
| ||||||
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cf
import numpy as np
import pandas as pd
import xarray as xr
var_name = "dis06"
ds = xr.open_dataset('../data/clim_2020111500.nc')
cmap = plt.cm.get_cmap('jet').copy()
cmap.set_under('white')
# set the coordinate reference system for EFAS
crs = ccrs.LambertAzimuthalEqualArea(central_longitude=10,central_latitude=52,false_easting=4321000,false_northing=3210000)
# define the filter for visualizing only discharge above that value
vmin = 20
# selecting a date
ds = ds[var_name].isel(time=1)
# Plot map discharge > 20 m/s
fig, ax = plt.subplots(1,1,subplot_kw={'projection': crs}, figsize=(20,20) )
ax.gridlines(crs=crs, linestyle="-")
ax.coastlines()
ax.add_feature(cf.BORDERS)
sc = ds[var_name].plot(ax=ax,cmap=cmap,vmin=vmin,add_colorbar=False)
ax.set_title(f'{ds[var_name].long_name}> {vmin} $m^{3}s^{-1}$')
cbar = plt.colorbar(sc, shrink=.5,)
cbar.set_label(ds[var_name].GRIB_name) |
Plot Discharge Timeseries
| Code Block | ||||||
|---|---|---|---|---|---|---|
| ||||||
import numpy as np import xarray as xr |
Plot Soil Wetness Index Map
| Code Block | ||||||
|---|---|---|---|---|---|---|
| ||||||
import xarray as xr
from matplotlib.pyplot import colorbar
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cf
import numpy as np
import pandas as pd
# Load the wilting point and field capacity datasets
thmin1 = xr.open_dataset('thmin1.nc')
thmin2 = xr.open_dataset('thmin2.nc')
thmin3 = xr.open_dataset('thmin3.nc')
thmax1 = xr.open_dataset('thmax1.nc')
thmax2 = xr.open_dataset('thmax2.nc')
thmax3 = xr.open_dataset('thmax3.nc')
# load the EFAS volumetric soil moisture dataset
ds = xr.open_dataset('eud_2019013000.nc')
# plotting function
def make_cmap(colors, position=None):
'''
make_cmap takes a list of tuples which contain RGB values.
The RGB values are [0 to 255]
make_cmap returns a cmap with equally spaced colors.
Arrange your tuples so that the first color is the lowest value for the
colorbar and the last is the highest.
position contains values from 0 to 1 to dictate the location of each color.
'''
import matplotlib as mpl
import numpy as np
import sys
bit_rgb = np.linspace(0,1,256)
if position == None:
position = np.linspace(0,1,len(colors))
else:
if len(position) != len(colors):
sys.exit("position length must be the same as colors")
elif position[0] != 0 or position[-1] != 1:
sys.exit("position must start with 0 and end with 1")
for i in range(len(colors)):
colors[i] = (bit_rgb[colors[i][0]],
bit_rgb[colors[i][1]],
bit_rgb[colors[i][2]])
cdict = {'red':[], 'green':[], 'blue':[]}
for pos, color in zip(position, colors):
cdict['red'].append((pos, color[0], color[0]))
cdict['green'].append((pos, color[1], color[1]))
cdict['blue'].append((pos, color[2], color[2]))
cmap = mpl.colors.LinearSegmentedColormap('my_colormap',cdict,256)
return cmap
# There are 3 layers of Volumetric Soil Moisture. We need to get the mean across all timesteps
mean_vsw = ds.vsw.mean('step') # Mean Volumetric Soil Moisture for all time steps
th1 = mean_vsw[0,:,:] # Mean Volumetric Soil Moisture for Layer 1
th2 = mean_vsw[1,:,:] # Mean Volumetric Soil Moisture for Layer 2
# We need to have the values for the first and second soil depths
# Layer 2 includes Layer 1, so we need to remove 2 from 1
# All layers are the same through out so we need to use step 1
sd1 = ds.sod[0,0,:,:].values
sd2 = ds.sod[0,1,:,:].values - ds.sod[0,0,:,:].values
# Compute the mean wilting point
data=(thmin1.wiltingpoint.values * sd1 + thmin2.wiltingpoint.values * sd2)/(sd1 + sd2)
thmin = xr.DataArray(data,dims=(ds.y.name,ds.x.name),name='thmin')
ds['thmin'] = thmin
# Compute the field capacity
data=(thmax2.fieldcapacity.values * sd1 + thmax2.fieldcapacity.values * sd2)/(sd1 + sd2)
thmax = xr.DataArray(data,dims=(ds.y.name,ds.x.name),name='thmax')
ds['thmax'] = thmax
# Compute the mean soil moisture for layer 1 and 2
data=((th1 * sd1) + (th2 * sd2))/(sd1 +sd2)
smtot = xr.DataArray(data,dims=(ds.y.name,ds.x.name),name='smtot')
ds['smtot'] = smtot
# Compute the soil wetness index
data=(((smtot - thmin)/(thmax - thmin)))
SM = xr.DataArray(data,dims=(ds.y.name,ds.x.name),name='SM')
ds['SM'] = SM
# Create a list of RGB colors for the plotting function
colors = [(64,0,3), (133,8,3), (249,0,0), (250,106,0), (249,210,0), (189,248,255), (92,138,255), (0,98,255), (0,0,255), (0,0,88)]
# Create an array or list of positions from 0 to 1.
position = [0, 0.2, 0.3,0.4,0.5,0.6,0.7,0.8,0.9,1]
cmap = make_cmap(colors, position=position)
# define the coordinate reference system for the EFAS data
crs = ccrs.LambertAzimuthalEqualArea(central_longitude=10,central_latitude=52,false_easting=4321000,false_northing=3210000)
fig, ax = plt.subplots(1,1,subplot_kw={'projection': crs}, figsize=(15,15) )
ax.gridlines(crs=crs, linestyle="-")
ax.coastlines()
ax.add_feature(cf.BORDERS)
sc = ds["SM"].plot(ax=ax,cmap=cmap,add_colorbar=False)
cbar = plt.colorbar(sc, shrink=.5,)
cbar.set_label("Soil Wetness Index") |
Plot Snow depth water equivalent map
| Code Block | ||||||
|---|---|---|---|---|---|---|
| ||||||
import xarray as xr
ds = xr.open_dataset('esd_2021013000.nc')
ds.sd.plot(col="step",col_wrap=4, robust=True) |
...
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 year 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.
...
| Code Block | ||||||
|---|---|---|---|---|---|---|
| ||||||
import cdsapi
from datetime import datetime, timedelta
def get_monthsdays(start =[2019,1,1],end=[2019,12,31]):
# reforecast time index
start, end = datetime(*start),datetime(*end)
days = [start + timedelta(days=i) for i in range((end - start).days + 1)]
monthday = [d.strftime("%B-%d").split("-") for d in days if d.weekday() in [0,3] ]
return monthday
if __name__ == '__main__':
c = cdsapi.Client()
# station coordinates (lat,lon)
COORDS = {
"Thames":[51.35,-0.45]
}
# select date index corresponding to the event
MONTHSDAYS = get_monthsdays(start =[2019,7,11],end=[2019,7,11])
YEAR = '2007'
LEADTIMES = ['%d'%(l) for l in range(24,1128,24)]
# loop over date index (just 1 in this case)
for md in MONTHSDAYS:
month = md[0].lower()
day = md[1]
# loop over station coordinates
for station in COORDS:
station_point_coord = COORDS[station]*2 # coordinates input for the area keyword
c.retrieve(
'cems-glofas-reforecast',
{
'system_version': 'version_2_2',
'variable': 'river_discharge_in_the_last_24_hours',
'format': 'grib',
'hydrological_model': 'htessel_lisflood',
'product_type': ['control_reforecast','ensemble_perturbed_reforecasts'],
'area':station_point_coord,
'hyear': YEAR,
'hmonth': month ,
'hday': day ,
'leadtime_hour': LEADTIMES,
},
f'glofas_reforecast_{station}_{month}_{day}.grib')
|
Plot ~20 years of reforecast time series from point coordinates
| Code Block | ||
|---|---|---|
| ||
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.
| Code Block | ||||
|---|---|---|---|---|
| ||||
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)}") |
