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")
|