...
| 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 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") |
Area cropping:
| 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")
|
...