Merge pull request #29 from noaa-ocs-modeling/bugfix/subset_rmhole

Merge fixes for remove hole functionality in subsetting script

environment.yml

@@ -17,7 +17,7 @@ dependencies:
   - utm
   - scipy<1.8
   - numba
-  - numpy>=1.21
+  - numpy>=1.21,<1.23 # numpy.typing, scipy 1.7.x support
   - matplotlib
   - requests
   - tqdm

ocsmesh/cli/subset_n_combine.py

@@ -115,31 +115,6 @@ def run(self, args):
             wind_speed, num_buffer_layers, rel_island_area_min,
             out_dir, crs, out_all)
 
-
-    def _remove_holes(self, poly):
-        if isinstance(poly, MultiPolygon):
-            return MultiPolygon([self._remove_holes(p) for p in poly.geoms])
-
-        if poly.interiors:
-            return Polygon(poly.exterior)
-
-        return poly
-
-
-    def _remove_holes_by_relative_size(self, poly, rel_size):
-        if isinstance(poly, MultiPolygon):
-            return MultiPolygon([
-                self._remove_holes_by_relative_size(p, rel_size) for p in poly.geoms])
-
-        if poly.interiors:
-            ref_area = poly.area
-            new_interiors = [
-                    intr for intr in poly.interiors
-                    if Polygon(intr).area / ref_area > rel_size]
-            return Polygon(poly.exterior, new_interiors)
-
-        return poly
-
     def _get_largest_polygon(self, mpoly):
         if isinstance(mpoly, Polygon):
             return mpoly
@@ -199,7 +174,7 @@ def _calculate_clipping_polygon(
                         upstream_size_max))
 
         # Islands are not of intereset when clipping high and low-res meshes
-        clip_poly = self._remove_holes(unary_union([
+        clip_poly = utils.remove_holes(unary_union([
             clip_poly_draft, *poly_upstreams]))
 
         return clip_poly
@@ -583,7 +558,7 @@ def _main(
 
         _logger.info("Calculate clipped polygons...")
         start = time()
-        poly_clip_hires_0 = self._remove_holes(
+        poly_clip_hires_0 = utils.remove_holes(
                 utils.get_mesh_polygons(jig_clip_hires_0))
         poly_clip_lowres_0 = utils.get_mesh_polygons(jig_clip_lowres_0)
         _logger.info(f"Done in {time() - start} sec")
@@ -614,7 +589,7 @@ def _main(
         poly_seam_5, jig_clip_lowres = self._add_overlap_to_polygon(jig_clip_lowres_0, poly_seam_4)
 
         # Cleanup buffer shape
-        poly_seam_6 = self._remove_holes_by_relative_size(
+        poly_seam_6 = utils.remove_holes_by_relative_size(
                 poly_seam_5, rel_island_area_min)
 
         poly_seam_7 = utils.drop_extra_vertex_from_polygon(poly_seam_6)
@@ -623,7 +598,7 @@ def _main(
 
         _logger.info("Calculate reclipped polygons...")
         start = time()
-        poly_clip_hires = self._remove_holes(
+        poly_clip_hires = utils.remove_holes(
                 utils.get_mesh_polygons(jig_clip_hires))
         poly_clip_lowres = utils.get_mesh_polygons(jig_clip_lowres)
         _logger.info(f"Done in {time() - start} sec")

ocsmesh/utils.py

@@ -19,7 +19,7 @@
         Polygon, MultiPolygon,
         box, GeometryCollection, Point, MultiPoint,
         LineString, LinearRing)
-from shapely.ops import polygonize, linemerge
+from shapely.ops import polygonize, linemerge, unary_union
 import geopandas as gpd
 import utm
 
@@ -1755,3 +1755,106 @@ def drop_extra_vertex_from_polygon(
         poly_seam_list.append(Polygon(extr, inters))
 
     return MultiPolygon(poly_seam_list)
+
+
+def remove_holes(
+    poly: Union[Polygon, MultiPolygon]
+) -> Union[Polygon, MultiPolygon]:
+    '''Remove holes from the input polygon(s)
+
+    Given input `Polygon` or `MultiPolygon`, remove all the geometric
+    holes and return a new shape object.
+
+    Parameters
+    ----------
+    poly : Polygon or MultiPolygon
+        The input shape from which the holes are removed
+
+    Returns
+    -------
+    Polygon or MultiPolygon
+        The resulting (multi)polygon after removing the holes
+
+    See Also
+    --------
+    remove_holes_by_relative_size :
+        Remove all the whole smaller than given size from the input shape
+
+    Notes
+    -----
+    For a `Polygon` with no holes, this function returns the original
+    object. For `MultiPolygon` with no holes, the return value is a
+    `unary_union` of all the underlying `Polygon`s.
+    '''
+
+    if isinstance(poly, MultiPolygon):
+        return unary_union([remove_holes(p) for p in poly.geoms])
+
+    elif not isinstance(poly, Polygon):
+        raise ValueError(
+            "The input must be either a `Polygon` or `MultiPolygon`:"
+            + f"\tType: {type(poly)}"
+        )
+
+    if poly.interiors:
+        return Polygon(poly.exterior)
+
+    return poly
+
+
+def remove_holes_by_relative_size(
+    poly: Union[Polygon, MultiPolygon],
+    rel_size:float = 0.1
+) -> Union[Polygon, MultiPolygon]:
+    '''Remove holes from the input polygon(s)
+
+    Given input `Polygon` or `MultiPolygon`, remove all the geometric
+    holes that are smaller than the input relative size `rel_size`
+    and return a new shape object.
+
+    Parameters
+    ----------
+    poly : Polygon or MultiPolygon
+        The input shape from which the holes are removed
+    rel_size : float, default=0.1
+        Maximum ratio of a hole area to the area of polygon
+
+    Returns
+    -------
+    Polygon or MultiPolygon
+        The resulting (multi)polygon after removing the holes
+
+    See Also
+    --------
+    remove_holes :
+        Remove all the whole from the input shape
+        
+    Notes
+    -----
+    For a `Polygon` with no holes, this function returns the original
+    object. For `MultiPolygon` with no holes, the return value is a
+    `unary_union` of all the underlying `Polygon`s.
+
+    If `rel_size=1` is specified the result is the same as 
+    `remove_holes` function, except for the additional cost of 
+    calculating the areas.
+    '''
+
+    if isinstance(poly, MultiPolygon):
+        return unary_union([
+            remove_holes_by_relative_size(p, rel_size) for p in poly.geoms])
+
+    elif not isinstance(poly, Polygon):
+        raise ValueError(
+            "The input must be either a `Polygon` or `MultiPolygon`:"
+            + f"\tType: {type(poly)}"
+        )
+
+    if poly.interiors:
+        ref_area = poly.area
+        new_interiors = [
+                intr for intr in poly.interiors
+                if Polygon(intr).area / ref_area > rel_size]
+        return Polygon(poly.exterior, new_interiors)
+
+    return poly