page under construction ------------------------------

e.g. Example scripts for EFAS and GloFAS


EFAS


Retrieve data from CDS API

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')


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')


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')


Retrieve static data

This is available only when marsurl adaptor will be in production. For the moment I am using the datasets from the test stack.



Plot map discharge


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

ds = xr.open_dataset('../data/clim_19910300.nc')


cmap = plt.cm.get_cmap('jet').copy()
cmap.set_under('white')

crs = ccrs.LambertAzimuthalEqualArea(central_longitude=10,central_latitude=52,false_easting=4321000,false_northing=3210000)

# 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["dis06"].plot(ax=ax,cmap=cmap,vmin=20,add_colorbar=False)
cbar = plt.colorbar(sc, shrink=.5,)
cbar.set_label(ds.dis06.GRIB_name)



Timeseries

import numpy as np
import xarray as xr



Soil moisture

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

# helper 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


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')
ds = xr.open_dataset('eud_2019013000.nc')


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

sd1 = ds.sod[0,0,:,:].values
sd2 = ds.sod[0,1,:,:].values - ds.sod[0,0,:,:].values

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

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

data=((th1 * sd1) + (th2 * sd2))/(sd1 +sd2)
smtot = xr.DataArray(data,dims=(ds.y.name,ds.x.name),name='smtot')
ds['smtot'] = smtot 

data=(((smtot - thmin)/(thmax - thmin)))
SM = xr.DataArray(data,dims=(ds.y.name,ds.x.name),name='SM')
ds['SM'] = SM




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)

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("SM")



Snow depth water equivalent




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. 

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')



Plot retrieved data:


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.

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')


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:

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)}")