Merge pull request #170 from felicio93/rivermapper_tri

Rivermapper tri

Author: felicio93

.github/workflows/functional_test.yml

@@ -59,7 +59,7 @@ jobs:
 
     - name: 'Upload conda environment'
       if: runner.os != 'Windows'
-      uses: actions/upload-artifact@v2
+      uses: actions/upload-artifact@v4
       with:
         name: conda-environment-${{ matrix.python-version }}
         path: environment_py${{ matrix.python-version }}.yml
@@ -68,7 +68,7 @@ jobs:
 
     - name: 'Upload test dem'
       if: runner.os != 'Windows'
-      uses: actions/upload-artifact@v2
+      uses: actions/upload-artifact@v4
       with:
         name: test-dem-${{ matrix.python-version }}
         path: /tmp/test_dem.tif
@@ -81,7 +81,7 @@ jobs:
 
     - name: 'Upload Geom build results'
       if: runner.os != 'Windows'
-      uses: actions/upload-artifact@v2
+      uses: actions/upload-artifact@v4
       with:
         name: geom-build-results-${{ matrix.python-version }}
         path: test_shape
@@ -94,7 +94,7 @@ jobs:
 
     - name: 'Upload Hfun build results'
       if: runner.os != 'Windows'
-      uses: actions/upload-artifact@v2
+      uses: actions/upload-artifact@v4
       with:
         name: hfun-build-results-${{ matrix.python-version }}
         path: test.2dm
@@ -107,7 +107,7 @@ jobs:
 
     - name: 'Upload remesh results'
       if: runner.os != 'Windows'
-      uses: actions/upload-artifact@v2
+      uses: actions/upload-artifact@v4
       with:
         name: remesh-bydem-results-${{ matrix.python-version }}
         path: remeshed.2dm

.github/workflows/functional_test_2.yml

@@ -60,7 +60,7 @@ jobs:
 
 
     - name: 'Upload test dem'
-      uses: actions/upload-artifact@v2
+      uses: actions/upload-artifact@v4
       with:
         name: test-dem-${{ matrix.python-version }}
         path: /tmp/test_dem.tif
@@ -71,7 +71,7 @@ jobs:
       run: source tests/cli/build_geom.sh
 
     - name: 'Upload Geom build results'
-      uses: actions/upload-artifact@v2
+      uses: actions/upload-artifact@v4
       with:
         name: geom-build-results-${{ matrix.python-version }}
         path: test_shape
@@ -82,7 +82,7 @@ jobs:
       run: source tests/cli/build_hfun.sh
 
     - name: 'Upload Hfun build results'
-      uses: actions/upload-artifact@v2
+      uses: actions/upload-artifact@v4
       with:
         name: hfun-build-results-${{ matrix.python-version }}
         path: test.2dm
@@ -93,7 +93,7 @@ jobs:
       run: source tests/cli/remesh_by_dem.sh
 
     - name: 'Upload remesh results'
-      uses: actions/upload-artifact@v2
+      uses: actions/upload-artifact@v4
       with:
         name: remesh-bydem-results-${{ matrix.python-version }}
         path: remeshed.2dm

ocsmesh/utils.py

@@ -21,11 +21,12 @@
     RectBivariateSpline, griddata)
 from scipy import sparse, constants
 from scipy.spatial import cKDTree
+from shapely import difference 
 from shapely.geometry import ( # type: ignore[import]
         Polygon, MultiPolygon,
         box, GeometryCollection, Point, MultiPoint,
         LineString, LinearRing)
-from shapely.ops import polygonize, linemerge, unary_union
+from shapely.ops import polygonize, linemerge, unary_union, triangulate
 import geopandas as gpd
 import pandas as pd
 import utm
@@ -3579,3 +3580,137 @@ def batched(iterable, n):
     iterator = iter(iterable)
     while batch := tuple(islice(iterator, n)):
         yield batch
+
+
+def delaunay_within(gdf):
+    '''
+    Creates  the initial delaunay triangules for 
+    a gpd composed of polygons (only).
+    Selects those delaunay triangules that fall within domain.
+
+    Parameters
+    ----------
+    gdf : gpd of polygons
+
+    Returns
+    -------
+    gdf : gpd of triangulated polygons
+    
+    Notes
+    -----
+
+    '''
+
+    tt=[]
+    for polygon in gdf['geometry']:
+        try:
+            tri = [triangle for triangle in \
+                   triangulate(polygon) if triangle.within(polygon)]
+            tt.append(MultiPolygon(tri))
+        except:
+            pass
+    shape_tri = gpd.GeoDataFrame(geometry=gpd.GeoSeries(tt)).set_crs(crs=4326)
+    shape_tri = shape_tri[~shape_tri.is_empty].dropna()
+    return shape_tri
+
+
+def triangulate_shp(gdf):
+    '''
+    Fills out the gaps left by the delaunay_within
+
+    Parameters
+    ----------
+    gdf : gpd of polygons
+
+    Returns
+    -------
+    gdf : gpd of triangulated polygons
+    
+    Notes
+    -----
+
+    '''
+
+    shape_tri = [delaunay_within(gdf)]
+    shape_diff = gpd.GeoDataFrame(geometry =
+                                  gpd.GeoSeries(
+                                      gdf.difference(shape_tri[0])))
+    shape_diff = shape_diff[~shape_diff.is_empty].dropna()#.explode()
+    shape_diff_len = len(shape_diff)
+    while shape_diff_len>0:
+        shape_diff_tri = delaunay_within(shape_diff)
+        shape_tri.append(shape_diff_tri)
+        shape_diff = gpd.GeoDataFrame(geometry=
+                                      gpd.GeoSeries(difference(
+                                          shape_diff.union_all(),
+                                          shape_diff_tri.union_all()
+                                      )))
+        shape_diff = shape_diff[~shape_diff.is_empty].dropna().explode()
+        if len(shape_diff) == shape_diff_len:
+            break
+        shape_diff_len = len(shape_diff)
+
+    shape_final = gpd.GeoDataFrame(pd.concat(
+                                    shape_tri, ignore_index=True)).explode()
+    shape_final = shape_final[shape_final.geometry.type == 'Polygon']
+    shape_final.reset_index(drop=True, inplace=True)
+
+    return shape_final
+
+
+def shptri_to_msht(triangulated_shp):
+    '''
+    Converts a triangulated shapefile to msh_t
+
+    Parameters
+    ----------
+    triangulated_shp : triangulated gpd
+
+    Returns
+    -------
+    jigsaw_msh_t
+    
+    Notes
+    -----
+
+    '''
+
+    coords = []
+    verts = []
+    for t in enumerate(triangulated_shp['geometry']):
+        if len(np.array(t[-1].boundary.xy).T) == 4:
+            coords.append(np.array(t[-1].boundary.xy).T[:-1])
+            verts.append(np.array([0,1,2])+3*t[0])
+
+    msht = msht_from_numpy(coordinates = np.vstack(coords),
+                           triangles = verts,
+                            )
+    cleanup_duplicates(msht)
+    cleanup_isolates(msht)
+    put_id_tags(msht)
+    
+    return msht
+
+
+def triangulate_rivermapper_poly(rm_poly):
+    '''
+    Creates triangulated mesh using the RiverMapperoutputs
+
+    Parameters
+    ----------
+    rm_poly : .shp (gpd) file with the RiverMapper outputs
+
+    Returns
+    -------
+    jigsaw_msh_t
+        River Mesh
+
+    Notes
+    -----
+
+    '''
+
+    rm_poly = triangulate_shp(rm_poly)
+    rm_mesh = shptri_to_msht(rm_poly)
+
+    return rm_mesh

tests/api/utils.py

@@ -116,6 +116,13 @@ def test_clean_concv(self):
 
 
 class RiverMapper(unittest.TestCase):
+    def test_triangulate_rivermapper_poly(self):
+        p = Path(__file__).parents[1] / "data" / "rm_poly.shp"
+        rm_poly = gpd.read_file(p)
+        trias =  utils.triangulate_rivermapper_poly(rm_poly)
+
+        self.assertEqual(len(trias.tria3), 56)
+
     def test_quadrangulate_rivermapper_arcs(self):
         p = Path(__file__).parents[1] / "data" / "diag_arcs_clipped.shp"
         arcs_shp = gpd.read_file(p)

tests/data/rm_poly.cpg

@@ -0,0 +1 @@
+UTF-8

tests/data/rm_poly.dbf

@@ -0,0 +1 @@
+GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

tests/data/rm_poly.prj

@@ -0,0 +1,27 @@
+<!DOCTYPE qgis PUBLIC 'http://mrcc.com/qgis.dtd' 'SYSTEM'>
+<qgis version="3.36.0-Maidenhead">
+  <identifier></identifier>
+  <parentidentifier></parentidentifier>
+  <language></language>
+  <type></type>
+  <title></title>
+  <abstract></abstract>
+  <links/>
+  <dates/>
+  <fees></fees>
+  <encoding></encoding>
+  <crs>
+    <spatialrefsys nativeFormat="Wkt">
+      <wkt></wkt>
+      <proj4>+proj=longlat +datum=WGS84 +no_defs</proj4>
+      <srsid>0</srsid>
+      <srid>0</srid>
+      <authid></authid>
+      <description></description>
+      <projectionacronym></projectionacronym>
+      <ellipsoidacronym></ellipsoidacronym>
+      <geographicflag>false</geographicflag>
+    </spatialrefsys>
+  </crs>
+  <extent/>
+</qgis>