page under construction ------------------------------
e.g. Example scripts for EFAS and GloFAS
Table of Content Zone | |
---|---|
|
Set up the environment
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 |
EFAS
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() # Climatology river discharge 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') # Forecast river discharge 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') # Forecast soil moisture c.retrieve( 'efas-forecast', { 'format': 'netcdf', 'originating_centre': 'ecmwf', 'product_type': 'high_resolution_forecast', 'variable': [ 'soil_depth', 'volumetric_soil_moisture', ], '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') # 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 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:
Note: the soil depth is now part of the auxiliary data. Before it was downloaded as you would download a model output.
Code Block | ||||||
---|---|---|---|---|---|---|
| ||||||
import cdsapi c = cdsapi.Client() c.retrieve( 'efas-forecast, { 'format': 'netcdf', 'variable': [ 'field_capacity', 'wilting_point', 'soil_depth' ], 'soil_level': [ '1', '2', '3', ], }, 'auxiliary.zip') |
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
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) |
GLOFAS
GloFAS Medium-range reforecast (download time series)
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.
| Plot retrieved data:
|
---|
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.
| 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 19th and 20th of July. Plot retrieved data:
|
---|