proper fix for corrupt grids (#461)

* proper fix for corrupt grids
* added testcase
* Bump version: 0.13.0 → 0.13.1
* updated whatsnew
* updated contributing
* updated notebook

.bumpversion.cfg

@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 0.13.0
+current_version = 0.13.1
 commit = True
 tag = True
 

CITATION.cff

@@ -8,5 +8,5 @@ authors:
     orcid: https://orcid.org/0000-0002-6349-818X
 title: "dfm_tools: A Python package for pre- and postprocessing D-FlowFM model input and output files"
 type: software
-version: 0.13.0
+version: 0.13.1
 doi: https://doi.org/10.5281/zenodo.7857393

dfm_tools/__init__.py

@@ -4,7 +4,7 @@
 
 __author__ = """Jelmer Veenstra"""
 __email__ = '[email protected]'
-__version__ = '0.13.0'
+__version__ = '0.13.1'
 
 from dfm_tools.deprecated import *
 from dfm_tools.errors import *

dfm_tools/meshkernel_helpers.py

@@ -145,46 +145,31 @@ def xugrid_remove_noncontiguous(grid):
     if remove_noncontiguous:
         xu_grid = xugrid_remove_noncontiguous(xu_grid)
 
-    #convert to dataset
+    #convert 0-based to 1-based indices for connectivity variables like face_node_connectivity
     xu_grid_ds = xu_grid.to_dataset()
+    xu_grid_ds = xr.decode_cf(xu_grid_ds) #decode_cf is essential since it replaces fillvalues with nans
+    ds_idx = xu_grid_ds.filter_by_attrs(start_index=0)
+    for varn_conn in ds_idx.data_vars:
+        xu_grid_ds[varn_conn] += 1 #from startindex 0 to 1 (fillvalues are now nans)
+        xu_grid_ds[varn_conn].attrs["start_index"] += 1
+        xu_grid_ds[varn_conn].encoding["_FillValue"] = -1 #can be any value <=0, but not 0 is currently the most convenient for proper xugrid plots.
 
-    #convert 0-based to 1-based grid for connectivity variables like face_node_connectivity #TODO: FM kernel needs 1-based grid, but it should read the attributes instead. Report this (#ug_get_meshgeom, #12, ierr=0. ** WARNING: Could not read mesh face x-coordinates)
-    # TODO: this now results in corrupt grid, workaround is using uds_to_1based_ds after bathy interpolation
-    # ds_idx = xu_grid_ds.filter_by_attrs(start_index=0)
-    # for varn_conn in ds_idx.data_vars:
-    #     xu_grid_ds[varn_conn] += 1
-    #     xu_grid_ds[varn_conn].attrs["_FillValue"] += 1
-    #     xu_grid_ds[varn_conn].attrs["start_index"] += 1
+    # convert to uds and add attrs and crs
+    xu_grid_uds = xu.UgridDataset(xu_grid_ds)
 
-    xu_grid_ds = xu_grid_ds.assign_attrs({#'Conventions': 'CF-1.8 UGRID-1.0 Deltares-0.10', #TODO: conventions come from xugrid, so this line is probably not necessary
+    xu_grid_uds = xu_grid_uds.assign_attrs({#'Conventions': 'CF-1.8 UGRID-1.0 Deltares-0.10', #TODO: conventions come from xugrid, so this line is probably not necessary
                                           'institution': 'Deltares',
                                           'references': 'https://www.deltares.nl',
                                           'source': f'Created with meshkernel {meshkernel.__version__}, xugrid {xu.__version__} and dfm_tools {__version__}',
                                           'history': 'Created on %s, %s'%(dt.datetime.now().strftime('%Y-%m-%dT%H:%M:%S%z'),getpass.getuser()), #TODO: add timezone
                                           })
     #TODO: xugrid overwrites these global attributes upon saving the network file: https://github.com/Deltares/xugrid/issues/111
 
-    xu_grid_uds = xu.UgridDataset(xu_grid_ds)
     add_crs_to_dataset(uds=xu_grid_uds,is_geographic=is_geographic,crs=crs)
 
     return xu_grid_uds
 
 
-def uds_to_1based_ds(uds):
-    ds = uds.ugrid.to_dataset()
-    ds_idx = ds.filter_by_attrs(start_index=0)
-    for varn_conn in ds_idx.data_vars:
-        ds[varn_conn] += 1
-        ds[varn_conn].attrs["_FillValue"] += 1
-        ds[varn_conn].attrs["start_index"] += 1
-
-    # print('_FillValue:', ds.mesh2d_face_nodes.attrs['_FillValue'])
-    # print('start_index:', ds.mesh2d_face_nodes.attrs['start_index'])
-    # print('min:', ds.mesh2d_face_nodes.to_numpy().min())
-    # print('max:', ds.mesh2d_face_nodes.to_numpy().max())
-    return ds
-
-
 def add_crs_to_dataset(uds:(xu.UgridDataset,xr.Dataset),is_geographic:bool,crs:(str,int)):
     """
     

docs/CONTRIBUTING.md

@@ -14,7 +14,7 @@
 
 - download and install Anaconda 64 bit Python 3.9 (or higher) from [anaconda.com](https://www.anaconda.com/distribution/#download-section) (miniconda should also be sufficient, but this is not yet tested). Install it with the recommended settings.
 - open anaconda prompt and navigate to dfm_tools checkout folder, e.g. ``C:\DATA\dfm_tools``
-- ``conda create --name dfm_tools_env python=3.9 git spyder -c conda-forge -y`` (``git`` and ``spyder``, you can also install a newer python version)
+- ``conda create --name dfm_tools_env python=3.9 git spyder -c conda-forge -y`` (``git`` and ``spyder`` are optional, you can also install a newer python version)
 - ``conda activate dfm_tools_env``
 - ``python -m pip install -e .[test]`` (pip developer mode, any updates to the local folder are immediately available in your python. It also installs all requirements via pip, ``[test]`` installs also the developer requirements)
 - ``conda deactivate``

docs/whats-new.md

@@ -1,3 +1,9 @@
+## 0.13.1 (2023-07-12)
+
+### Fix
+- proper fix for FM-compatible network file (1-based face_node_connectivity) by [@veenstrajelmer](https://github.com/veenstrajelmer) in [#461](https://github.com/Deltares/dfm_tools/pull/461)
+
+
 ## 0.13.0 (2023-07-12)
 
 ### Feat

setup.cfg

@@ -1,6 +1,6 @@
 [metadata]
 name = dfm_tools
-version = 0.13.0
+version = 0.13.1
 author = Jelmer Veenstra
 author_email = [email protected]
 description = dfm_tools are pre- and post-processing tools for Delft3D FM

tests/examples_workinprogress/workinprogress_meshkernel_creategrid.py

@@ -151,11 +151,9 @@
 
 #write xugrid grid to netcdf
 netfile = 'englishchannel_net.nc'
-#xu_grid_uds.ugrid.to_netcdf(netfile) # TODO: this currently fails, below is a workaround: https://github.com/Deltares/xugrid/issues/119
-xu_grid_ds = dfmt.uds_to_1based_ds(xu_grid_uds)
-xu_grid_ds.to_netcdf(netfile)
+xu_grid_uds.ugrid.to_netcdf(netfile)
 
 #TODO: update https://github.com/Deltares/dfm_tools/issues/217
 
-#TODO: there is (was) a link missing, maybe due to ldb?
+#TODO: small flow links in resulting grid
 

tests/examples_workinprogress/workinprogress_modelbuilder.py

@@ -100,9 +100,7 @@
 
 #write xugrid grid to netcdf
 netfile  = os.path.join(dir_output, f'{model_name}_net.nc')
-#xu_grid_uds.ugrid.to_netcdf(netfile) # TODO: this currently fails, below is a workaround: https://github.com/Deltares/xugrid/issues/119
-xu_grid_ds = dfmt.uds_to_1based_ds(xu_grid_uds)
-xu_grid_ds.to_netcdf(netfile)
+xu_grid_uds.ugrid.to_netcdf(netfile)
 
 
 #%% generate plifile from grid extent 

tests/test_meshkernel_helpers.py

@@ -5,9 +5,13 @@
 @author: veenstra
 """
 
+import os
 import pytest
 import xugrid as xu
 import dfm_tools as dfmt
+import meshkernel
+import xarray as xr
+import numpy as np
 
 
 @pytest.mark.unittest
@@ -84,3 +88,53 @@ def test_meshkernel_delete_withgdf():
 
     assert len(mk.mesh2d_get().face_nodes) == 17364
 
+
+@pytest.mark.systemtest
+def test_meshkernel_to_UgridDataset():
+    """
+    generate grid with meshkernel. Then convert with `dfmt.meshkernel_to_UgridDataset()` from 0-based to 1-based indexing to make FM-compatible network.
+    assert if _FillValue, start_index, min and max are the expected values, this ensures FM-compatibility
+    """
+    is_geographic = False #TODO: polygon refinement does not work for spherical grids: https://github.com/Deltares/MeshKernelPy/issues/78
+    crs = 'EPSG:28992' #arbitrary non-spherical epsg code
+
+    # create basegrid
+    lon_min, lon_max, lat_min, lat_max = -6, 2, 48.5, 51.2
+    dxy = 0.5
+    make_grid_parameters = meshkernel.MakeGridParameters(origin_x=lon_min,
+                                                         origin_y=lat_min,
+                                                         upper_right_x=lon_max,
+                                                         upper_right_y=lat_max,
+                                                         block_size_x=dxy,
+                                                         block_size_y=dxy)
+    mk = meshkernel.MeshKernel(is_geographic=is_geographic)
+    mk.curvilinear_make_uniform_on_extension(make_grid_parameters)
+    mk.curvilinear_convert_to_mesh2d() #convert to ugrid/mesh2d
+
+    # refine with polygon
+    pol_x = np.array([-5,-4,0,-5], dtype=np.double)
+    pol_y = np.array([49,51,49.5,49], dtype=np.double)
+    geometry_list = meshkernel.GeometryList(pol_x, pol_y)
+    mrp = meshkernel.MeshRefinementParameters()
+    mk.mesh2d_refine_based_on_polygon(polygon=geometry_list, mesh_refinement_params=mrp)
+
+    #convert to xugrid and write to netcdf
+    xu_grid_uds = dfmt.meshkernel_to_UgridDataset(mk=mk, crs=crs)
+    netfile = 'test_startindex_net.nc'
+    xu_grid_uds.ugrid.to_netcdf(netfile)
+
+    # plot
+    # fig,ax = plt.subplots()
+    # xu_grid_uds.grid.plot(ax=ax)
+    # ax.plot(pol_x,pol_y,'r-')
+
+    #assert output grid
+    ds_out = xr.open_dataset(netfile,decode_cf=False).load()
+    ds_out.close()
+    os.remove(netfile)
+    assert ds_out.mesh2d_face_nodes.attrs['_FillValue'] == -1
+    assert ds_out.mesh2d_face_nodes.attrs['start_index'] == 1
+    assert 0 not in ds_out.mesh2d_face_nodes.to_numpy()
+    assert ds_out.mesh2d_face_nodes.to_numpy().min() == -1
+    assert ds_out.mesh2d_face_nodes.to_numpy().max() == 135
+