Versions Compared

Key

  • This line was added.
  • This line was removed.
  • Formatting was changed.

...

  • Locally (full control on the process)


GloFAS

...

Example script to crop and extract time series from different GloFAS products:

Code Block
languagepy
titleArea cropping: GloFAS Medium-range reforecast 
collapsetrue
## === retrieve GloFAS Medium-Range Reforecast ===

## === subset India, Pakistan, Nepal and Bangladesh region === 


import cdsapi
from datetime import datetime, timedelta


def get_monthsdays():

    start, end = datetime(2019, 1, 1), datetime(2019, 12, 31)
    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

MONTHSDAYS = get_monthsdays()

if __name__ == '__main__':
    c = cdsapi.Client()
    
    # user inputs
	BBOX = [40.05 ,59.95, 4.95, 95.05] # North West South East
    YEARS  = ['%d'%(y) for y in range(1999,2019)]
    LEADTIMES = ['%d'%(l) for l in range(24,1128,24)]
    
    # submit request
    for md in MONTHSDAYS:

        month = md[0].lower()
        day = md[1]

        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': BBOX,# < - subset
                'hyear': YEARS,
                'hmonth': month ,
                'hday': day ,
                'leadtime_hour': LEADTIMES,
            },
            f'glofas_reforecast_{month}_{day}.grib')

...

)


Code Block
languagepy
titleArea cropping: GloFAS Seasonal Forecast
collapsetrue
## === retrieve GloFAS Seasonal Forecast ===

## === subset South America/Amazon region === 

import cdsapi


if __name__ == '__main__':
    c = cdsapi.Client()

    YEARS  = ['%d'%(y) for y in range(2020,2022)]


    MONTHS = ['%02d'%(m) for m in range(1,13)]

    LEADTIMES = ['%d'%(l) for l in range(24,2976,24)]
    
    for year in YEARS:

        for month in MONTHS:
            
            c.retrieve(
                'cems-glofas-seasonal',
                {   
                'variable': 'river_discharge_in_the_last_24_hours',
                'format': 'grib',
                'year': year,
                'month': '12' if year == '2020' else month,
                'leadtime_hour': LEADTIMES,
				'area': [ 10.95, -90.95, -30.95, -29.95 ]

                },
                f'glofas_seasonal_{year}_{month}.grib')


Code Block
languagepy
titleArea cropping: GloFAS Seasonal Reforecast

...

collapsetrue
Code Block
languagepy
## === retrieve GloFAS Seasonal Reforecast ===

## === subset South America/Amazon region === 

import cdsapi

if __name__ == '__main__':


    c = cdsapi.Client()

    YEARS  = ['%d'%(y) for y in range(1981,2021)]

    MONTHS = ['january', 'february', 'march', 'april', 'may', 'june', 'july', 'august', 'september', 'october', 'november', 'december']

    LEADTIMES = ['%d'%(l) for l in range(24,2976,24)]
    
    for year in YEARS:
        for month in MONTHS:

            c.retrieve(
                'cems-glofas-seasonal-reforecast',
                {
                    'system_version': 'version_2_2',
                    'variable':'river_discharge_in_the_last_24_hours',
                    'format':'grib',
                    'hydrological_model':'htessel_lisflood',
                    'hyear': year,
                    'hmonth': month,
                    'leadtime_hour': LEADTIMES,
                    'area': [ 10.95, -90.95, -30.95, -29.95 ]
                }, 
                f'glofas_seasonal_reforecast_{year}_{month}.grib')
Code Block
languagepy
titleExtract timeseries: GloFAS Medium-range

...

Reforecast
collapsetrue



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



...