Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Improve merge overlap meshes #174

Merged
merged 3 commits into from
Sep 12, 2024
Merged

Improve merge overlap meshes #174

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
a56b4db
changed merge_overlapping_meshes + tests
felicio93 Sep 12, 2024
6e4eb9f
changed merge_overlapping_meshes + tests 02
felicio93 Sep 12, 2024
fc7e694
changed merge_overlapping_meshes + tests 03
felicio93 Sep 12, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion ocsmesh/geom/collector.py
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Loading