page under construction ------------------------------
e.g. Example scripts for EFAS and GloFAS
Table of Content Zone | |
---|---|
|
...
Also for this exercise we are going to retrieve data from CDS. This time,
Retrieve data
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.
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 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 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
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:
...
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 an historic flood 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.
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
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.
...
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)}") |