...
| Code Block | ||||||
|---|---|---|---|---|---|---|
|
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 |
|---|
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'], |
...
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:
| 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)}") |
...
'area':station_point_coord,
'hyear': YEAR,
'hmonth': month ,
'hday': day ,
'leadtime_hour': LEADTIMES,
},
f'glofas_reforecast_{station}_{month}_{day}.grib')
|