Image Added

Table of Content Zone

Table of Contents

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

language	bash
theme	Emacs
title	install required packages

# 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

Code Block

language	py
title	Historical river discharge
collapse	true

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

language	py
title	Forecast river discharge
collapse	true

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

language	py
title	Forecast volumetric soil moisture and soil depth
collapse	true

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

language	py
title	Forecast snow depth water equivalent
collapse	true

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

...

title	MARS url adaptor

...

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.

Code Block

title
language	py

...

Retrieving a single auxiliary

...

file
collapse	true

import cdsapi

c = cdsapi.Client()

c.retrieve(
    'efas-forecast',
    {
        '

...

variable': '

...

upstream_area_v4_0',
        '

...

model_levels':

...

'surface_level',

...

format'

...

: '

...

netcdf',

...

area': [

...

72, -35, 24,
            74,
        ],
    },
    'upstream.nc')

Code Block

language	py
title	Retrieve multiple auxiliary files
collapse	true

import cdsapi

c = cdsapi.Client()

c.retrieve(
    'efas-forecast',
    {
        'variable': ['soil_depth_v4_0','field_capacity_v4_0','wilting_point_v4_0'],
        'model_levels': 'soil_levels',
        'soil_level': '2',
        'format': 'netcdf',
        'area': [
            72, -35, 24,
            74,
        ],
    },
    'auxiliary.zip')

And then unzip the auxiliary.zip

Code Block

language	bash
theme	Emacs
title	unzip

$ 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

Code Block

language	py
title	Discharge map
collapse	true

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 = 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,

...

Image Removed

Plot Soil Wetness Index Map

...

title	Note for when marsurl is ready in production

In the code block below the soild depth will need to be loaded from the netCDF file:

xr.open_dataset('./sd_1.nc')

xr.open_dataset('./sd_2.nc')

xr.open_dataset('./sd_3.nc')

)
cbar.set_label(ds[var_name].GRIB_name)

Image Added

Plot Soil Wetness Index Map

...

Code Block

language	py
title	Soil Wetness Index Map
collapse	true

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,)
cbar.set_label("Soil Wetness Index")

Image Modified

Plot Snow depth water equivalent map

Code Block

language	py
title	Snow Depth Water Equivalent
collapse	true

import xarray as xr
ds = xr.open_dataset('./esd_2021013000.nc')

ds.sd.plot(col="step",col_wrap=4, robust=True)

 Image Modified

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

language	py
title	Retrieve ~20 years of reforecasts for two locations
collapse	true

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

language	py
title	Retrieve reforecast flood event on the Thames in July 2007
collapse	true

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

language	py

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

Image ModifiedImage Modified

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.

Image Modified

Code Block

language	py
title	Plot

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

Image Modified

Page tree

Versions Compared

Old Version 42

New Version Current

Key

Set up the Python environment

Visualize EFAS products

Retrieve data

Model outputs

Auxiliary data

Plot map discharge

Plot Soil Wetness Index Map

Plot Soil Wetness Index Map

Code Block
language py
title Soil Wetness Index Map
collapse true
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

Plot Snow depth water equivalent map

Visualize GLOFAS products

Retrieve data

Plot ~20 years of reforecast time series from point coordinates

Plot reforecast ensemble members' time series for a historic flood event

Page tree

Page History

Versions Compared

Old Version 42

New Version Current

Key

Set up the Python environment

Visualize EFAS products

Retrieve data

Model outputs

Auxiliary data

Plot map discharge

Plot Soil Wetness Index Map

Plot Soil Wetness Index Map

Code BlocklanguagepytitleSoil Wetness Index Mapcollapsetrueimport 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

Plot Snow depth water equivalent map

Visualize GLOFAS products

Retrieve data

Plot ~20 years of reforecast time series from point coordinates

Plot reforecast ensemble members' time series for a historic flood event

Code Block
language py
title Soil Wetness Index Map
collapse true
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