Skip to content

Commit

Permalink
Merge pull request #174 from felicio93/improve_merge_overlap_meshes
Browse files Browse the repository at this point in the history 
Improve merge overlap meshes
  • Loading branch information
@felicio93
felicio93 authored Sep 12, 2024
2 parents 451d41f + fc7e694 commit a484d9e
Show file tree
Hide file tree
Showing 3 changed files with 80 additions and 27 deletions.
2 changes: 1 addition & 1 deletion ocsmesh/geom/collector.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -225,7 +225,7 @@ def __init__(
clip_shape = ops.transform(
transformer.transform, clip_shape)
try:
in_item.clip(clip_shape)
raster.clip(clip_shape)
except ValueError as err:
# This raster does not intersect shape
_logger.debug(err)
Expand Down
99 changes: 76 additions & 23 deletions ocsmesh/utils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@
RectBivariateSpline, griddata)
from scipy import sparse, constants
from scipy.spatial import cKDTree
from shapely import difference
from shapely import intersection,difference
from shapely.geometry import ( # type: ignore[import]
Polygon, MultiPolygon,
box, GeometryCollection, Point, MultiPoint,
Expand Down Expand Up @@ -2524,7 +2524,9 @@ def triangulate_polygon(

def triangulate_polygon_s(
shape: Union[Polygon, MultiPolygon, gpd.GeoDataFrame, gpd.GeoSeries],
min_int_ang=30
aux_pts:Union[np.array,Point,MultiPoint,
gpd.GeoDataFrame,gpd.GeoSeries]=None,
min_int_ang=30,
) -> None:
'''Triangulate the input shape smoothly
Expand All @@ -2541,8 +2543,20 @@ def triangulate_polygon_s(
Generated triangulation
'''
#First triangulation. Smooth but adds extra nodes at boundary
if aux_pts is not None:
if isinstance(aux_pts, (np.ndarray, list, tuple)):
aux_pts = MultiPoint(aux_pts)
if isinstance(aux_pts, (Point, MultiPoint)):
aux_pts = gpd.GeoSeries(aux_pts)
elif isinstance(aux_pts, gpd.GeoDataFrame):
aux_pts = aux_pts.geometry
elif not isinstance(aux_pts, gpd.GeoSeries):
raise ValueError("Wrong input type_t for <aux_pts>!")

aux_pts = aux_pts.get_coordinates().values

min_int_ang = 'q'+str(min_int_ang)
mesh1 = triangulate_polygon(shape,opts=['p',min_int_ang])
mesh1 = triangulate_polygon(shape,opts=['p',min_int_ang],aux_pts=aux_pts)
#Find all boundary modes
nb = get_boundary_edges(mesh1)
nb = np.sort(list(set(nb.ravel())))
Expand Down Expand Up @@ -2720,7 +2734,8 @@ def create_mesh_from_mesh_diff(domain: Union[jigsaw_msh_t,
mesh_2: jigsaw_msh_t,
crs=CRS.from_epsg(4326),
min_int_ang=None,
buffer_domain = 0.001
buffer_domain = 0.001,
hfun_mesh = None
) -> jigsaw_msh_t:
'''
Create a triangulation for the area correspondent to
Expand Down Expand Up @@ -2763,35 +2778,64 @@ def create_mesh_from_mesh_diff(domain: Union[jigsaw_msh_t,
if isinstance(domain, (list)):
domain = pd.concat([gpd.GeoDataFrame\
(geometry=[get_mesh_polygons(i)\
.buffer(buffer_domain,join_style=2)],\
# .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.union_all()],crs=crs)
domain_buffer = gpd.GeoDataFrame(geometry=[i[-1].geometry.buffer(\
buffer_domain,join_style=2) for i in domain.iterrows()],crs=crs)
domain_buffer = domain_buffer.dissolve().explode(index_parts=True)
domain_buffer.crs = domain_buffer.estimate_utm_crs()
domain_buffer =domain_buffer.loc[domain_buffer['geometry'].area == \
max(domain_buffer['geometry'].area)]
domain_buffer.crs = crs
domain_buffer = gpd.GeoDataFrame(geometry=[domain_buffer.union_all()],
crs=crs)

mesh_1_poly = get_mesh_polygons(mesh_1)
mesh_2_poly = get_mesh_polygons(mesh_2)

poly_buffer = domain.union_all().difference(
poly_buffer = domain_buffer.union_all().difference(
gpd.GeoDataFrame(
geometry=[
get_mesh_polygons(mesh_1),
get_mesh_polygons(mesh_2),
mesh_1_poly,
mesh_2_poly,
],
crs = crs
).union_all()
)
gdf_full_buffer = gpd.GeoDataFrame(
geometry = [poly_buffer],crs=crs)

gap_poly = domain.union_all().\
difference(mesh_1_poly).\
difference(mesh_2_poly)

if hfun_mesh is not None:
hfun_mesh = clip_mesh_by_shape(hfun_mesh,gap_poly,fit_inside=True)
boundary = np.unique(get_boundary_edges(hfun_mesh))
all_nodes = hfun_mesh.vert2['coord']
hfun_nodes = np.delete(all_nodes, boundary, axis=0)
hfun_nodes = MultiPoint(hfun_nodes)
hfun_nodes = gpd.GeoDataFrame(geometry=
gpd.GeoSeries(intersection(
hfun_nodes,
gap_poly.buffer(-0.0001),
)))
if hfun_mesh is None:
hfun_nodes=None

if min_int_ang is None:
msht_buffer = triangulate_polygon(gdf_full_buffer)
msht_buffer = triangulate_polygon(gdf_full_buffer,
aux_pts=hfun_nodes)
else:
msht_buffer = triangulate_polygon_s(gdf_full_buffer,
min_int_ang=min_int_ang)
min_int_ang=min_int_ang,
aux_pts=hfun_nodes)
msht_buffer.crs = crs
msht_buffer=clip_mesh_by_shape(msht_buffer,gap_poly.buffer(buffer_domain))

return msht_buffer

Expand Down Expand Up @@ -2948,7 +2992,8 @@ def fix_small_el(mesh_w_problem: jigsaw_msh_t,
buffer_size=buffer_size)

fixed_mesh = merge_overlapping_meshes([mesh_w_problem.msh_t,patch],
adjacent_layers=adjacent_layers
adjacent_layers=adjacent_layers,
min_int_ang=30
)

#cleaning up mesh:
Expand All @@ -2969,6 +3014,7 @@ def merge_overlapping_meshes(all_msht: list,
buffer_size: float = 0.005,
buffer_domain: float = 0.01,
min_int_ang: int = 30,
hfun_mesh = None,
crs=CRS.from_epsg(4326)
) -> jigsaw_msh_t:
'''
Expand Down Expand Up @@ -3013,15 +3059,19 @@ def merge_overlapping_meshes(all_msht: list,
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]])

if hfun_mesh is None:
hfun_mesh = msht_combined

buff_mesh = create_mesh_from_mesh_diff([msht_combined,msht],
carved_mesh,msht,
min_int_ang=min_int_ang,
buffer_domain=buffer_domain,
hfun_mesh = hfun_mesh
)
msht_combined = merge_neighboring_meshes(buff_mesh,
carved_mesh,
msht
Expand All @@ -3031,7 +3081,10 @@ def merge_overlapping_meshes(all_msht: list,
)
del carved_mesh,buff_mesh,domain,msht

msht_combined = cleanup_folded_bound_el(msht_combined)
# msht_combined = cleanup_folded_bound_el(msht_combined)
cleanup_duplicates(msht_combined)
cleanup_isolates(msht_combined)
put_id_tags(msht_combined)

return msht_combined

Expand Down Expand Up @@ -3688,7 +3741,7 @@ def shptri_to_msht(triangulated_shp):
cleanup_duplicates(msht)
cleanup_isolates(msht)
put_id_tags(msht)

return msht


Expand Down
6 changes: 3 additions & 3 deletions tests/api/utils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -250,7 +250,7 @@ def test_create_mesh_from_mesh_diff(self):
patch.msh_t,
carved_mesh)

self.assertEqual(len(msht_buffer.tria3), 49)
self.assertEqual(len(msht_buffer.tria3), 43)

def test_merge_neighboring_meshes(self):
p0 = Path(__file__).parents[1] / "data" / "test_mesh_1.2dm"
Expand Down Expand Up @@ -278,7 +278,7 @@ def test_fix_small_el(self):

fixed_mesh = utils.fix_small_el(mesh,mesh_for_patch)

self.assertEqual(len(fixed_mesh.tria3), 1130876)
self.assertEqual(len(fixed_mesh.tria3), 1130233)

def test_merge_overlapping_meshes(self):
p = Path(__file__).parents[1] / "data" / "test_mesh_1.2dm"
Expand All @@ -288,7 +288,7 @@ def test_merge_overlapping_meshes(self):

smooth = utils.merge_overlapping_meshes([mesh.msh_t,patch.msh_t])

self.assertEqual(len(smooth.tria3), 1130935)
self.assertEqual(len(smooth.tria3), 1130295)


class FinalizeMesh(unittest.TestCase):
Expand Down

0 comments on commit a484d9e

Please sign in to comment.