...
Code Block |
---|
language | py |
---|
title | Extract timeseries: GloFAS Medium-range Reforecast |
---|
collapse | true |
---|
|
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') |
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)}") |
Image Removed
Image Removed
...