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

When transforming from lat/lon (source coordinates) to projected LAEA (target coordinates), you need to consider that the number of decimal places of the source coordinates affects the target coordinates precision:

An interval of 0.001 degrees corresponds to about 100 metres in LAEA.

An interval of 0.00001 degrees corresponds to about 1 metre in LAEA.

grdc,time,y,x,latitude,longitude,valid_time,dis06
6118010,2007-03-04,2827500.0,3372500.0,47.82327782578349,-2.7316459165737292,2007-03-04 00:00:00,9.68457
6118010,2007-03-04,2827500.0,3372500.0,47.82327782578349,-2.7316459165737292,2007-03-04 06:00:00,9.5390625
6118010,2007-03-04,2827500.0,3372500.0,47.82327782578349,-2.7316459165737292,2007-03-04 12:00:00,9.783691
6118010,2007-03-04,2827500.0,3372500.0,47.82327782578349,-2.7316459165737292,2007-03-04 18:00:00,9.83252
6118010,2007-03-07,2827500.0,3372500.0,47.82327782578349,-2.7316459165737292,2007-03-07 00:00:00,11.904785
6118010,2007-03-07,2827500.0,3372500.0,47.82327782578349,-2.7316459165737292,2007-03-07 06:00:00,12.73584
6118010,2007-03-07,2827500.0,3372500.0,47.82327782578349,-2.7316459165737292,2007-03-07 12:00:00,12.832031
6118010,2007-03-07,2827500.0,3372500.0,47.82327782578349,-2.7316459165737292,2007-03-07 18:00:00,13.336914
6118015,2007-03-04,2842500.0,3402500.0,48.001715448932764,-2.3682021184138047,2007-03-04 00:00:00,11.508301
6118015,2007-03-04,2842500.0,3402500.0,48.001715448932764,-2.3682021184138047,2007-03-04 06:00:00,11.17334
6118015,2007-03-04,2842500.0,3402500.0,48.001715448932764,-2.3682021184138047,2007-03-04 12:00:00,11.09082
6118015,2007-03-04,2842500.0,3402500.0,48.001715448932764,-2.3682021184138047,2007-03-04 18:00:00,11.20752
6118015,2007-03-07,2842500.0,3402500.0,48.001715448932764,-2.3682021184138047,2007-03-07 00:00:00,14.920898
6118015,2007-03-07,2842500.0,3402500.0,48.001715448932764,-2.3682021184138047,2007-03-07 06:00:00,15.956055
6118015,2007-03-07,2842500.0,3402500.0,48.001715448932764,-2.3682021184138047,2007-03-07 12:00:00,16.127441
6118015,2007-03-07,2842500.0,3402500.0,48.001715448932764,-2.3682021184138047,2007-03-07 18:00:00,15.900879
6118020,2007-03-04,2802500.0,3447500.0,47.71290346585623,-1.6892419697226784,2007-03-04 00:00:00,13.782227

import xarray as xr
from pyproj import Transformer, CRS
import numpy as np

# Rhine's basin bounding box
bbox = [50.972204, 46.296530, 5.450796, 11.871059]  # N,S,W,E

ds = xr.open_dataset("efas_reforecast.grib", engine="cfgrib")
uparea = xr.open_dataset("ec_uparea4.0.nc")

# replace x, y
ds["x"] = uparea["x"]
ds["y"] = uparea["y"]

# define reprojection parameters
laea_proj = CRS.from_proj4(
    "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs"
)
transformer = Transformer.from_crs("epsg:4326", laea_proj, always_xy=True)

we = bbox[2:]
ns = bbox[:2]

we_xy, ns_xy = transformer.transform(we, ns)

we_xy = [np.floor(we_xy[0]), np.ceil(we_xy[1])]
ns_xy = [np.ceil(ns_xy[0]), np.floor(ns_xy[1])]

ds_cropped = ds.sel(
    x=slice(we_xy[0], we_xy[1]), y=slice(ns_xy[0], ns_xy[1])
)

ds_cropped.to_netcdf("efas_forecast_cropped.nc")