diff --git a/dfm_tools/hydrolib_helpers.py b/dfm_tools/hydrolib_helpers.py index 175d46041..3375b9e05 100644 --- a/dfm_tools/hydrolib_helpers.py +++ b/dfm_tools/hydrolib_helpers.py @@ -41,12 +41,12 @@ def get_ncbnd_construct(): "axis":"Z", } - ncbnd_construct = {"varn_depth":"depth", - "dimn_depth":"depth", + ncbnd_construct = {"varn_depth":"z", + "dimn_depth":"z", "varn_pointx":"lon", "varn_pointy":"lat", "varn_pointname":"station_id", - "dimn_point":"location", + "dimn_point":"node", "attrs_pointx":attrs_pointx, "attrs_pointy":attrs_pointy, "attrs_depth":attrs_depth, @@ -270,7 +270,7 @@ def DataFrame_to_TimModel(tim_pd, refdate:(dt.datetime, pd.Timestamp, str)): return timmodel -def ForcingModel_to_plipointsDataset(forcingmodel:hcdfm.ForcingModel, npoints=None) -> xr.Dataset: +def ForcingModel_to_plipointsDataset(forcingmodel:hcdfm.ForcingModel, npoints=None, convertnan=False) -> xr.Dataset: if not isinstance(forcingmodel, hcdfm.ForcingModel): raise TypeError('ForcingModel_to_plipointsDataset expects type hcdfm.ForcingModel, not type {type(forcingobj)}') @@ -279,7 +279,7 @@ def ForcingModel_to_plipointsDataset(forcingmodel:hcdfm.ForcingModel, npoints=No varn_pointname = ncbnd_construct['varn_pointname'] plipointsDataset_list = [] for forcinglike in forcingmodel.forcing[:npoints]: - ds_onepoint = forcinglike_to_Dataset(forcinglike) + ds_onepoint = forcinglike_to_Dataset(forcinglike, convertnan=convertnan) ds_onepoint = ds_onepoint.expand_dims(dimn_point) for datavar in ds_onepoint.data_vars: longname = datavar diff --git a/dfm_tools/interpolate_grid2bnd.py b/dfm_tools/interpolate_grid2bnd.py index 205f52054..df848b365 100644 --- a/dfm_tools/interpolate_grid2bnd.py +++ b/dfm_tools/interpolate_grid2bnd.py @@ -545,16 +545,15 @@ def maybe_convert_fews_to_dfmt(ds): ncbnd_construct = get_ncbnd_construct() dimn_point = ncbnd_construct['dimn_point'] varn_pointname = ncbnd_construct['varn_pointname'] - dimn_depth = ncbnd_construct['dimn_depth'] - varn_depth = ncbnd_construct['varn_depth'] - # potential FEWS converts + # convert station data_vars to coords to avoid dfmt issues for var_to_coord in [varn_pointname,'station_names']: if var_to_coord in ds.data_vars: ds = ds.set_coords(var_to_coord) + # assign timeseries_id cf_role to let FM read the station names ds[varn_pointname] = ds[varn_pointname].assign_attrs({'cf_role': 'timeseries_id'}) - # rename data_vars to long_name (e.g. renames so to salinitybnd) + # rename data_vars to long_name (e.g. renames FEWS so to salinitybnd) for datavar in ds.data_vars: if datavar in ['ux','uy']: #TODO: keeping these is consistent with hardcoded behaviour in dfm_tools elsewhere, but not desireable continue @@ -562,19 +561,12 @@ def maybe_convert_fews_to_dfmt(ds): longname = ds[datavar].attrs['long_name'] ds = ds.rename_vars({datavar:longname}) - # rename dims/vars - if 'node' in ds.dims: - ds = ds.rename_dims({'node':dimn_point}) - if 'z' in ds.dims: - ds = ds.rename_dims({'z':dimn_depth}) - if 'z' in ds.variables: - ds = ds.rename_vars({'z':varn_depth}) - # transpose dims #TODO: the order impacts the model results: https://issuetracker.deltares.nl/browse/UNST-7402 # dfmt (arbitrary) dimension ordering is node/time/z # required to reorder to FEWS time/node/z order for comparable results # also time/z/node will result in unexpected results - ds = ds.transpose("time", dimn_point, ...) + if "time" in ds.dims: # check if time dimension is present (astronomic does not have time) + ds = ds.transpose("time", dimn_point, ...) # convert station names to string format (keep attrs and encoding) # also needed to properly export, since we cannot encode it at dtype S1 properly otherwise diff --git a/tests/examples/preprocess_hydrolib_readbc.py b/tests/examples/preprocess_hydrolib_readbc.py index 5200f0839..adbe79658 100644 --- a/tests/examples/preprocess_hydrolib_readbc.py +++ b/tests/examples/preprocess_hydrolib_readbc.py @@ -44,6 +44,8 @@ fig, ax = plt.subplots(figsize=(12, 6)) forcingobj = m.forcing[1] forcing_xr = dfmt.forcinglike_to_Dataset(forcingobj, convertnan=True) + forcing_xr_all = dfmt.ForcingModel_to_plipointsDataset(m, convertnan=True) + forcing_xr2 = forcing_xr_all.sel(node=1) data_vars = list(forcing_xr.data_vars) if forcingobj.function=='t3d': forcing_ts = dfmt.Dataset_to_T3D(forcing_xr)