Make UGRID in QGIS presentable

[visr]

If we load the results from model.to_xugrid we see https://github.com/Deltares/Ribasim/pull/1369.

Links: look ok, but we only write the direct link and not the linestring that may be following the waterway. (#2178) Nodes: look bad, see #1369. Developments in QGIS could help us here, but for this issues I'm more interested in quick wins by modifying to_xugrid. Perhaps @Huite has more ideas here, but the main one I'm thinking about it to convert Basin / area to a mesh if it exists, and present the results on the mesh.

Having them look good would also be an argument to being able to write them directly from the core, #2198.

[Huite]

Linestrings are easy; you should be able to convert directly using xugrid's from_geodataframe, and then you just need a link to linestring indexing array.

convert Basin / area to a mesh if it exists

Do you have concrete ideas about this? There are two flavors to come to mind:

1. Convert each polygon intact into a UGRID face, i.e. 1:1 mapping.
2. Triangulate each polygon first (cheaply via e.g. earcut), i.e. 1:N mapping.

The first option has some downsides:

• The face_node_connectivity is likely quite irregular and will require a lot of fill values (some polygons might have 10 000 vertices, others maybe only 10).
• UGRID does not support holes and the like.

The second option basically requires https://github.com/JuliaGeometry/EarCut.jl or https://juliageometry.github.io/MeshesDocs/stable/algorithms/discretization/#HeldTriangulation

(Mostly out of personal interest, I tried translating the mapbox C++ implementation to numba accelerated Python. This didn't work out so well; numba seems to have bad ergonomics and performance if you try to hack numpy structured array into an array based doubly-linked list... I'm sure a Julia version would be fine though.)

In terms of QGIS rendering performance, options are probably equivalent since QGIS will triangulate anyway.

[Huite]

Out of curiosity, I just tried the straightforward approach and I'm actually quite pleased!

import xugrid as xu
gdf = xu.data.provinces_nl()
grid = xu.Ugrid2d.from_geodataframe(gdf)
grid.plot()

You see ugly stuff appearing due to Baarle-Hertog in Noord-Brabant, but this is easily circumvented by taking only the exteriors.

import xugrid as xu
import shapely
import numpy as np

gdf = xu.data.provinces_nl()
exterior_polygons = shapely.polygonize(gdf.exterior.to_numpy())
grid = xu.Ugrid2d.from_shapely(np.array(exterior_polygons.geoms))# %%
grid.plot()

Fill values aren't too bad in this case...

(grid.face_node_connectivity == -1).sum(axis=1)
array([380, 488, 436,   0, 430, 136, 248, 413, 168, 294, 437, 385])

[visr]

Oh nice, yeah then we can try option 1 first and consider option 2 if that isn't good enough. I was thinking to first implement in Python in model.to_xugrid, and if we do #2198 we can then use the same structure in Julia.

[Huite]

Plotting via xugrid works fine because xugrid constructs matplotlib polygons, which may be concave.

import xarray as xr
uda = xu.UgridDataArray(xr.DataArray(np.arange(12), dims=(grid.face_dimension,), name="provinces"), grid=grid)
uda.ugrid.plot()
# %%

uda.ugrid.to_netcdf("qgis-test-ugrid.nc")

But I'm getting a blank screen in QGIS. I think in this case, I'm guessing the triangulation fails. Most like because QGIS is expecting convex faces and so it's running a simple and fast triangulation. Cleaning it up with a negative buffer and simplification doesn't help.

[Huite]

I just ran into https://github.com/Deltares/xugrid/issues/344

Which is easily fixed by .astype(int) which I'll add a test for an push, but otherwise this works fine:

gdf = xu.data.provinces_nl()
uda = xu.earcut_triangulate_polygons(gdf, column="level").astype(float)
uda.ugrid.to_netcdf("qgis-test-ugrid.nc")

[Huite]

This earcut method should be quite robust (that's the sales pitch too), so it seems unlikely that this creates new demands on the basin geometries. A downside is that if the polygons are very detailed, you need to duplicate many values; consequently, minimal compression on the netCDF is probably desirable? Also worthwhile to check what the typical number of vertices is for the input; in my memory QGIS/MDAL renders fairly smoothly to about a million faces or so but you notice the strain when it gets bigger than that.