modflowpy / flopy

A Python package to create, run, and post-process MODFLOW-based models.
https://flopy.readthedocs.io
Other
517 stars 313 forks source link

VoronoiGrid creates "polygons" with less than 3 vertex pairs #2073

Closed jlarsen-usgs closed 9 months ago

jlarsen-usgs commented 9 months ago

Describe the bug In certain cases such as fine triangular refinement the VoronoiGrid class creates modflow grid nodes that have either 1) empty iverts, 2) a single ivert (point feature), or 3) 2 iverts (line feature). While the VertexGrid class of FloPy can digest this information, it is invalid for geospatial processing such as raster resampling or creating GeoDataFrames

To Reproduce See the notebook Sagehen_structured_2_voronoi_workshop.ipynb

Additional context A simple check at the end of the tri2vor method resolves this issue. After the code that finds unique iverts and sorts them, a few lines of code can be added to check for nodes with less than 3 iverts, remove them, and remove their xy centers from tri.verts. Here is the proposed solution

        # existing code
        vor_verts = np.array(vor_verts)
        for icell in range(len(vor_iverts)):
            iverts_cell = vor_iverts[icell]
            vor_iverts[icell] = list(
                get_sorted_vertices(iverts_cell, vor_verts)
            )

        # new code to solve issue
        # remove empty polygons/iverts, point, and line freatures and their associated xy centers
        points = list(tri.verts)
        pop_list = []
        for icell, ivlist in enumerate(vor_iverts):
            if len(ivlist) < 3:
                pop_list.append(icell)

        for icell in pop_list[::-1]:
            vor_iverts.pop(icell)
            points.pop(icell)

        points = np.array(points)
langevin-usgs commented 9 months ago

Hey @jlarsen-usgs, it makes sense that we should filter out these invalid cells. I think @wpbonelli also had a fix for the case where unnecessary vertices were placed along a cell edge. Would be great if you submitted a PR with this change. Would also be good in the future if we added some tests for these degenerate cases.

jlarsen-usgs commented 9 months ago

I'll submit a PR with this change today.

wpbonelli commented 9 months ago

The duplicate edge vertex fix went in here since it is relevant to PRT but it is not independently tested yet. I will plan to add some tests.