Skip to content

Commit

Permalink
Add gridtype attribute to REGIS contruction (#294)
Browse files Browse the repository at this point in the history
* Add gridtype attribute to REGIS contruction

* delr and delc required to compute affine

Attempting to obtain ahn with just the layer_model from REGIS instead of the model_ds. Additional properties are needed from layer_model. Review needed! Does this also work for irregular grids?
  • Loading branch information
bdestombe authored Nov 24, 2023
1 parent ebcf075 commit 5db5905
Show file tree
Hide file tree
Showing 2 changed files with 14 additions and 1 deletion.
14 changes: 13 additions & 1 deletion nlmod/dims/resample.py
Original file line number Diff line number Diff line change
Expand Up @@ -494,10 +494,22 @@ def structured_da_to_ds(da, ds, method="average", nodata=np.NaN):
logger.info(f"No crs in da. Setting crs equal to ds: {ds.rio.crs}")
da = da.rio.write_crs(ds.rio.crs)
if ds.gridtype == "structured":
# assume an average delr and delc to calculate the affine
if "extent" in ds.attrs:
# xmin, xmax, ymin, ymax
dx = (ds.attrs["extent"][1] - ds.attrs["extent"][0]) / len(ds.x)
dy = (ds.attrs["extent"][3] - ds.attrs["extent"][2]) / len(ds.y)
elif "delr" in ds.attrs and "delc" in ds.attrs:
dx = ds.attrs["delr"]
dy = ds.attrs["delc"]
else:
raise ValueError(
"No extent or delr and delc in ds. Cannot determine affine."
)
da_out = da.rio.reproject(
dst_crs=ds.rio.crs,
shape=(len(ds.y), len(ds.x)),
transform=get_affine(ds),
transform=get_affine(ds, sx=dx, sy=-dy),
resampling=resampling,
nodata=nodata,
)
Expand Down
1 change: 1 addition & 0 deletions nlmod/read/regis.py
Original file line number Diff line number Diff line change
Expand Up @@ -177,6 +177,7 @@ def get_regis(
variables = variables + ("sdh", "sdv")
ds = ds[list(variables)]

ds.attrs["gridtype"] = "structured"
ds.attrs["extent"] = extent
for datavar in ds:
ds[datavar].attrs["grid_mapping"] = "crs"
Expand Down

0 comments on commit 5db5905

Please sign in to comment.