From 626bdf7b394201e846a9ab6f2c73c6e704f4113e Mon Sep 17 00:00:00 2001 From: felicio93 Date: Thu, 19 Dec 2024 13:28:19 -0500 Subject: [PATCH] added functions for cleaning RiverMapper invalid geoms --- ocsmesh/utils.py | 154 ++++++++++++++++++++++++++++++++++++++++++++- tests/api/utils.py | 4 +- 2 files changed, 153 insertions(+), 5 deletions(-) diff --git a/ocsmesh/utils.py b/ocsmesh/utils.py index 55f79d5..c6fe5ec 100644 --- a/ocsmesh/utils.py +++ b/ocsmesh/utils.py @@ -27,6 +27,7 @@ box, GeometryCollection, Point, MultiPoint, LineString, LinearRing) from shapely.ops import polygonize, linemerge, unary_union, triangulate +from shapely.validation import make_valid import geopandas as gpd import pandas as pd import utm @@ -3745,13 +3746,13 @@ def shptri_to_msht(triangulated_shp): return msht -def triangulate_rivermapper_poly(rm_poly): +def triangulate_poly(rm_poly): ''' - Creates triangulated mesh using the RiverMapperoutputs + Creates triangulated meshes from gpf with poly and multipoly Parameters ---------- - rm_poly : .shp (gpd) file with the RiverMapper outputs + rm_poly : .shp (gpd) (e.g.file with the RiverMapper outputs) Returns ------- @@ -3767,3 +3768,150 @@ def triangulate_rivermapper_poly(rm_poly): rm_mesh = shptri_to_msht(rm_poly) return rm_mesh + +def validate_poly(gdf): + ''' + Goes over all polygons in a gpf and applied the make_valid func + + Parameters + ---------- + gdf : .shp (gpd) (e.g.file with the RiverMapper outputs) + + Returns + ------- + gdf_valid, gdf_invalid: valid and invalid polygon in the gpdf + + Notes + ----- + + ''' + gdf = [make_valid(p) for p in gdf['geometry'].to_list()] + gdf = gpd.GeoDataFrame(geometry = + gpd.GeoSeries(gdf), + crs=4326).explode() + gdf_valid = gdf[gdf.geometry.geom_type.isin(['Polygon', + 'MultiPolygon'] + )].drop_duplicates() + gdf_invalid = gdf[~gdf.geometry.geom_type.isin(['Polygon', + 'MultiPolygon'] + )].drop_duplicates() + + return gdf_valid, gdf_invalid + +def find_polyneighbors(target_gdf, ref_gdf): + ''' + Finds all polygons in target_gdf that + border polygons in ref_gdf + + Parameters + ---------- + target_gdf: .shp (gpd) + ref_gdf: .shp (gpd) + + Returns + ------- + sjoin: .shp (gpd) with all target_gdf that + border polygons in ref_gdf + Notes + ----- + ''' + intersects = gpd.sjoin( + target_gdf, ref_gdf,how='inner', predicate='intersects' + ) + crosses = gpd.sjoin( + target_gdf, ref_gdf, how='inner', predicate='crosses' + ) + touches = gpd.sjoin( + target_gdf, ref_gdf, how='inner', predicate='touches' + ) + overlaps = gpd.sjoin( + target_gdf, ref_gdf, how='inner', predicate='overlaps' + ) + sjoin = pd.concat([intersects,crosses,touches,overlaps]) + sjoin = sjoin.drop_duplicates() + + return sjoin + +def validate_RMmesh(RMmesh): + ''' + Takes a mesh triangulated from a polygon and + removes invalid elements (and their neighbors) + + Parameters + ---------- + RMmesh: triangualed msh_t + + Returns + ------- + jigsaw_msh_t + River Mesh + + Notes + ----- + ''' + + tri = RMmesh.tria3['index'] + xy = RMmesh.vert2['coord'] + poly=[] + for t in tri: + p=[xy[t][0], xy[t][1], xy[t][2], xy[t][0], + ] + poly.append(p) + + RMmesh_poly = gpd.GeoDataFrame(gpd.GeoSeries( + [Polygon(p) for p in poly]), + columns=['geometry'], + crs=4326) + RMmesh_invalid = RMmesh_poly.loc[RMmesh_poly['geometry'].area < 1e-15] + + polyneighbors = find_polyneighbors(RMmesh_poly, RMmesh_invalid) + + el2remove_idx = np.concatenate((RMmesh_invalid.index.to_numpy(), + polyneighbors.index.to_numpy())) + new_tri = np.delete(tri, el2remove_idx, axis=0) + + msht_valid = msht_from_numpy( + coordinates = xy, + triangles = new_tri, + ) + cleanup_duplicates(msht_valid) + cleanup_isolates(msht_valid) + put_id_tags(msht_valid) + + return msht_valid + +def triangulate_rivermapper_poly(rm_poly): + ''' + Creates triangulated mesh using the RiverMapperoutputs + 1) Finds slivers (gdf_invalid) on RM shapefile + 2) Finds all polygons next to slivers + 3) Removes all polygons next to slivers + 4) Triangulate the slivers-free RM shapefile + 5) Finds/Removes non-matching vertices and their neighbor el + + Parameters + ---------- + rm_poly: .shp (gpd) file with the RiverMapper outputs + + Returns + ------- + jigsaw_msh_t + River Mesh + + Notes + ----- + + ''' + + # Finds slivers (gdf_invalid) on RM shapefile: + gdf_valid, gdf_invalid = validate_poly(rm_poly) + # Finds all polygons next to slivers: + polyneighbors = find_polyneighbors(gdf_valid, gdf_invalid) + # Removes all polygons next to slivers: + fixed_rm = gdf_valid[~gdf_valid.index.isin(polyneighbors.index)] + # Triangulate the slivers-free RM shapefile: + RMmesh = triangulate_poly(fixed_rm) + # Finds/Removes non-matching vertices and their neighbor el: + validated_RMmesh = validate_RMmesh(RMmesh) + + return validated_RMmesh diff --git a/tests/api/utils.py b/tests/api/utils.py index ff35a60..fcb358b 100644 --- a/tests/api/utils.py +++ b/tests/api/utils.py @@ -116,10 +116,10 @@ def test_clean_concv(self): class RiverMapper(unittest.TestCase): - def test_triangulate_rivermapper_poly(self): + def test_triangulate_poly(self): p = Path(__file__).parents[1] / "data" / "rm_poly.shp" rm_poly = gpd.read_file(p) - trias = utils.triangulate_rivermapper_poly(rm_poly) + trias = utils.triangulate_poly(rm_poly) self.assertEqual(len(trias.tria3), 56)