Calculate mesh polygons using shapely 2 new functionality

ocsmesh/utils.py

@@ -21,6 +21,7 @@
         Polygon, MultiPolygon,
         box, GeometryCollection, Point, MultiPoint,
         LineString, LinearRing)
+from shapely import build_area, get_parts
 from shapely.ops import polygonize, linemerge, unary_union
 import geopandas as gpd
 import utm
@@ -193,84 +194,10 @@ def get_boundary_segments(mesh):
 def get_mesh_polygons(mesh):
 
 
-    # TODO: Copy mesh?
-    target_mesh = mesh
-    result_polys = []
+    lines = get_boundary_segments(mesh)
+    geom_areas = build_area(lines)
 
-    # 2-pass find, first find using polygons that intersect non-boundary
-    # vertices, then from the rest of the mesh find polygons that
-    # intersect any vertex
-    for find_pass in range(2):
-
-
-        coords = target_mesh.vert2['coord']
-
-        if len(coords) == 0:
-            continue
-
-        boundary_edges = get_boundary_edges(target_mesh)
-
-        lines = get_boundary_segments(target_mesh)
-
-        poly_gen = polygonize(lines)
-        polys = list(poly_gen)
-        polys = sorted(polys, key=lambda p: p.area, reverse=True)
-
-
-        bndry_verts = np.unique(boundary_edges)
-
-        if find_pass == 0:
-            non_bndry_verts = np.setdiff1d(
-                    np.arange(len(coords)), bndry_verts)
-            pnts = MultiPoint(coords[non_bndry_verts])
-        else:
-            pnts = MultiPoint(coords[bndry_verts])
-
-
-        # NOTE: This logic requires polygons to be sorted by area
-        pass_valid_polys = []
-        while len(pnts.geoms) > 0:
-
-
-            idx = np.random.randint(len(pnts.geoms))
-            pnt = pnts.geoms[idx]
-
-            polys_gdf = gpd.GeoDataFrame(
-                {'geometry': polys, 'list_index': range(len(polys))})
-
-
-            res_gdf = polys_gdf[polys_gdf.intersects(pnt)]
-            if len(res_gdf) == 0:
-                # How is this possible for the remaining points not to
-                # intersect with the remaining polys?!
-                pnts = MultiPoint([pt for i, pt in enumerate(pnts.geoms) if i != idx])
-                if pnts.is_empty:
-                    break
-
-                continue
-
-            poly = res_gdf.geometry.iloc[0]
-            polys.pop(res_gdf.iloc[0].list_index)
-
-
-
-            pass_valid_polys.append(poly)
-            pnts = pnts.difference(poly)
-            if pnts.is_empty:
-                break
-            if isinstance(pnts, Point):
-                pnts = MultiPoint([pnts])
-
-
-        result_polys.extend(pass_valid_polys)
-        target_mesh = clip_mesh_by_shape(
-            target_mesh,
-            shape=MultiPolygon(pass_valid_polys),
-            inverse=True, fit_inside=True)
-
-
-
-    return MultiPolygon(result_polys)
+    return MultiPolygon([p for g in geom_areas for p in get_parts(g)])
 
 
 def repartition_features(linestring, max_verts):

pyproject.toml

@@ -25,7 +25,7 @@ dependencies = [
     "jigsawpy", "matplotlib", "netCDF4", "numba",
     "numpy>=1.21", # introduce npt.NDArray
     "pyarrow", "rtree", "pyproj>=3.0", "rasterio", "scipy",
-    "shapely", "typing_extensions", "utm",
+    "shapely>=2", "typing_extensions", "utm",
     ]
 dynamic = ["version"]