Edge case: Export of mesh cell polygons

[thwllms]

In testing rashdf on the Coal model geometry from the Kanawha pilot, I got this error:

python ./cli.py mesh_cell_polygons /mnt/c/Users/USTW720127/Downloads/ras-jan-2024/Coal/Coal.g01.hdf /mnt/c/temp/coal-mesh-cell-polygons.shp                               venv-rashdf  15:08:57
Traceback (most recent call last):
  File "/home/thomaswilliams/dev/rashdf/src/./cli.py", line 129, in <module>
    main()
  File "/home/thomaswilliams/dev/rashdf/src/./cli.py", line 125, in main
    export(args)
  File "/home/thomaswilliams/dev/rashdf/src/./cli.py", line 107, in export
    gdf: GeoDataFrame = func()
  File "/home/thomaswilliams/dev/rashdf/src/rashdf/geom.py", line 129, in mesh_cell_polygons
    np.vectorize(
  File "/home/thomaswilliams/dev/rashdf/venv-rashdf/lib/python3.10/site-packages/numpy/lib/function_base.py", line 2372, in __call__
    return self._call_as_normal(*args, **kwargs)
  File "/home/thomaswilliams/dev/rashdf/venv-rashdf/lib/python3.10/site-packages/numpy/lib/function_base.py", line 2365, in _call_as_normal
    return self._vectorize_call(func=func, args=vargs)
  File "/home/thomaswilliams/dev/rashdf/venv-rashdf/lib/python3.10/site-packages/numpy/lib/function_base.py", line 2455, in _vectorize_call
    outputs = ufunc(*inputs)
  File "/home/thomaswilliams/dev/rashdf/src/rashdf/geom.py", line 130, in <lambda>
    lambda face_id_list: polygonize(
  File "/home/thomaswilliams/dev/rashdf/venv-rashdf/lib/python3.10/site-packages/shapely/geometry/base.py", line 997, in __getitem__
    raise IndexError("index out of range")
IndexError: index out of range

After a bunch of debugging I figured out that this was caused by overlapping cell edges on 3 cells (out of 20k).

Each of the offending cells is on the perimeter of the 2D flow area. The Shapely polygonize function does not allow lines that overlap like this.

When exporting cell polygons from RAS Mapper, the result is a Shapefile that replicates these overlaps as Polygon features (not MultiPolygon). I think we need to replicate that.

Even though polygonize doesn't like this sort of topology issue, the Shapely Polygon object seems to allow it.

Here's how I think we should approach this for a given cell:

1. Attempt to use polygonize to create a cell Polygon object from its faces
2. If polygonize results in a GeometryCollection with len(result.geoms) == 0 (which is what causes the error above), then we need to build a Polygon object manually
3. Given the faces of a cell, put the coordinates in order:
   1. Start with the first face in the sequence, adding its coordinates to a list
   2. Get the end point of the first face
   3. Check start points of each of the remaining faces and see if they match the end point above
      • If there's a match, add the matching face's coordinates to the list
      • If there isn't a match, check the end points of the remaining faces and add the matching face's reversed coordinates to the list
   4. Continue with the matched face and the remaining faces. Repeat until all faces' coordinates are added to the list.
4. Convert the list of ordered coordinates to a Polygon object.