mesh generation of Harvey with BEST track failed with pygeos.GEOSException: TopologyException

[FariborzDaneshvar-NOAA]

Note added by @SorooshMani-NOAA: This issue occurs on the enhance/subset_triangle branch

@SorooshMani-NOAA this is the error I'm getting during the mesh generation step for the BEST track of harvey:

Traceback (most recent call last):
  File "/opt/conda/envs/ocsmesh/lib/python3.9/runpy.py", line 197, in _run_module_as_main
    return _run_code(code, main_globals, None,
  File "/opt/conda/envs/ocsmesh/lib/python3.9/runpy.py", line 87, in _run_code
    exec(code, run_globals)
  File "/scripts/hurricane_mesh.py", line 546, in <module>
    main(args, [hurrmesh_client, subset_client])
  File "/scripts/hurricane_mesh.py", line 97, in main
    clients_dict[cmd].run(args)
  File "/opt/conda/envs/ocsmesh/lib/python3.9/site-packages/ocsmesh/cli/subset_n_combine.py", line 113, in run
    self._main(
  File "/opt/conda/envs/ocsmesh/lib/python3.9/site-packages/ocsmesh/cli/subset_n_combine.py", line 657, in _main
    utils.finalize_mesh(jig_combined_mesh)
  File "/opt/conda/envs/ocsmesh/lib/python3.9/site-packages/ocsmesh/utils.py", line 428, in finalize_mesh
    boundary_polys = get_mesh_polygons(mesh)
  File "/opt/conda/envs/ocsmesh/lib/python3.9/site-packages/ocsmesh/utils.py", line 244, in get_mesh_polygons
    res_gdf = polys_gdf[polys_gdf.intersects(pnt)]
  File "/opt/conda/envs/ocsmesh/lib/python3.9/site-packages/geopandas/base.py", line 1559, in intersects
    return _binary_op("intersects", self, other, align)
  File "/opt/conda/envs/ocsmesh/lib/python3.9/site-packages/geopandas/base.py", line 59, in _binary_op
    data, index = _delegate_binary_method(op, this, other, align, *args, **kwargs)
  File "/opt/conda/envs/ocsmesh/lib/python3.9/site-packages/geopandas/base.py", line 43, in _delegate_binary_method
    data = getattr(a_this, op)(other, *args, **kwargs)
  File "/opt/conda/envs/ocsmesh/lib/python3.9/site-packages/geopandas/array.py", line 587, in intersects
    return self._binary_method("intersects", self, other)
  File "/opt/conda/envs/ocsmesh/lib/python3.9/site-packages/geopandas/array.py", line 566, in _binary_method
    return getattr(vectorized, op)(left._data, right, **kwargs)
  File "/opt/conda/envs/ocsmesh/lib/python3.9/site-packages/geopandas/_vectorized.py", line 768, in intersects
    return _binary_method("intersects", data, other)
  File "/opt/conda/envs/ocsmesh/lib/python3.9/site-packages/geopandas/_vectorized.py", line 284, in _binary_method
    return getattr(pygeos, op)(left, right, **kwargs)
  File "/opt/conda/envs/ocsmesh/lib/python3.9/site-packages/pygeos/decorators.py", line 80, in wrapped
    return func(*args, **kwargs)
  File "/opt/conda/envs/ocsmesh/lib/python3.9/site-packages/pygeos/predicates.py", line 764, in intersects
    return lib.intersects(a, b, **kwargs)
pygeos.GEOSException: TopologyException: side location conflict at -92.000761999999995 30.268000000000004. This can occur if the input geometry is invalid.
ERROR conda.cli.main_run:execute(49): `conda run python -m hurricane_mesh harvey 2017 subset_n_combine /lustre/static_data/grid/stofs3d_atl_v2.1_eval.gr3 /lustre/static_data/grid/WNAT_1km.14 /lustre/hurricanes/harvey_2017_82b5052c-bd22-4ec3-8aed-1732c675c0e2/windswath --rasters /lustre/static_data/dem/gebco/gebco_2021_sub_ice_topo_n0.0_s-90.0_w-180.0_e-90.0.tif /lustre/static_data/dem/gebco/gebco_2021_sub_ice_topo_n0.0_s-90.0_w-90.0_e0.0.tif /lustre/static_data/dem/gebco/gebco_2021_sub_ice_topo_n0.0_s-90.0_w0.0_e90.0.tif /lustre/static_data/dem/gebco/gebco_2021_sub_ice_topo_n0.0_s-90.0_w90.0_e180.0.tif /lustre/static_data/dem/gebco/gebco_2021_sub_ice_topo_n90.0_s0.0_w-180.0_e-90.0.tif /lustre/static_data/dem/gebco/gebco_2021_sub_ice_topo_n90.0_s0.0_w-90.0_e0.0.tif /lustre/static_data/dem/gebco/gebco_2021_sub_ice_topo_n90.0_s0.0_w0.0_e90.0.tif /lustre/static_data/dem/gebco/gebco_2021_sub_ice_topo_n90.0_s0.0_w90.0_e180.0.tif --out /lustre/hurricanes/harvey_2017_82b5052c-bd22-4ec3-8aed-1732c675c0e2/mesh` failed. (See above for error)

working directory on NHC_COLAB_2: /lustre/hurricanes/harvey_2017_7515c27c-22f2-49e4-b369-ec7fb89e1ee8

Here is the windswath file for this run.
windswath.zip

[yichengt900]

I also encounter a different error Negative elem. areas when running BEST track Florence case using the mesh generated from the new ocsmesh. The workflow crashed at the model spinup stage due to the Negative elem. areas error.

In the error log file, it indicated AQUIRE_HGRID: negative area at 1491512 (see the figure below, the red circle indicates that tiny element; the file is located here: /lustre/hurricanes/florence_2018_nowave_noriver/setup/ensemble.dir/spinup/outputs/nonfatal_000029).

Linking issue: https://github.com/noaa-ocs-modeling/SurgeTeamCoordination/issues/108

[SorooshMani-NOAA]

The Harvey issue seems to be due to triangulation of the buffer region, for some reason the triangulation function also triangulates a "hole" in the buffer zones which causes conflict with existing high-res section of the mesh.

This needs further inspection of the triangulation function to figure out why this happens.

[SorooshMani-NOAA]

Further investigation shows that the feature in previous comment is actually two polygons in a multipolygon that are touching at two points. it's not a single hole touching the sides of a single polygon. This is still valid in shapely, but it causes problem for triangle triangulation package:

[SorooshMani-NOAA]

The main issue seems to be addressed now in the branch (132e114 and 6dbd530), but now I see some strange triangulations like:

This used to work fine, at least partially!

[SorooshMani-NOAA]

@FariborzDaneshvar-NOAA the original issue is resolved. Please test and let me know if you still run into errors. The new issue with strange mesh patches is something I'm still investigating

[SorooshMani-NOAA]

For some reason some of the buffer patches do not have hfun defined!

[SorooshMani-NOAA]

The issue looks to be related to how jigsaw meshes the input! Specifically at:

OCSMesh/ocsmesh/cli/subset_n_combine.py

Lines 289 to 293 in 6dbd530

	mesh_domain_rep = JigsawDriver(
	geom=Geom(buffer_domain, crs=utm),
	hfun=hfun_rep,
	initial_mesh=False
	).run(sieve=0)

The geometry looks like this:

And the hfun like:

But the mesh ends up ignoring the patches:

Note that this mesh was generated in a separate iPython session too (from the saved domain and hfun) and still has this issue. The mesh was generated using ocsmesh's JigsawDriver, maybe let's try to also directly use Jigsaw and see if it differs!

test_files.zip

[SorooshMani-NOAA]

After further investigation, by trying to modify the buffer geom by buffer operation, etc. I still get similar issues. This should be something related to how jigsaw handles separate input polygons.

[SorooshMani-NOAA]

It doesn't make any sense! For the buffer zone polygons (like above) I have this issue with Jigsaw that it doesn't mesh some of them. If I send a subset of these polygons, I get different set of them unmeshed! I also thought it might be that shapes are too complicated and cause truncation errors, so I buffered the shapes by 10000 meters and simplifyd the shapes by tolerance of 5000 meters, but I still get similar issue. I even tried using a background mesh as size function (so that I know the issue is not the size function):

Then I thought it's just Jigsaw cannot handle multiple polygons with holes, but then trying buffered points for polygon it meshes fine (I used the same background mesh as size function):

m_lo = ocsmesh.Mesh.open('path/to/bg/mesh.2dm', crs=4326)
mlo_hfun = ocsmesh.Hfun(m_lo)
mlo_hfun.size_from_mesh()

pts = points([[-92+i, 20] for i in range(10)])
geom = ocsmesh.Geom(MultiPoint(pts).buffer(0.35).difference(MultiPoint(pts).buffer(0.05)), crs=4326)

dr = ocsmesh.driver.JigsawDriver(geom=geom, hfun=mlo_hfun)
m = dr.run(sieve=0)
m.write('multicircle.2dm', format='2dm', overwrite=True)

[SorooshMani-NOAA]

Not addressed as a part of the merge!