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

'auxiliary.zip')

And then unzip the auxiliary.zip

$ unzip auxiliary.zip

Archive:  auxiliary.zip
 extracting: sd_2_4.0.nc
 extracting: thmax_2_4.0.nc
 extracting: thmin_2_4.0.nc

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

    # 

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