# 

activate 

the local

 environment

. 
conda 

activate myenv

# 

install the 

required 

packages
conda 

install 

cartopy matplotlib

# 

and the 

cdsapi
pip install cdsapi

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')

'

        ],
        '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')

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')

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)

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

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')

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')

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