Merge pull request #299 from Deltares/remove-fill-value

Use a FILL_VALUE=-1 everywhere

docs/changelog.rst

@@ -6,6 +6,27 @@ All notable changes to this project will be documented in this file.
 The format is based on `Keep a Changelog`_, and this project adheres to
 `Semantic Versioning`_.
 
+Unreleased
+----------
+
+Fixed
+~~~~~
+
+- Release 0.12.0 changed the return type of the face node connectivity of
+  :attr:`xugrid.Ugrid2d.voronoi_topology` from a `scipy.sparse.coo_matrix` to
+  an ordinary `np.array` of integers (and similarly for internal voronoi
+  tesselations); this dense array had fill (hard-coded) values of -1,
+  potentially differing from the grid's fill value. This lead to a number of
+  errors for methods relying on voronoi tesselations (such as contour plots)
+  if the fill value of the grid was not -1. Internally, a ``FILL_VALUE = -1``
+  is now used everywhere in connectivity arrays, and fill values are no longer
+  passed for internal methods; a value of -1 is always assumed. When converting
+  the grid (back) to a dataset with :meth:`xugrid.Ugrid1d.to_dataset` or
+  :meth:`xugrid.Ugrid2d.to_dataset`, the fill value is set back to its original
+  value; the fill value is also set when calling
+  :meth:`xugrid.UgridDataArrayAccessor.to_netcdf` or
+  :meth:`xugrid.UgridDatasetAccessor.to_netcdf`.
+
 [0.12.0] 2024-09-03
 -------------------
 

examples-dev/voronoi.py

@@ -19,11 +19,14 @@
 modules, should you not want to rely on more complex dependencies such as
 ``xugrid`` and ``xarray``.
 """
+# %%
 import matplotlib.pyplot as plt
 import matplotlib.tri as mtri
 import numpy as np
 from matplotlib.collections import LineCollection, PolyCollection
 
+from xugrid.constants import FILL_VALUE  # equals -1
+
 # %%
 # From xugrid, we need only import the ``connectivity`` and ``voronoi``
 # modules. The functions in these modules depend only on ``numpy`` and
@@ -63,12 +66,12 @@ def generate_disk(partitions: int, depth: int):
     return np.column_stack((x, y)), triang.triangles
 
 
-def edge_plot(vertices, edge_nodes, ax, fill_value=-1, **kwargs):
+def edge_plot(vertices, edge_nodes, ax, **kwargs):
     n_edge = len(edge_nodes)
     edge_coords = np.empty((n_edge, 2, 2), dtype=float)
     node_0 = edge_nodes[:, 0]
     node_1 = edge_nodes[:, 1]
-    valid = (node_0 != fill_value) & (node_1 != fill_value)
+    valid = (node_0 != FILL_VALUE) & (node_1 != FILL_VALUE)
     node_0 = node_0[valid]
     node_1 = node_1[valid]
     edge_coords[:, 0, 0] = vertices[node_0, 0]
@@ -81,10 +84,10 @@ def edge_plot(vertices, edge_nodes, ax, fill_value=-1, **kwargs):
     return primitive
 
 
-def face_plot(vertices, face_nodes, ax, fill_value=-1, **kwargs):
+def face_plot(vertices, face_nodes, ax, **kwargs):
     vertices = vertices[face_nodes]
     # Replace fill value; PolyCollection ignores NaN.
-    vertices[face_nodes == fill_value] = np.nan
+    vertices[face_nodes == FILL_VALUE] = np.nan
     collection = PolyCollection(vertices, **kwargs)
     primitive = ax.add_collection(collection)
     ax.autoscale()
@@ -107,12 +110,12 @@ def comparison_plot(
         sharex=True,
     )
 
-    edges0, _ = connectivity.edge_connectivity(faces0, -1)
+    edges0, _ = connectivity.edge_connectivity(faces0)
     edge_plot(vertices0, edges0, ax0, colors="black")
     ax0.scatter(*centroids0.T, color="red")
     ax0.scatter(*vertices0.T, color="black")
 
-    edges1, _ = connectivity.edge_connectivity(faces1, -1)
+    edges1, _ = connectivity.edge_connectivity(faces1)
     edge_plot(vertices0, edges0, ax1, colors="black")
     edge_plot(vertices1, edges1, ax1, colors="red")
 
@@ -126,13 +129,11 @@ def comparison_plot(
 # %%
 # Let's start by generating a simple unstructured mesh and use only its
 # centroids to generate a voronoi tesselation.
-#
-# Note: ``-1`` functions as the fill value in this example.
 
 vertices, faces = generate_disk(5, 2)
 centroids = vertices[faces].mean(axis=1)
 
-node_face_connectivity = connectivity.invert_dense_to_sparse(faces, -1)
+node_face_connectivity = connectivity.invert_dense_to_sparse(faces)
 voronoi_vertices, voronoi_faces, face_index, _ = voronoi.voronoi_topology(
     node_face_connectivity,
     vertices,
@@ -158,17 +159,14 @@ def comparison_plot(
 # The ``voronoi_topology`` is capable of preserving the exterior exactly, but
 # this requires more topological information.
 
-edge_node_connectivity, face_edge_connectivity = connectivity.edge_connectivity(
-    faces, -1
-)
-edge_face_connectivity = connectivity.invert_dense(face_edge_connectivity, -1)
+edge_node_connectivity, face_edge_connectivity = connectivity.edge_connectivity(faces)
+edge_face_connectivity = connectivity.invert_dense(face_edge_connectivity)
 voronoi_vertices, voronoi_faces, face_index, _ = voronoi.voronoi_topology(
     node_face_connectivity,
     vertices,
     centroids,
     edge_face_connectivity=edge_face_connectivity,
     edge_node_connectivity=edge_node_connectivity,
-    fill_value=-1,
     add_exterior=True,
     add_vertices=True,
 )
@@ -190,18 +188,15 @@ def comparison_plot(
 faces = connectivity.renumber(new)
 centroids = vertices[faces].mean(axis=1)
 
-node_face_connectivity = connectivity.invert_dense_to_sparse(faces, -1)
-edge_node_connectivity, face_edge_connectivity = connectivity.edge_connectivity(
-    faces, -1
-)
-edge_face_connectivity = connectivity.invert_dense(face_edge_connectivity, -1)
+node_face_connectivity = connectivity.invert_dense_to_sparse(faces)
+edge_node_connectivity, face_edge_connectivity = connectivity.edge_connectivity(faces)
+edge_face_connectivity = connectivity.invert_dense(face_edge_connectivity)
 voronoi_vertices, voronoi_faces, face_index, _ = voronoi.voronoi_topology(
     node_face_connectivity,
     vertices,
     centroids,
     edge_face_connectivity=edge_face_connectivity,
     edge_node_connectivity=edge_node_connectivity,
-    fill_value=-1,
     add_exterior=True,
     add_vertices=True,
 )
@@ -217,18 +212,15 @@ def comparison_plot(
 # of the original mesh altogether. We still add an orthogonal projection of
 # every centroid to exterior edges.
 
-node_face_connectivity = connectivity.invert_dense_to_sparse(faces, -1)
-edge_node_connectivity, face_edge_connectivity = connectivity.edge_connectivity(
-    faces, -1
-)
-edge_face_connectivity = connectivity.invert_dense(face_edge_connectivity, -1)
+node_face_connectivity = connectivity.invert_dense_to_sparse(faces)
+edge_node_connectivity, face_edge_connectivity = connectivity.edge_connectivity(faces)
+edge_face_connectivity = connectivity.invert_dense(face_edge_connectivity)
 voronoi_vertices, voronoi_faces, face_index, _ = voronoi.voronoi_topology(
     node_face_connectivity,
     vertices,
     centroids,
     edge_face_connectivity=edge_face_connectivity,
     edge_node_connectivity=edge_node_connectivity,
-    fill_value=-1,
     add_exterior=True,
     add_vertices=False,
 )
@@ -240,18 +232,15 @@ def comparison_plot(
 # exactly. Alternatively, we can choose to skip the exterior vertex if it
 # creates a concave face:
 
-node_face_connectivity = connectivity.invert_dense_to_sparse(faces, -1)
-edge_node_connectivity, face_edge_connectivity = connectivity.edge_connectivity(
-    faces, -1
-)
-edge_face_connectivity = connectivity.invert_dense(face_edge_connectivity, -1)
+node_face_connectivity = connectivity.invert_dense_to_sparse(faces)
+edge_node_connectivity, face_edge_connectivity = connectivity.edge_connectivity(faces)
+edge_face_connectivity = connectivity.invert_dense(face_edge_connectivity)
 voronoi_vertices, voronoi_faces, face_index, _ = voronoi.voronoi_topology(
     node_face_connectivity,
     vertices,
     centroids,
     edge_face_connectivity=edge_face_connectivity,
     edge_node_connectivity=edge_node_connectivity,
-    fill_value=-1,
     add_exterior=True,
     add_vertices=True,
     skip_concave=True,
@@ -267,44 +256,41 @@ def comparison_plot(
     vertices,
     centroids,
 )
-edges0, _ = connectivity.edge_connectivity(faces0, -1)
+edges0, _ = connectivity.edge_connectivity(faces0)
 
 nodes1, faces1, face_index1, _ = voronoi.voronoi_topology(
     node_face_connectivity,
     vertices,
     centroids,
     edge_face_connectivity=edge_face_connectivity,
     edge_node_connectivity=edge_node_connectivity,
-    fill_value=-1,
     add_exterior=True,
     add_vertices=False,
 )
-edges1, _ = connectivity.edge_connectivity(faces1, -1)
+edges1, _ = connectivity.edge_connectivity(faces1)
 
 nodes2, faces2, _, _ = voronoi.voronoi_topology(
     node_face_connectivity,
     vertices,
     centroids,
     edge_face_connectivity=edge_face_connectivity,
     edge_node_connectivity=edge_node_connectivity,
-    fill_value=-1,
     add_exterior=True,
     add_vertices=True,
 )
-edges2, _ = connectivity.edge_connectivity(faces2, -1)
+edges2, _ = connectivity.edge_connectivity(faces2)
 
 nodes3, faces3, face_index3, node_map3 = voronoi.voronoi_topology(
     node_face_connectivity,
     vertices,
     centroids,
     edge_face_connectivity=edge_face_connectivity,
     edge_node_connectivity=edge_node_connectivity,
-    fill_value=-1,
     add_exterior=True,
     add_vertices=True,
     skip_concave=True,
 )
-edges3, _ = connectivity.edge_connectivity(faces3, -1)
+edges3, _ = connectivity.edge_connectivity(faces3)
 
 fig, axes = plt.subplots(
     nrows=1,
@@ -342,10 +328,10 @@ def comparison_plot(
 # fourth option, since it includes some vertices of the original mesh, which
 # are connected to multiple faces.
 
-triangles0, face_triangles0 = connectivity.triangulate(faces0, -1)
+triangles0, face_triangles0 = connectivity.triangulate(faces0)
 triangulation0 = mtri.Triangulation(nodes0[:, 0], nodes0[:, 1], triangles0)
 
-triangles1, face_triangles1 = connectivity.triangulate(faces1, -1)
+triangles1, face_triangles1 = connectivity.triangulate(faces1)
 triangulation1 = mtri.Triangulation(nodes1[:, 0], nodes1[:, 1], triangles1)
 
 
@@ -422,3 +408,5 @@ def comparison_plot(
 
 ax0.set_xlim(-1.5, 1.5)
 ax0.set_ylim(-1.5, 1.5)
+
+# %%