...
This section provides scripts for both cases and for both CEMS-Flood products, GloFAS and EFAS.
Set up Python environment
If you have not done it yet, create a Python virtual environment.
Prepare and retrieve data
...
to update once cropping works....
Time series extraction:
| Code Block | ||||
|---|---|---|---|---|
| ||||
...
| 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 static file that contains the projected coordinates and replace it them in EFAS, as described in the code block below. |
| Code Block | ||||||
|---|---|---|---|---|---|---|
| ||||||
import xarray as xr
import pandas as pd
from pyproj import Transformer,CRS
parameter = "dis06"
ds = xr.open_dataset("efas_reforecast.grib", engine="cfgrib")
df = pd.read_csv("GRDC.csv")
uparea = xr.open_dataset("ec_uparea4.0.nc") # the upstream area
# 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)
total = len(df)
rows = []
count = 0
for lon, lat, id in zip(df.long, df.lat, df.grdc_no):
x1, y1 = transformer.transform(lon, lat)
extracted = ds.sel(x=x1, y=y1, number = 1, method="nearest")[parameter]
df_temp = extracted.drop_vars(["surface", "number"]).to_dataframe().reset_index()
df_temp["grdc"] = str(id)
df_temp = df_temp.drop("step", axis=1)
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") |
...
| Code Block | ||||||
|---|---|---|---|---|---|---|
| ||||||
import xarray as xr import rioxarray as rio from pyproj import Transformer, CRS import numpy as np # Rhine's basin bounding box bbox coordinates in WGS84. coords = [505.972204450796, 4611.296530871059, 546.450796296530, 1150.871059972204] # NW,E,S,W,EN 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")# 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.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, t_srs) ds = d.rio.clip(bbox) |
...
GloFAS
Remote processing
Time series extraction:
...