...
Example script to crop and extract time series from different GloFAS products:
CDS API
Time series extraction:
Code Block | ||||||
---|---|---|---|---|---|---|
| ||||||
## === 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(],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 MONTHSDAYS = get_monthsdays() if __name__ == '__main__': c = cdsapi.Client() # station # user inputs BBOXcoordinates (lat,lon) COORDS = [40.05 ,59.95, 4.95, 95.05] # North West South East YEARS = ['%d'%(y) for y in range(1999,2019)] { "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 # submit request this case) for md in MONTHSDAYS: month = md[0].lower() day = md[1] # loop over station coordinates for station c.retrieve( 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': ''variable': 'river_discharge_in_the_last_24_hours', 'format': 'grib', 'hydrological_model': 'htessel_lisflood', 'product_type': ['control_reforecast', 'ensemble_perturbed_reforecasts'], 'area': BBOX,# < - subset station_point_coord, 'hyear': YEARSYEAR, 'hmonth': month , 'hday': day , 'leadtime_hour': LEADTIMES, }, f'glofas_reforecast_{station}_{month}_{day}.grib') |
Area cropping:
Code Block | ||||||
---|---|---|---|---|---|---|
| ||||||
## === retrieve GloFAS SeasonalMedium-Range ForecastReforecast === ## === subset South America/Amazon India, Pakistan, Nepal and Bangladesh region === import cdsapi from datetime import datetime, timedelta if def __name__ == '__main__': get_monthsdays(): start, cend = cdsapi.Client() datetime(2019, 1, 1), datetime(2019, 12, 31) YEARS days = ['%d'%(y) for ystart + timedelta(days=i) for i in range(2020,2022)] MONTHS = ['%02d'%(m) for m in range(1,13)] (end - start).days + 1)] LEADTIMESmonthday = ['%d'%(l) for l in range(24,2976,24)] d.strftime("%B-%d").split("-") for d in days if d.weekday() in [0,3] ] for year in YEARS: return monthday MONTHSDAYS = get_monthsdays() if __name__ == '__main__': c for month in MONTHS:= cdsapi.Client() # user inputs BBOX = [40.05 ,59.95, 4.95, 95.05] # North West South East c.retrieve( YEARS = ['%d'%(y) for y in range(1999,2019)] LEADTIMES = 'cems-glofas-seasonal', ['%d'%(l) for l in range(24,1128,24)] # { submit request for md in MONTHSDAYS: month = md[0].lower() day 'variable': 'river_discharge_in_the_last_24_hours',= md[1] c.retrieve( 'format': 'gribcems-glofas-reforecast', { 'year 'system_version': year'version_2_2', 'monthvariable': '12' if year == '2020' else monthriver_discharge_in_the_last_24_hours', 'leadtime_hourformat': LEADTIMES'grib', 'area': [ 10.95, -90.95, -30.95, -29.95 ] }'hydrological_model': 'htessel_lisflood', f'glofas_seasonal_{year}_{month}.grib') | ||||||
Code Block | ||||||
| ||||||
## === 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:product_type': 'control_reforecast', 'area': BBOX,# < - subset 'hyear': YEARS, 'hmonth': month , 'hday': day , 'leadtime_hour': LEADTIMES, }, for month in MONTHS: f'glofas_reforecast_{month}_{day}.grib') |
Code Block | ||||||
---|---|---|---|---|---|---|
| ||||||
## === 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: c.retrieve( 'cems-glofas-seasonal-reforecast', { 'system_version': 'version_2_2', for month in MONTHS: 'variable':'river_discharge_in_the_last_24_hours', c.retrieve( 'format':'grib'cems-glofas-seasonal', { 'hydrological_model':'htessel_lisflood', 'hyearvariable': year'river_discharge_in_the_last_24_hours', 'format': 'grib', 'hmonth 'year': monthyear, 'month': '12' if year == 'leadtime_hour': LEADTIMES2020' else month, 'leadtime_hour': LEADTIMES, 'area': [ 10.95, -90.95, -30.95, -29.95 ] }, f'glofas_seasonal_reforecast_{year}_{month}.grib') |
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 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 = ## === 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,11282976,24)] # loop over date index (just 1 in this case) for mdyear in MONTHSDAYSYEARS: month = md[0].lower() day = md[1] # loop over station coordinates for stationmonth in COORDSMONTHS: station_point_coord = COORDS[station]*2 # coordinates input for the area keyword 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', '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') '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') |
Local machine
EFAS
Warning | ||
---|---|---|
| ||
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. |
...
to update once cropping works....
Local machine
Time series extraction:
We are going to extract EFAS reforecast's timeseries at locations defined by latitude and longitude coordinates from a tiny subset of the GRDC dataset.
...
Code Block | ||||
---|---|---|---|---|
| ||||
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 |
...
Area cropping:
Code Block | ||||||
---|---|---|---|---|---|---|
| ||||||
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") |
...