Support a start_index of 1

docs/api.rst

@@ -195,6 +195,8 @@ UGRID1D Topology
     Ugrid1d.sizes
     Ugrid1d.attrs
     Ugrid1d.coords
+    Ugrid1d.fill_value
+    Ugrid1d.start_index
 
     Ugrid1d.n_node
     Ugrid1d.node_dimension
@@ -258,6 +260,8 @@ UGRID2D Topology
     Ugrid2d.sizes
     Ugrid2d.attrs
     Ugrid2d.coords
+    Ugrid2d.fill_value
+    Ugrid2d.start_index
 
     Ugrid2d.n_node
     Ugrid2d.node_dimension

docs/changelog.rst

@@ -26,6 +26,29 @@ Fixed
   value; the fill value is also set when calling
   :meth:`xugrid.UgridDataArrayAccessor.to_netcdf` or
   :meth:`xugrid.UgridDatasetAccessor.to_netcdf`.
+
+Added
+~~~~~
+
+- :class:`xugrid.Ugrid1d` and :class:`xugrid.Ugrid2d` now take an optional
+  ``start_index`` which controls the start index for the UGRID connectivity
+  arrays.
+- :attr:`xugrid.Ugrid1d.fill_value`, :attr:`xugrid.Ugrid1d.start_index`,
+  :attr:`xugrid.Ugrid2d.fill_value`, and :attr:`xugrid.Ugrid2d.start_index`,
+  have been added to get and set the fill value and start index for the UGRID
+  connectivity arrays. (Internally, every array is 0-based, and has a fill
+  value of -1.)
+
+Changed
+~~~~~~~
+
+- :class:`xugrid.Ugrid1d` and :class:`xugrid.Ugrid2d` will generally preserve
+  the fill value and start index of grids when roundtripping from and to xarray
+  Dataset. An exception is when the start index or fill value varies per
+  connectivity: ``xugrid`` will enforce a single start index and a single fill
+  value per grid. In case of inconsistent values across connectivity arrays,
+  the values associated with the core connectivity are used: for Ugrid2d, this
+  is the face node connectivity.
 
 [0.12.0] 2024-09-03
 -------------------

tests/test_ugrid1d.py

@@ -110,6 +110,11 @@ def test_ugrid1d_properties():
     assert np.array_equal(coords[grid.node_dimension], grid.node_coordinates)
     assert np.array_equal(coords[grid.edge_dimension], grid.edge_coordinates)
 
+    with pytest.raises(ValueError, match="start_index must be 0 or 1, received: 2"):
+        grid.start_index = 2
+    grid.start_index = 1
+    assert grid._start_index == 1
+
 
 def test_ugrid1d_egde_bounds():
     grid = grid1d()

tests/test_ugrid2d.py

@@ -167,6 +167,11 @@ def test_ugrid2d_properties():
     assert np.array_equal(coords[grid.edge_dimension], grid.edge_coordinates)
     assert np.array_equal(coords[grid.face_dimension], grid.face_coordinates)
 
+    with pytest.raises(ValueError, match="start_index must be 0 or 1, received: 2"):
+        grid.start_index = 2
+    grid.start_index = 1
+    assert grid._start_index == 1
+
 
 def test_validate_edge_node_connectivity():
     # Full test at test_connectivity
@@ -332,6 +337,24 @@ def test_ugrid2d_dataset_no_mutation():
     assert ds.identical(reference)
 
 
+@pytest.mark.parametrize("edge_start_index", [0, 1])
+@pytest.mark.parametrize("face_start_index", [0, 1])
+def test_ugrid2d_from_dataset__different_start_index(
+    face_start_index, edge_start_index
+):
+    grid = grid2d()
+    ds = grid.to_dataset(optional_attributes=True)  # include edge_nodes
+    faces = ds["mesh2d_face_nodes"].to_numpy()
+    faces[faces != -1] += face_start_index
+    ds["mesh2d_face_nodes"].attrs["start_index"] = face_start_index
+    ds["mesh2d_edge_nodes"] += edge_start_index
+    ds["mesh2d_edge_nodes"].attrs["start_index"] = edge_start_index
+    new = xugrid.Ugrid2d.from_dataset(ds)
+    assert new.start_index == face_start_index
+    assert np.array_equal(new.face_node_connectivity, grid.face_node_connectivity)
+    assert np.array_equal(new.edge_node_connectivity, grid.edge_node_connectivity)
+
+
 def test_ugrid2d_from_meshkernel():
     # Setup a meshkernel Mesh2d mimick
     class Mesh2d(NamedTuple):

tests/test_ugrid_dataset.py

@@ -1191,19 +1191,48 @@ def test_fm_fillvalue_startindex_isel():
     # xugrid 0.6.0 raises "ValueError: Invalid edge_node_connectivity"
     uds.isel({uds.grid.face_dimension: [1]})
 
+
+def test_alternative_fill_value_start_index():
+    uds = get_ugrid_fillvaluem999_startindex1_uds()
     # Check internal fill value. Should be FILL_VALUE
     grid = uds.ugrid.grid
+    assert grid.face_node_connectivity.dtype == "int64"
+    assert grid.start_index == 1
+    assert grid.fill_value == -999
     assert (grid.face_node_connectivity != -999).all()
     gridds = grid.to_dataset()
     # Should be set back to the origina fill value.
-    assert (gridds["mesh2d_face_nodes"] != xugrid.constants.FILL_VALUE).all()
+    faces = gridds["mesh2d_face_nodes"]
+    assert faces.attrs["start_index"] == 1
+    uniq = np.unique(faces)
+    assert uniq[0] == -999
+    assert uniq[1] == 1
 
     # And similarly for the UgridAccessors.
     ds = uds.ugrid.to_dataset()
-    assert (ds["mesh2d_face_nodes"] != xugrid.constants.FILL_VALUE).all()
+    faces = ds["mesh2d_face_nodes"]
+    assert faces.attrs["start_index"] == 1
+    uniq = np.unique(faces)
+    assert uniq[0] == -999
+    assert uniq[1] == 1
 
     ds_uda = uds["mesh2d_facevar"].ugrid.to_dataset()
-    assert (ds_uda["mesh2d_face_nodes"] != xugrid.constants.FILL_VALUE).all()
+    faces = ds_uda["mesh2d_face_nodes"]
+    assert faces.attrs["start_index"] == 1
+    uniq = np.unique(faces)
+    assert uniq[0] == -999
+    assert uniq[1] == 1
+
+    # Alternative value
+    grid.start_index = 0
+    grid.fill_value = -2
+    gridds = grid.to_dataset()
+    # Should be set back to the origina fill value.
+    faces = gridds["mesh2d_face_nodes"]
+    assert faces.attrs["start_index"] == 0
+    uniq = np.unique(faces)
+    assert uniq[0] == -2
+    assert uniq[1] == 0
 
 
 def test_fm_facenodeconnectivity_fillvalue():

xugrid/ugrid/ugrid1d.py

@@ -51,6 +51,9 @@ class Ugrid1d(AbstractUgrid):
         UGRID topology attributes. Should not be provided together with
         dataset: if other names are required, update the dataset instead.
         A name entry is ignored, as name is given explicitly.
+    start_index: int, 0 or 1, default is 0.
+        Start index of the connectivity arrays. Must match the start index
+        of the provided face_node_connectivity and edge_node_connectivity.
     """
 
     def __init__(
@@ -65,11 +68,13 @@ def __init__(
         projected: bool = True,
         crs: Any = None,
         attrs: Dict[str, str] = None,
+        start_index: int = 0,
     ):
         self.node_x = np.ascontiguousarray(node_x)
         self.node_y = np.ascontiguousarray(node_y)
         self.fill_value = fill_value
-        self.edge_node_connectivity = edge_node_connectivity
+        self.start_index = start_index
+        self.edge_node_connectivity = edge_node_connectivity - self.start_index
         self.name = name
         self.projected = projected
 
@@ -143,8 +148,10 @@ def from_dataset(cls, dataset: xr.Dataset, topology: str = None):
         node_y_coordinates = ds[y_index].astype(FloatDType).to_numpy()
 
         edge_nodes = connectivity["edge_node_connectivity"]
+        fill_value = ds[edge_nodes].encoding.get("_FillValue", -1)
+        start_index = ds[edge_nodes].attrs.get("start_index", 0)
         edge_node_connectivity = cls._prepare_connectivity(
-            ds[edge_nodes], FILL_VALUE, dtype=IntDType
+            ds[edge_nodes], fill_value, dtype=IntDType
         ).to_numpy()
 
         indexes["node_x"] = x_index
@@ -154,13 +161,14 @@ def from_dataset(cls, dataset: xr.Dataset, topology: str = None):
         return cls(
             node_x_coordinates,
             node_y_coordinates,
-            FILL_VALUE,
+            fill_value,
             edge_node_connectivity,
             name=topology,
             dataset=dataset[ugrid_vars],
             indexes=indexes,
             projected=projected,
             crs=None,
+            start_index=start_index,
         )
 
     def _clear_geometry_properties(self):
@@ -229,7 +237,7 @@ def to_dataset(
         data_vars = {
             self.name: 0,
             edge_nodes: xr.DataArray(
-                data=self.edge_node_connectivity,
+                data=self._adjust_connectivity(self.edge_node_connectivity),
                 attrs=edge_nodes_attrs,
                 dims=(self.edge_dimension, "two"),
             ),
@@ -241,7 +249,7 @@ def to_dataset(
 
         dataset = xr.Dataset(data_vars, attrs=attrs)
         if self._dataset:
-            dataset.update(self._dataset)
+            dataset = dataset.merge(self._dataset, compat="override")
         if other is not None:
             dataset = dataset.merge(other)
         if node_x not in dataset or node_y not in dataset:
@@ -569,17 +577,18 @@ def topology_subset(
         new_edges = connectivity.renumber(edge_subset)
         node_x = self.node_x[node_index]
         node_y = self.node_y[node_index]
-        grid = self.__class__(
+        grid = Ugrid1d(
             node_x,
             node_y,
-            self.fill_value,
+            FILL_VALUE,
             new_edges,
             name=self.name,
             indexes=self._indexes,
             projected=self.projected,
             crs=self.crs,
             attrs=self._attrs,
         )
+        self._propagate_properties(grid)
         if return_index:
             indexes = {
                 self.node_dimension: pd.Index(node_index),