...
If you have not done it yet, create a Python virtual environment.
Activate the conda environment and install the additional Phython package https://corteva.github.io/rioxarray
Code Block | ||||
---|---|---|---|---|
| ||||
conda install rioxarray |
Prepare and retrieve data
...
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.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_clipped = dds.rio.clip(bbox) |
...
GloFAS
Remote processing
...