page under construction ------------------------------
...
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 |
...
EFAS
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 | |||
---|---|---|---|
|
...
| |||
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') |
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:
...
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
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 = 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 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,)
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
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
In order to plot a map, we are going to use
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
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") |
Plot Snow depth water equivalent map
...
language | py |
---|---|
title | Snow Depth Water Equivalent |
collapse | true |
GLOFAS
...
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
|
...
Plot retrieved data:
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)}") |