page under construction ------------------------------
e.g. Example scripts for EFAS and GloFAS
Table of Content Zone | |
---|---|
|
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
Code Block | |||||
---|---|---|---|---|---|
|
...
|
...
# |
...
activate |
...
the local |
...
environment |
...
. conda |
...
activate myenv |
...
# |
...
install the |
...
required |
...
packages conda |
...
install |
...
cartopy matplotlib # |
...
EFAS
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() Â c.retrieve( 'efas-forecast', { 'format': 'netcdf', 'originating_centre': 'ecmwf', 'product_type': 'ensemble_perturbed_forecasts', 'variable': 'river_discharge_in_the_last_6_hours', 'model_levels': 'surface_level', 'year': '2020', 'month': '11', 'day': '15', 'time': '00:00', 'leadtime_hour': [ '6', '12', '18', '24', '30', '36', '42', '48', '54', '60', '66', '72', ], }, 'eue_2020111500.nc') |
...
Code Block | ||||||
---|---|---|---|---|---|---|
| ||||||
import cdsapi c = cdsapi.Client() Â Â c.retrieve( 'efas-forecast', { 'format': 'netcdf', 'originating_centre': 'ecmwf', 'product_type': 'high_resolution_forecast', 'variable': [ 'soil_depth |
...
' |
...
],
'model_levels': 'soil_levels',
'year': '2019',
'month': '01',
'day': '30',
'time': '00:00',
'leadtime_hour': [
'6', '12', '18',
'24', '30', '36',
'42', '48', '54',
'60', '66', '72',
],
'soil_level': [
'1', '2', '3',
],
},
'eud_2019013000.nc') |
Code Block |
---|
Auxiliary data
From the CDS Static files link, download the following NetCDFs:
thmin1.nc, thmin2.nc, thmin3.nc, thmax1.nc, thmax2.nc, thmax3.nc
When the marsurl will be in production one can download the auxiliary data simply requesting them through a CDS API request:
...
|
...
| |||
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
Info | |||||||||||||||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
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
|
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('./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 = |
And then unzip the auxiliary.zip
Code Block | ||||||
---|---|---|---|---|---|---|
| ||||||
$ unzip auxiliary.zip
Archive: auxiliary.zip
extracting: thmax_1.nc
extracting: thmin_1.nc
extracting: thmax_2.nc
extracting: thmin_2.nc
extracting: thmax_3.nc
extracting: thmin_3.nc
|
Plot map discharge
...
language | py |
---|---|
title | Discharge map |
collapse | true |
...
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
...
language | py |
---|---|
title | Discharge timeseries forecasts |
collapse | true |
...
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('./thmin_1_4.0.nc') thmin2 = xr.open_dataset('./thmin_2_4.0.nc') thmin3 = xr.open_dataset('./thmin_3_4.0.nc') thmax1 = xr.open_dataset(' |
...
./thmax_1_4.0.nc') |
...
thmax2 = xr.open_dataset(' |
...
./thmax_2_4.0.nc') |
...
thmax3 = xr.open_dataset(' |
...
./thmax_3_4.0.nc') |
...
sd1 = xr.open_dataset(' |
...
./sd_1_4.0.nc') |
...
sd2 = xr.open_dataset( |
...
'./sd_2_4.0.nc') |
...
sd3 = xr.open_dataset(' |
...
./sd_3_4.0.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 |
...
Plot Snow depth water equivalent map
...
language | py |
---|---|
title | Snow Depth Water Equivalent |
collapse | true |
GLOFAS
...
,)
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
...
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
# get date index
MONTHSDAYS = get_monthsdays()
if __name__ == '__main__':
c = cdsapi.Client()
# set station coordinates (lat,lon)
COORDS = {
"Thames":[51.35,-0.45],
"Po":[44.85, 11.65]
}
# select all years
YEARS = ['%d'%(y) for y in range(1999,2019)]
# select all leadtime hours
LEADTIMES = ['%d'%(l) for l in range(24,1128,24)]
# loop over date index
for md in MONTHSDAYS:
month = md[0].lower()
day = md[1]
# loop over 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',
'area':station_point_coord,
'hyear': YEARS,
'hmonth': month ,
'hday': day ,
'leadtime_hour': LEADTIMES,
},
f'glofas_reforecast_{station}_{month}_{day}.grib') |
...
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.
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)}")
GloFAS Medium-range reforecast (download time series for an event)
This 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.
...
language | sass |
---|
...
| |||||
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.
Plot retrieved data:
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)}") |