...
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. |
Remote processing
Time series extraction:
Code Block | ||||||
---|---|---|---|---|---|---|
| ||||||
...
Code Block | ||||||
---|---|---|---|---|---|---|
| ||||||
## === retrieve EFAS Seasonal Forecast === ## === subset Switzerland 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,5160,24)] for year in YEARS: for month in MONTHS: c.retrieve( 'efas-seasonal', { 'variable': 'river_discharge_in_the_last_24_hours', 'format': 'grib', 'model_levels': 'surface_level', 'year': year, 'month': '12' if year == '2020' else month, 'leadtime_hour': LEADTIMES, 'area': [ 47.9, 5.8, 45.7, 10.6 ] }, f'efas_seasonal_{year}_{month}.grib') |
Local processing
Time series extraction:
Info | ||
---|---|---|
| ||
EFAS x and y coordinates, when converted from GRIB to NetCDF, are not projected coordinates but matrix indexes (i, j), It is necessary to download the upstream area that contains the projected coordinates and replace them in EFAS, as described in the code block below. |
...
Code Block | ||||||
---|---|---|---|---|---|---|
| ||||||
import xarray as xr import rioxarray as rio from pyproj import Transformer, CRS import numpy as np # Rhine's basin bounding box coordinates in WGS84. coords = [5.450796, 11.871059, 46.296530, 50.972204] # W,E,S,N # source/target reference systems proj EPSG_4326 = '+proj=longlat +datum=WGS84 +no_defs' EPSG_3035 = '+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs' # EFAS # read EFAS reforecast and the upstream area ds = xr.open_dataset("efas_reforecast.grib", engine="cfgrib") uparea = xr.open_dataset("ec_uparea4.0.nc") # replace x, y coordinates ds["x"] = uparea["x"] ds["y"] = uparea["y"] # add reference system to EFAS dataset ds = ds.rio.write_crs(EPSG_3035) # Function to convert 4 coordinates into a Polygon def bbox_from_wesn(coords, s_srs, t_srs): w, e, s, n = coords transformer = Transformer.from_crs(s_srs, t_srs, always_xy=True) # topleft topleft = transformer.transform(w,n) #bottomleft bottomleft = transformer.transform(w, s) #topright topright = transformer.transform(e,n) # bottomright bottomright = transformer.transform(e,s) bbox = [ { 'type': 'Polygon', 'coordinates': [[ topleft, bottomleft, bottomright, topright, topleft ]] } ] return bbox bbox = bbox_from_wesn(coords, s_srs=EPSG_4326, t_srs=EPSG_3035) ds_clipped = ds.rio.clip(bbox) ds_clipped.to_netcdf("efas_reforecast_rhyne.nc") |
...
GloFAS
Remote processing
Time series extraction:
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 = ['%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') |
...
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: 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') |
Local processing
Time series extraction:
Code Block | ||||||
---|---|---|---|---|---|---|
| ||||||
import xarray as xr import pandas as pd parameter = "dis24" ds = xr.open_dataset("glofas_historical.grib", engine="cfgrib",backend_kwargs={'time_dims':['time']}) df = pd.read_csv("GRDC.csv") total = len(df) rows = [] count = 0 for lon, lat, id in zip(df.long, df.lat, df.grdc_no): extracted = ds.sel(longitude=lon, latitude=lat, method="nearest")[parameter] df_temp = extracted.drop_vars(["surface"]).to_dataframe().reset_index() df_temp["grdc"] = str(id) df_temp = df_temp.set_index(["grdc", "time"]) rows.append(df_temp) count += 1 print(f"progress: {count/total*100} %") out = pd.concat(rows) out.to_csv("extracted.csv", index="grdc") |
...