Skip to content

Commit

Permalink
Merge pull request #57 from Deltares/map_results
Browse files Browse the repository at this point in the history
fix reading map files with different edge dimension names
  • Loading branch information
DirkEilander authored Feb 23, 2022
2 parents 07d6187 + 979bb12 commit 73534c5
Show file tree
Hide file tree
Showing 3 changed files with 24 additions and 15 deletions.
2 changes: 1 addition & 1 deletion .github/workflows/docs.yml
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,7 @@ jobs:
with:
github_token: ${{ secrets.GITHUB_TOKEN }}
publish_dir: ./docs/_build/html
exclude_assets: '.buildinfo,_sources/*'
exclude_assets: '.buildinfo,_sources/*,_examples/*.ipynb'
destination_dir: ./${{ env.DOC_VERSION }}
keep_files: false
full_commit_message: Deploy ${{ env.DOC_VERSION }} to GitHub Pages
3 changes: 2 additions & 1 deletion docs/changelog.rst
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ All notable changes to this project will be documented in this page.
The format is based on `Keep a Changelog`_, and this project adheres to
`Semantic Versioning`_.

v0.2.1 (18 February 2022)
v0.2.1 (23 February 2022)
-------------------------

Deprecated
Expand All @@ -20,6 +20,7 @@ Bugfix
^^^^^^
- bugfix **setup_p_forcing** to ensure the data is 1D when passed to set_forcing_1d method
- bugfix **setup_p_forcing_from_grid** when aggregating with a multi polygon region.
- bugfix **read_results** with new corner_x/y instead of edge_x/y dimensions in sfincs_map.nc

New
^^^
Expand Down
34 changes: 21 additions & 13 deletions hydromt_sfincs/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -756,17 +756,26 @@ def read_sfincs_map_results(
"""

ds_map = xr.open_dataset(fn_map, chunks={"time": chunksize}, **kwargs)
ds_map = ds_map.set_coords(("x", "y", "edge_x", "edge_y"))
dvars = list(ds_map.data_vars.keys())
edge_dims = [var for var in dvars if (var.endswith("_x") or var.endswith("_y"))]
ds_map = ds_map.set_coords(["x", "y"] + edge_dims)
crs = ds_map["crs"].item() if ds_map["crs"].item() != 0 else crs

if ds_map["inp"].attrs.get("rotation") != 0:
logger.warning("Cannot parse rotated maps. Skip reading sfincs.map.nc")
return xr.Dataset(), xr.Dataset()

# split general+face and edge vars
dvars = list(ds_map.data_vars.keys())
face_vars = [v for v in dvars if "edge_m" not in ds_map[v].dims and v not in drop]
edge_vars = [v for v in dvars if "edge_m" in ds_map[v].dims and v not in drop]
face_vars = list(ds_map.data_vars.keys())
edge_vars = []
if len(edge_dims) == 2:
edge = edge_dims[0][:-2]
face_vars = [
v for v in face_vars if f"{edge}_m" not in ds_map[v].dims and v not in drop
]
edge_vars = [
v for v in face_vars if f"{edge}_m" in ds_map[v].dims and v not in drop
]

# read face vars
face_coords = {
Expand Down Expand Up @@ -795,21 +804,20 @@ def read_sfincs_map_results(
ds_edge = xr.Dataset()
if len(edge_vars) > 0:
edge_coords = {
"edge_x": xr.IndexVariable(
"edge_x", ds_map["edge_x"].isel(edge_n=0).values
f"{edge}_x": xr.IndexVariable(
f"{edge}_x", ds_map[f"{edge}_x"].isel(edge_n=0).values
),
"edge_y": xr.IndexVariable(
"edge_y", ds_map["edge_y"].isel(edge_m=0).values
f"{edge}_y": xr.IndexVariable(
f"{edge}_y", ds_map[f"{edge}_y"].isel(edge_m=0).values
),
}
ds_edge = (
ds_map[edge_vars]
.drop_vars(["edge_x", "edge_y"])
.rename({"edge_n": "edge_y", "edge_m": "edge_x"})
.drop_vars([f"{edge}_x", f"{edge}_y"])
.rename({f"{edge}_n": f"{edge}_y", f"{edge}_m": f"{edge}_x"})
.assign_coords(edge_coords)
.transpose(..., "edge_y", "edge_x")
)
ds_edge.raster.set_spatial_dims(x_dim="edge_x", y_dim="edge_y")
.transpose(..., f"{edge}_y", f"{edge}_x")
).rename({f"{edge}_x": "x", f"{edge}_y": "y"})
ds_edge.raster.set_crs(crs)

return ds_face, ds_edge
Expand Down

0 comments on commit 73534c5

Please sign in to comment.