944 cmems nc to ini produces incoherent depth value (#955)

* improved test for cmems_nc_to_ini

* fixed incorrect depth definition in 3D initial netcdf files

* updated whatsnew

dfm_tools/modelbuilder.py

@@ -3,12 +3,14 @@
 import pandas as pd
 import dfm_tools as dfmt
 import datetime as dt
+import xarray as xr
 import hydrolib.core.dflowfm as hcdfm
 from hydrolib.core.dimr.models import DIMR, FMComponent, Start
 from hydrolib.core.utils import get_path_style_for_current_operating_system
 from dfm_tools.hydrolib_helpers import get_ncbnd_construct
 from dfm_tools.interpolate_grid2bnd import (ext_add_boundary_object_per_polyline,
                                             open_dataset_extra,
+                                            ds_apply_conversion_dict,
                                             )
 
 __all__ = [
@@ -97,6 +99,11 @@ def preprocess_ini_cmems_to_nc(**kwargs):
 
 
 def cmems_nc_to_ini(ext_old, dir_output, list_quantities, tstart, dir_pattern, conversion_dict=None):
+    """
+    This makes quite specific 3D initial input files based on CMEMS data. delft3dfm is quite picky, 
+    so it works with CMEMS files because they have a depth variable with standard_name='depth'.
+    If this is not present, or dimensions are ordered differently, there will be an error or incorrect model results.
+    """
 
     if conversion_dict is None:
         conversion_dict = dfmt.get_conversion_dict()
@@ -119,35 +126,20 @@ def cmems_nc_to_ini(ext_old, dir_output, list_quantities, tstart, dir_pattern, c
         if quan_bnd=="salinitybnd":
             # 3D initialsalinity/initialtemperature fields are silently ignored
             # initial 3D conditions are only possible via nudging 1st timestep via quantity=nudge_salinity_temperature
-            data_xr = open_dataset_extra(dir_pattern=dir_pattern_one, quantity="salinitybnd",
-                                         tstart=tstart_round, tstop=tstop_round,
-                                         conversion_dict=conversion_dict)
+            data_xr = xr.open_mfdataset(dir_pattern_one)
             ncvarname_tem = get_ncvarname(quantity="temperaturebnd", conversion_dict=conversion_dict)
             dir_pattern_tem = dir_pattern.format(ncvarname=ncvarname_tem)
-            data_xr_tem = open_dataset_extra(dir_pattern=dir_pattern_tem, quantity="temperaturebnd",
-                                             tstart=tstart_round, tstop=tstop_round,
-                                             conversion_dict=conversion_dict)
-            data_xr["temperaturebnd"] = data_xr_tem["temperaturebnd"]
-            data_xr = data_xr.rename_vars({"salinitybnd":"so", "temperaturebnd":"thetao"})
+            data_xr_tem = xr.open_mfdataset(dir_pattern_tem)
+            data_xr["thetao"] = data_xr_tem["thetao"]
             quantity = "nudge_salinity_temperature"
             varname = None
         else:
-            data_xr = open_dataset_extra(dir_pattern=dir_pattern_one, quantity=quan_bnd,
-                                         tstart=tstart_round, tstop=tstop_round,
-                                         conversion_dict=conversion_dict)
+            data_xr = xr.open_mfdataset(dir_pattern_one)
+            data_xr = ds_apply_conversion_dict(data_xr=data_xr, conversion_dict=conversion_dict, quantity=quan_bnd)
             quantity = f'initial{quan_bnd.replace("bnd","")}'
             varname = quantity
             data_xr = data_xr.rename_vars({quan_bnd:quantity})
 
-
-        # open_dataset_extra converted depths from positive down to positive up, including update of the "positive" attribute
-        # TODO: this correctly updated attr negatively impacts model results when using netcdf inifields, so we revert it here
-        # https://issuetracker.deltares.nl/browse/UNST-7455
-        ncbnd_construct = get_ncbnd_construct()
-        varn_depth = ncbnd_construct['varn_depth']
-        if varn_depth in data_xr.coords:
-            data_xr[varn_depth].attrs['positive'] = 'down'
-
         # subset two times. interp to tstart would be the proper way to do it, 
         # but FM needs two timesteps for nudge_salinity_temperature and initial waq vars
         data_xr = data_xr.sel(time=slice(tstart_round, tstop_round))
@@ -160,7 +152,6 @@ def cmems_nc_to_ini(ext_old, dir_output, list_quantities, tstart, dir_pattern, c
         data_xr = data_xr.ffill(dim='latitude').bfill(dim='latitude')
         data_xr = data_xr.ffill(dim='longitude').bfill(dim='longitude')
         data_xr = data_xr.ffill(dim='depth').bfill(dim='depth')
-
 
         print('writing file')
         file_output = os.path.join(dir_output,f"{quantity}_{tstart_str}.nc")

docs/whats-new.md

@@ -15,6 +15,7 @@
 - also apply convert_360to180 to longitude variable in `dfmt.open_dataset_curvilinear()` in [#913](https://github.com/Deltares/dfm_tools/pull/913)
 - prevented adding of time dimension and dropped 0-sized cells in `dfmt.open_dataset_curvilinear()` in [#929](https://github.com/Deltares/dfm_tools/pull/929)
 - dropping lat/lon/verts variables to maintain single source of truth in `dfmt.open_dataset_curvilinear()` in [#931](https://github.com/Deltares/dfm_tools/pull/931)
+- fix for 3D initial fields in `dfmt.cmems_nc_to_ini()` to avoid top layer values over the entire depth in [#955](https://github.com/Deltares/dfm_tools/pull/955)
 
 
 ## 0.24.0 (2024-07-12)

tests/test_interpolate_grid2bnd.py

@@ -58,16 +58,20 @@ def cmems_dataset_notime():
                        [35.792107, np.nan,    np.nan],
                        [35.822628, np.nan,    np.nan]],
 
-                      [[35.781425, np.nan,    np.nan],
-                       [35.792107, np.nan,    np.nan],
-                       [35.789055, np.nan,    np.nan]]])
+                      [[np.nan, np.nan,    np.nan],
+                       [np.nan, np.nan,    np.nan],
+                       [np.nan, np.nan,    np.nan]]])
     ds['so'] = xr.DataArray(so_np,dims=('depth','latitude','longitude'))
     lons = [-9.6, -9.5, -9.4]
     lats = [42.9, 43.0, 43.1]
     depths = [-0.494025, -1.541375, -2.645669, -3.819495, -5.078224]
+
+    depth_attrs = {'positive': 'up'}
+
     ds['longitude'] = xr.DataArray(lons, dims=('longitude'))
     ds['latitude'] = xr.DataArray(lats, dims=('latitude'))
-    ds['depth'] = xr.DataArray(depths, dims=('depth'))
+    ds['depth'] = xr.DataArray(depths, dims=('depth')).assign_attrs(depth_attrs)
+
     return ds
 
 

tests/test_modelbuilder.py

@@ -11,7 +11,7 @@
 import pandas as pd
 import hydrolib.core.dflowfm as hcdfm
 import xarray as xr
-from dfm_tools.hydrolib_helpers import get_ncbnd_construct
+import numpy as np
 from dfm_tools.modelbuilder import get_quantity_list, get_ncvarname
 
 
@@ -81,11 +81,11 @@ def test_cmems_nc_to_ini_midnight_centered(tmp_path):
     assert "thetao" in ds_out.data_vars
     assert ds_out.so.isnull().sum().load() == 0
 
-    ncbnd_construct = get_ncbnd_construct()
-    varn_depth = ncbnd_construct['varn_depth']
-    assert varn_depth in ds_out.coords
-    # the below is inconsistent since depth is actually defined positive up, but FM requires this for inifields: https://issuetracker.deltares.nl/browse/UNST-7455
-    assert ds_out[varn_depth].attrs['positive'] == 'down'
+    # check whether the cmems depth definition comes trough
+    assert 'depth' in ds_out.coords
+    depths_expected = np.array([-0.494025, -1.541375, -2.645669, -3.819495, -5.078224])
+    assert np.allclose(ds_out['depth'].to_numpy(), depths_expected)
+    assert ds_out['depth'].attrs['positive'] == 'up'
 
 
 @pytest.mark.unittest