Page History

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)