Smooth merge

ocsmesh/utils.py

@@ -2391,7 +2391,6 @@ def shape_to_msh_t_2(
     return msht
 
 
-
 def triangulate_polygon(
     shape: Union[Polygon, MultiPolygon, gpd.GeoDataFrame, gpd.GeoSeries],
     aux_pts: Union[np.array, Point, MultiPoint, gpd.GeoDataFrame, gpd.GeoSeries] = None,
@@ -2489,6 +2488,40 @@ def triangulate_polygon(
     return msht_tri
 
 
+def triangulate_polygon_s(
+    shape: Union[Polygon, MultiPolygon, gpd.GeoDataFrame, gpd.GeoSeries],
+    min_int_ang=30
+) -> None:
+    '''Triangulate the input shape smoothly
+
+    Parameters
+    ----------
+    shape : Polygon or MultiPolygon or GeoDataFrame or GeoSeries
+        Input polygon to triangulate
+    min_int_ang : integer, default=30
+        Minimal internal angle for triangulation
+
+    Returns
+    -------
+    jigsaw_msh_t
+        Generated triangulation
+    '''
+    #First triangulation. Smooth but adds extra nodes at boundary
+    min_int_ang = 'q'+str(min_int_ang)
+    mesh1 = triangulate_polygon(shape,opts=['p',min_int_ang])
+    #Find all boundary modes
+    nb = get_boundary_edges(mesh1)
+    nb = np.sort(list(set(nb.ravel())))
+    #Select all internal modes
+    all_pts = mesh1.vert2['coord']
+    arr = np.delete(all_pts, nb, axis=0)
+    del mesh1,all_pts,nb
+    #Second triangulation. Not smooth ('p'), but intenal nodes are passed as aux_pts
+    mesh2 = triangulate_polygon(shape=shape,aux_pts=arr,opts='p')
+
+    return mesh2
+
+
 def filter_el_by_area(mesh,l_limit=0,u_limit=1e-7):
     '''
     Get the elements that have area within limits
@@ -2532,7 +2565,7 @@ def filter_el_by_area(mesh,l_limit=0,u_limit=1e-7):
 def create_patch_mesh(mesh_w_problem,
                       sel_el_dict,
                       mesh_for_patch,
-                      buffer_size=1000,
+                      buffer_size=0.01,
                       crs=CRS.from_epsg(4326)
                       ):
     '''
@@ -2574,9 +2607,9 @@ def create_patch_mesh(mesh_w_problem,
 
     clip_tri = gpd.GeoDataFrame(pd.concat(all_gdf, ignore_index=True))
     clip_tri = clip_tri.dissolve()
-    clip_tri = clip_tri.to_crs(crs=3857)
+    clip_tri.crs = clip_tri.estimate_utm_crs()
     clip_tri.geometry = clip_tri.buffer(buffer_size)
-    clip_tri = clip_tri.to_crs(crs=crs)
+    clip_tri.crs = crs
 
     patch = clip_mesh_by_shape(
                 mesh_for_patch.msh_t,
@@ -2596,7 +2629,8 @@ def clip_mesh_by_mesh(mesh_to_be_clipped: jigsaw_msh_t,
                       fit_inside: bool = False,
                       check_cross_edges: bool =True,
                       adjacent_layers: int = 2,
-                      crs=CRS.from_epsg(4326)
+                      crs=CRS.from_epsg(4326),
+                      buffer_size = None,
                       )-> jigsaw_msh_t:
     '''
     Clip a mesh ("mesh_to_be_clipped") based on the
@@ -2624,6 +2658,11 @@ def clip_mesh_by_mesh(mesh_to_be_clipped: jigsaw_msh_t,
     gdf_clipper = [get_mesh_polygons(mesh_clipper)]
     gdf_clipper = gpd.GeoDataFrame(geometry=gdf_clipper, crs=crs)
 
+    if buffer_size != None:
+        gdf_clipper.crs = gdf_clipper.estimate_utm_crs()
+        gdf_clipper.geometry = gdf_clipper.buffer(buffer_size)
+        gdf_clipper.crs = crs
+
     clipped_mesh = clip_mesh_by_shape(
                 mesh_to_be_clipped,
                 shape=gdf_clipper.unary_union,
@@ -2636,23 +2675,35 @@ def clip_mesh_by_mesh(mesh_to_be_clipped: jigsaw_msh_t,
     return clipped_mesh
 
 
-def create_mesh_from_mesh_diff(mesh_for_domain,
-                               mesh_1,
-                               mesh_2,
-                               crs=CRS.from_epsg(4326)):
+def create_mesh_from_mesh_diff(domain: Union[jigsaw_msh_t,
+                                             Polygon,
+                                             MultiPolygon,
+                                             gpd.GeoDataFrame,
+                                             gpd.GeoSeries,
+                                             list],
+                               mesh_1: jigsaw_msh_t,
+                               mesh_2: jigsaw_msh_t,
+                               crs=CRS.from_epsg(4326),
+                               min_int_ang=None,
+                               buffer_domain = 0.001                            
+) -> jigsaw_msh_t:
     '''
     Create a triangulation for the area correspondent to
     the difference between "mesh_1" and "mesh_2" 
     for the extent defined by "mesh_for_domain"
 
     Parameters
     ----------
-    mesh_for_domain : jigsawpy.msh_t.jigsaw_msh_t
-        mesh used as the extent
+    domain : jigsawpy.msh_t or Polygon or MultiPolygon or GeoDataFrame or GeoSeries
+        extent of entire domain (mesh_1+mesh_2+gap)
     mesh_1 : jigsawpy.msh_t.jigsaw_msh_t
         mesh_1
     mesh_2 : jigsawpy.msh_t.jigsaw_msh_t
         mesh_2
+    min_int_ang : 
+        Minimal internal angle for triangulation
+        default = None, i.e. direct connection between vertices
+        recommended = 30
 
     Returns
     -------
@@ -2663,10 +2714,31 @@ def create_mesh_from_mesh_diff(mesh_for_domain,
     -----
     '''
 
-    gdf_mesh = [get_mesh_polygons(mesh_for_domain)]
-    gdf_mesh = gpd.GeoDataFrame(geometry=gdf_mesh,crs=crs)
-
-    poly_buffer = gdf_mesh.unary_union.difference(
+    if isinstance(domain, (gpd.GeoDataFrame)):
+        domain = domain
+    if isinstance(domain, (gpd.GeoSeries)):
+        domain = gpd.GeoDataFrame(geometry=gpd.GeoSeries(domain))
+    if isinstance(domain, Polygon):
+        domain = MultiPolygon([domain])
+        domain = gpd.GeoDataFrame(index=[0], crs=crs, geometry=[domain])
+    if isinstance(domain, jigsaw_msh_t):
+        domain = [get_mesh_polygons(domain)]
+        domain = gpd.GeoDataFrame(geometry=domain,crs=crs)
+    if isinstance(domain, (list)):
+        domain = pd.concat([gpd.GeoDataFrame\
+                            (geometry=[get_mesh_polygons(i)\
+                                       .buffer(buffer_domain,join_style=2)],\
+                                        crs=crs) for i in domain])
+    if not isinstance(domain, (gpd.GeoDataFrame)):
+        raise ValueError("Input shape must be a gpd.GeoDataFrame!")
+
+    domain = domain.dissolve().explode(index_parts=True)
+    domain.crs = domain.estimate_utm_crs()
+    domain = domain.loc[domain['geometry'].area == max(domain['geometry'].area)]
+    domain.crs = crs
+    domain = gpd.GeoDataFrame(geometry=[domain.unary_union],crs=crs)
+
+    poly_buffer = domain.unary_union.difference(
         gpd.GeoDataFrame(
             geometry=[
                 get_mesh_polygons(mesh_1),
@@ -2678,15 +2750,18 @@ def create_mesh_from_mesh_diff(mesh_for_domain,
     gdf_full_buffer = gpd.GeoDataFrame(
         geometry = [poly_buffer],crs=crs)
 
-    msht_buffer = triangulate_polygon(gdf_full_buffer)
-    msht_buffer.crs = gdf_mesh.crs
+    if min_int_ang == None:
+        msht_buffer = triangulate_polygon(gdf_full_buffer)
+    else:
+        msht_buffer = triangulate_polygon_s(gdf_full_buffer,min_int_ang=min_int_ang)
+    msht_buffer.crs = crs
 
     return msht_buffer
 
 
 def merge_neighboring_meshes(*all_msht):
     '''
-    Get combine meshes whose boundaries match
+    Combine meshes whose boundaries match
     Adapted from:
     https://github.com/SorooshMani-NOAA/river-in-mesh/blob/main/river_in_mesh/mesh.py
 
@@ -2790,11 +2865,13 @@ def merge_neighboring_meshes(*all_msht):
     return msht_combined
 
 
-def fix_small_el(mesh_w_problem,
-                 mesh_for_patch,
-                 l_limit=0,u_limit=1e-7,
-                 buffer_size=1000,
-                adjacent_layers=2):
+def fix_small_el(mesh_w_problem: jigsaw_msh_t,
+                 mesh_for_patch: jigsaw_msh_t,
+                 l_limit: float = 0.,
+                 u_limit: float = 1e-7,
+                 buffer_size: float = 0.01,
+                 adjacent_layers: int =2
+) -> jigsaw_msh_t:
     '''
     Fix spurious small elements (<u_limit=1e-7)
 
@@ -2807,15 +2884,15 @@ def fix_small_el(mesh_w_problem,
 
     Returns
     -------
-    ocsmesh.mesh.mesh.EuclideanMesh2D
+    jigsaw_msh_t
         fixed mesh with no small area elements
 
     Notes
     -----
     Other optimal attributes were determined based on testing
     and they can be changed as needed:
          l_limit=0,u_limit=1e-7,
-         buffer_size=1000,
+         buffer_size=0.001,
          adjacent_layers=2    
 
     '''
@@ -2825,25 +2902,15 @@ def fix_small_el(mesh_w_problem,
                                  l_limit=l_limit,
                                  u_limit=u_limit)
 
-    #creating the 3 mesh to be merged:
+    #creating the mesh to be merged:
     patch = create_patch_mesh(mesh_w_problem,
                               small_el,
                               mesh_for_patch,
                               buffer_size=buffer_size)
-    carved_mesh=clip_mesh_by_mesh(mesh_w_problem.msh_t,
-                                  patch,
-                                  inverse=True,
-                                  fit_inside=False,
-                                  check_cross_edges=True,
-                                  adjacent_layers=adjacent_layers)
-    msht_buffer = create_mesh_from_mesh_diff(mesh_w_problem.msh_t,
-                                             patch,
-                                             carved_mesh)
-
-    # merged mesh:
-    fixed_mesh = merge_neighboring_meshes(patch,
-                                          carved_mesh,
-                                          msht_buffer)
+
+    fixed_mesh = merge_overlapping_meshes([mesh_w_problem.msh_t,patch],
+                                          adjacent_layers=adjacent_layers
+                                          )
 
     #cleaning up mesh:
     fixed_mesh = clip_mesh_by_mesh(fixed_mesh,
@@ -2856,3 +2923,73 @@ def fix_small_el(mesh_w_problem,
     fixed_mesh = cleanup_folded_bound_el(fixed_mesh)
 
     return fixed_mesh
+
+
+def merge_overlapping_meshes(all_msht: list,
+                             adjacent_layers: int = 0,
+                             buffer_size: float = 0.005,
+                             buffer_domain: float = 0.01,
+                             min_int_ang: int = 30,
+                             crs=CRS.from_epsg(4326)
+) -> jigsaw_msh_t:
+    '''
+    Combine meshes whose boundaries match
+
+    Parameters
+    ----------
+    *all_msht : jigsaw_msh_t 
+        list of meshes to be merged,
+        the later on the list the higher the priority
+    adjacent_layers : int
+        # of elements of the low priority mesh that will be
+        deleted passed the overlap
+    buffer_size : float
+        size of the buffer that will be added around the
+        high priority mesh when carving out the low priority mesh
+        (sometimes this necessary if the mesh is too narrow)
+    buffer_domain : float
+        size of the buffer that will be added
+        around the entire mesh domain. 
+        This is necessary for preventing slivers at the mesh boundary intersections
+    min_int_ang : int
+        Minimal internal angle for triangulation
+        default = 30, i.e. triangles will have internal angles of at least 30 deg
+        For no smoothness min_int_ang = None
+
+    Returns
+    -------
+    jigsaw_msh_t
+        final merged mesh
+
+    Notes
+    -----
+    '''
+
+    msht_combined = all_msht[0]
+    for msht in all_msht[1:]:
+        carved_mesh = clip_mesh_by_mesh(msht_combined,
+                                        msht,
+                                        adjacent_layers=adjacent_layers,
+                                        buffer_size=buffer_size
+                                        )
+        buff_mesh = create_mesh_from_mesh_diff([msht_combined,msht],
+                                               carved_mesh,msht,
+                                               min_int_ang=min_int_ang,
+                                               buffer_domain=buffer_domain
+                                               )
+        domain = pd.concat([gpd.GeoDataFrame(geometry=\
+                                             [get_mesh_polygons(i)],crs=crs) for \
+                                                i in [msht_combined,msht]])
+
+        msht_combined = merge_neighboring_meshes(buff_mesh,
+                                                 carved_mesh,
+                                                 msht
+                                                 )
+        msht_combined = clip_mesh_by_shape(msht_combined,
+                                           domain.unary_union
+                                           )
+        del carved_mesh,buff_mesh,domain,msht
+
+    msht_combined = cleanup_folded_bound_el(msht_combined)
+
+    return msht_combined

tests/api/utils.py

@@ -67,7 +67,7 @@ def test_create_patch_mesh(self):
                                 filtered,
                                 mesh_for_patch)
 
-        self.assertEqual(len(patch.tria3), 73)
+        self.assertEqual(len(patch.tria3), 101)
 
     def test_clip_mesh_by_mesh(self):
         p = Path(__file__).parents[1] / "data" / "test_mesh_1.2dm"
@@ -92,7 +92,7 @@ def test_create_mesh_from_mesh_diff(self):
                                                        patch.msh_t,
                                                        carved_mesh)
 
-        self.assertEqual(len(msht_buffer.tria3), 48)
+        self.assertEqual(len(msht_buffer.tria3), 49)
 
     def test_merge_neighboring_meshes(self):
         p0 = Path(__file__).parents[1] / "data" / "test_mesh_1.2dm"
@@ -110,7 +110,7 @@ def test_merge_neighboring_meshes(self):
                                                      carved_mesh,
                                                      msht_buffer.msh_t)
 
-        self.assertEqual(len(merged_mesh.tria3), 1130251)
+        self.assertEqual(len(merged_mesh.tria3), 1130280)
 
     def test_fix_small_el(self):
         p = Path(__file__).parents[1] / "data" / "test_mesh_1.2dm"
@@ -120,7 +120,17 @@ def test_fix_small_el(self):
 
         fixed_mesh = utils.fix_small_el(mesh,mesh_for_patch)
 
-        self.assertEqual(len(fixed_mesh.tria3), 1130213)
+        self.assertEqual(len(fixed_mesh.tria3), 1130876)
+
+    def test_merge_overlapping_meshes(self):
+        p = Path(__file__).parents[1] / "data" / "test_mesh_1.2dm"
+        mesh = Mesh.open(p, crs=4326)
+        p3 = Path(__file__).parents[1] / "data" / "patch.2dm"
+        patch = Mesh.open(p3, crs=4326)
+
+        smooth = utils.merge_overlapping_meshes([mesh.msh_t,patch.msh_t])
+
+        self.assertEqual(len(smooth.tria3), 1130935)
 
 
 class FinalizeMesh(unittest.TestCase):