Unexpected negative values in triangles

[veenstrajelmer]

Plotting contour/contourf for the westernscheldt example model raises "ValueError: repeats may not contain negative values."

Example file and code:

import matplotlib.pyplot as plt
plt.close('all')
import xugrid as xu

file_nc = r'p:\dflowfm\maintenance\JIRA\05000-05999\05477\c103_ws_3d_fourier\DFM_OUTPUT_westerscheldt01_0subst\westerscheldt01_0subst_map.nc'

uds = xu.open_dataset(file_nc)

fig, ax = plt.subplots()
pc = uds['mesh2d_flowelem_bl'].ugrid.plot(ax=ax, linewidth=0.5, edgecolors='face', cmap='jet')
fig, ax = plt.subplots()
pc = uds['mesh2d_flowelem_bl'].ugrid.plot.contourf(ax=ax, levels=11, cmap='jet')

This is where it seems to go wrong: https://github.com/Deltares/xugrid/blob/2125805f399b5bdb4ead59363f12f8813d9c94bf/xugrid/ugrid/connectivity.py#L402 'face_node_connectivity' contains no negative values ncol_per_row.min() = 1 n_triangle_per_row = -1

[Huite]

While debugging this, I came across at least one bug with the creation of a dual grid (the centroidal voronoi tesselation) of only existing centroids (rather that also the projection of the voronoi rays with the mesh exterior edges). I've fixed that, but it doesn't solve the problem.

While trying to make a smaller example, I ran into a different issue -- which maybe points to the real issue.

This works fine:

path = "westerscheldt01_0subst_map.nc"
uds = xu.open_dataset(path)
uda = uds["mesh2d_flowelem_bl"]
uda_small = uda.ugrid.sel(x=slice(35_000.0, 45_000.0))
uda_small.ugrid.plot(cmap="viridis")

But I cannot plot the grid as a line collection:

uda_small.ugrid.grid.plot()

IndexError: index 1004 is out of bounds for axis 0 with size 1004

The reason appears to be that the edge_node_connectivity of the subselection contains higher node indexes than there are nodes in the grid. There are four additional nodes, as two edges appear to be misplaced (off in the middle of nowhere rather than closing off a cell).

The problem disappears if I discard the existing edge_node_connectivity array:

path = "westerscheldt01_0subst_map.nc"
ds = xr.open_dataset(path).drop_vars("mesh2d_edge_nodes")
grid = xu.Ugrid2d.from_dataset(ds)
da = ds["mesh2d_flowelem_bl"]
uda = xu.UgridDataArray(da, grid)
uda.ugrid.plot.contourf(cmap="viridis")

So this has me wondering whether there's a problem with the edge node connectivity in the dataset. Alternatively, I'm making a mistake when deriving the face edge connectivity with an existing edge node connectivity. I think I'll start by writing a validation function for the edge node connectivity.

[Huite]

Took rather a lot of digging, but I found the issue:

There's a weird edge in the edge_node_connectivity, which isn't associated with any face.

fig, ax = plt.subplots()
uda.ugrid.plot(ax=ax)
grid.plot(ax=ax, edgecolor="k")
ax.set_xlim(57_000.0, 61_000.0)
ax.set_ylim(383_000.0, 387_000.0)
ax.set_aspect(1.0)

I use the face_edge_connectivity for a number of things, which requires me to derive a new edge_node_connectivity. A little additional logic is then required to preserve the old edge definition... but that only works if the old and the new match in size! I'll have a think a bit more whether I should error here (which isn't very helpful for a user), or try and salvage.

[Huite]

Trying to solve this internally creates a lot of headache. I've added some checking for the edge_node_connectivity now (9a5bb600490c5807389da6f239ab510ad86c7bc8), it'll raise a ValueError stating which edges are the problem (in this case, 16640). Then, it's a relatively easy manual fix:

# Get rid of the invalid edge.
path = "westerscheldt01_0subst_map.nc"
ds = xr.open_dataset(path)
edge_select = np.full(ds.dims["mesh2d_nEdges"], True, dtype=bool)
edge_select[16640] = False
ds = ds.isel(mesh2d_nEdges=edge_select)

# Then create the UgridDataset.
uds = xu.UgridDataset(ds)
uda = uds["mesh2d_flowelem_bl"]
uda.ugrid.plot.contourf()

But you probably want to investigate upstream why this edge is in here...

[veenstrajelmer]

The resulting error upon opening the model output is "ValueError: edge_node_connectivity is invalid, the following edges are not associated with any face: [16640]". This at least shows that the mapfile is invalid and since this is the only usecase up to now, it is not necessary to implement an automatic fix. Could be reconsidered later.

[veenstrajelmer]

One could argue that opening this dataset would result in a warning instead of an error, since most of the operations on the uds still work even if not all edges are associated with a face.

Furthermore, it would be convenient if the associated array can be requested from a dataset, so hardcoding is not necessary. For instance a check_edges_associated function that has ds as input instead of self and values like is the case in edge_node_connectivity: https://github.com/Deltares/xugrid/blob/5d7bb9cf389ed1da3dbd8fd9a6431c4399959552/xugrid/ugrid/ugrid2d.py#L443-L453

The user can then simply use this code to properly read the dataset:

# Get rid of the invalid edge.
path = "westerscheldt01_0subst_map.nc"
ds = xr.open_dataset(path)
edge_select = xu.check_edges_associated(ds) #the edge_node_connectivity func
if not edge_select.all():
    ds = ds.isel(mesh2d_nEdges=edge_select)

# Then create the UgridDataset.
uds = xu.UgridDataset(ds)
uda = uds["mesh2d_flowelem_bl"]
uda.ugrid.plot.contourf()

Update: The function remove_unassociated_edges() was recently implemented in dfm_tools (https://github.com/Deltares/dfm_tools/pull/470). We might want to consider moving such a function into xarray, although silently dropping edges in xugrid might not be desireable.