add the fix to resolve the island flowline issue

pyflowline/algorithms/auxiliary/find_index_in_list.py

@@ -83,7 +83,6 @@ def find_vertex_on_edge(aVertex_in, pEdge_in):
 
     return iFlag_exist, npoint , aIndex_order
 
-
 def find_edge_in_list(aEdge_in, pEdge_in):
     """[find the index of an edge in a list]
 
@@ -109,7 +108,6 @@ def find_edge_in_list(aEdge_in, pEdge_in):
 
     return iFlag_exist, lIndex
 
-
 def find_flowline_in_list(aFlowline_in, pFlowline_in):
     """[find the index of a flowline in a list]
 
@@ -170,7 +168,6 @@ def check_if_duplicates(aList_in):
     #return iFlag_unique
     return int(len(aList_in) == len(set(aList_in)))
 
-
 def add_unique_vertex(aVertex_in, pVertex_in, dThreshold_in = 1.0E-6):
     """[add a vertex to a list if it is not already included]
 

pyflowline/algorithms/index/define_stream_order.py

@@ -69,7 +69,6 @@ def define_stream_order(aFlowline_in, aConfluence_in):
                 pBound = (x - 1E-5, y - 1E-5, x + 1E-5, y + 1E-5)
                 index_confluence.insert(i, pBound)
 
-
             while aFlowline_in[0].iStream_order < 0:
                 for i in range(nConfleunce):
                     if aFlag_confluence_treated[i] == 1:

pyflowline/algorithms/intersect/intersect_flowline_with_mesh.py

@@ -13,8 +13,10 @@
 if iFlag_cython is not None:
     from tinyr import RTree
     iFlag_use_rtree = 1
+    from pyflowline.algorithms.cython.kernel import find_vertex_in_list
 else:
-    iFlag_use_rtree =0
+    iFlag_use_rtree = 0
+    from pyflowline.algorithms.auxiliary.find_vertex_in_list import find_vertex_in_list
     pass
 
 def intersect_flowline_with_mesh(iMesh_type_in, sFilename_mesh_in, sFilename_flowline_in, sFilename_output_in):
@@ -29,7 +31,7 @@ def intersect_flowline_with_mesh(iMesh_type_in, sFilename_mesh_in, sFilename_flo
 
 
     pDriver_geojson = ogr.GetDriverByName( "GeoJSON")
-    aCell=list()
+    #aCell=list()
     aCell_intersect=list()
 
     pDataset_mesh = pDriver_geojson.Open(sFilename_mesh_in, 0)
@@ -179,10 +181,10 @@ def intersect_flowline_with_mesh(iMesh_type_in, sFilename_mesh_in, sFilename_flo
 
                     #replace flowline length if there is an actual flowline
                     aCell_intersect.append(pCell)
-                    aCell.append(pCell)
+                    #aCell.append(pCell)
                 else:
                     pCell.iFlag_intersected = 0
-                    aCell.append(pCell)
+                    #aCell.append(pCell)
                     pass
     else:
 
@@ -292,14 +294,84 @@ def intersect_flowline_with_mesh(iMesh_type_in, sFilename_mesh_in, sFilename_flo
 
                     #replace flowline length if there is an actual flowline
                     aCell_intersect.append(pCell)
-                    aCell.append(pCell)
+                    #aCell.append(pCell)
                 else:
                     pCell.iFlag_intersected = 0
-                    aCell.append(pCell)
+                    #aCell.append(pCell)
                     pass
 
             else:
                 pass
 
 
-    return  aCell, aCell_intersect, aFlowline_intersect_all
+
+    #quality control
+    #find the flowline that has the largest stream segment index
+    pVertex_outlet = None
+    iSegment_max = 0
+    iIndex_flowline = -1
+    for i in range(nfeature_flowline):
+        pFeature_flowline = pLayer_flowline.GetFeature(i)
+        iStream_segment = pFeature_flowline.GetField("stream_segment")
+        if (iStream_segment > iSegment_max):
+            iSegment_max = iStream_segment
+            iIndex_flowline = i
+
+    if (iIndex_flowline != -1):
+        #find all the flowlines and cells that have the same stream segment
+
+        #get the first vertex of the flowline
+        pFeature_flowline = pLayer_flowline.GetFeature(iIndex_flowline)
+        pGeometry_flowline = pFeature_flowline.GetGeometryRef()
+        aCoords = list()
+        for i in range(0, pGeometry_flowline.GetPointCount()):
+            pt = pGeometry_flowline.GetPoint(i)
+            aCoords.append( [ pt[0], pt[1]])
+        dummy1= np.array(aCoords)
+        pFlowline = convert_gcs_coordinates_to_flowline(dummy1)
+        pVertex_start = pFlowline.pVertex_start
+        pVertex_outlet = pFlowline.pVertex_end
+        aFlowline_last = list()
+        nFlowline_intersect = len(aFlowline_intersect_all)
+        for i in range(nFlowline_intersect):
+            pFlowline = aFlowline_intersect_all[i]
+            if ( pFlowline.iStream_segment == iSegment_max):
+                aFlowline_last.append(pFlowline)
+
+
+        aFlowline_stay = list()
+        #now check the connectivity of the flowlines because the flowlines are not necessarily connected
+        iFlag_done = 0
+        while (iFlag_done ==0):
+            nFlowline_last = len(aFlowline_last)
+            for i in range(nFlowline_last):
+                pFlowline = aFlowline_last[i]
+                pVertex_start_dummy = pFlowline.pVertex_start
+                pVertex_end_dummy = pFlowline.pVertex_end
+                if pVertex_start_dummy == pVertex_start:
+                    #found the head
+                    iFlag_found = 1
+                    pVertex_start = pVertex_end_dummy
+                    aFlowline_stay.append(i)
+                    pVertex_outlet = pVertex_end_dummy
+                    break
+
+            if iFlag_found == 1:
+                iFlag_done = 0
+            else:
+                iFlag_done = 1
+
+        #check the size of stay
+        #nFlowline_stay = len(aFlowline_stay)
+        #if (nFlowline_stay == nFlowline_intersect):
+        #    print('All flowlines in the last segment are connected')
+        #else:
+        #    for i in range(nFlowline_last):
+        #        if i not in aFlowline_stay:
+        #            pFlowline = aFlowline_last[i]
+        #            #remove this flowline from the result
+        #            pass
+
+
+    #return  aCell, aCell_intersect, aFlowline_intersect_all, pVertex_outlet
+    return   aCell_intersect, aFlowline_intersect_all, pVertex_outlet

pyflowline/classes/basin.py

@@ -554,7 +554,7 @@ def basin_flowline_simplification(self):
             ptimer.start()
             nFlowline_before = len(aFlowline_basin_simplified)
             #this flowline is ordered from downstream to upstream
-            aFlowline_basin_simplified = correct_flowline_direction(aFlowline_basin_simplified,  pVertex_outlet )
+            aFlowline_basin_simplified = correct_flowline_direction(aFlowline_basin_simplified,  pVertex_outlet)
             nFlowline_after = len(aFlowline_basin_simplified)
             ptimer.stop()
             #update outlet
@@ -736,7 +736,7 @@ def basin_reconstruct_topological_relationship(self, iMesh_type, sFilename_mesh)
             print('Basin ',  self.sBasinID, 'Start flowline and mesh intersection')
             sys.stdout.flush()
             ptimer.start()
-            aCell, aCell_intersect_basin, aFlowline_intersect_all = intersect_flowline_with_mesh(iMesh_type, sFilename_mesh, \
+            aCell_intersect_basin, aFlowline_intersect_all, pVertex_outlet_checked = intersect_flowline_with_mesh(iMesh_type, sFilename_mesh, \
                 sFilename_flowline_in, sFilename_flowline_intersect_out)
             ptimer.stop()
 
@@ -748,13 +748,16 @@ def basin_reconstruct_topological_relationship(self, iMesh_type, sFilename_mesh)
             print('Error in flowline and mesh intersection.')
 
         #not an ideal setup, but it could be improved
-        if self.iFlag_simplification_done ==1:
-            pVertex_outlet_initial = self.pVertex_outlet
-        else:
-            point= dict()
-            point['dLongitude_degree'] = self.dLongitude_outlet_degree
-            point['dLatitude_degree'] = self.dLatitude_outlet_degree
-            pVertex_outlet_initial=pyvertex(point)
+        #if self.iFlag_simplification_done == 1:
+        pVertex_outlet_initial = pVertex_outlet_checked #self.pVertex_outlet
+        #else:
+            #point= dict()
+            #point['dLongitude_degree'] = self.dLongitude_outlet_degree
+            #point['dLatitude_degree'] = self.dLatitude_outlet_degree
+            #pVertex_outlet_initial= pVertex_outlet_checked #pyvertex(point)
+
+        print('Outlet ID initial', pVertex_outlet_initial.lID)
+        print('Outlet initial location', pVertex_outlet_initial.dLongitude_degree, pVertex_outlet_initial.dLatitude_degree)
 
         #from this point, aFlowline_basin is conceptual
         #segment based