0

I was trying to re-construct the topology in stream networks, but I found some loops in the raw datasets after the Feature To Line tool.

The red lines are background flowline, the blue/aqua lines are the flowline with loops.

Loop in stream networks

Upon implementing the algorithm, I found out that there are some remaining loop even though I have removed many loops in the first run. (As a result, it seems to require an iterative approach.)

Nested loop

Is there any tool or algorithm to automatically remove one segment when loop occurs?

My initial thought is to use the start and end point location. If two line features share the same start-end or reverse, then one of them can be considered as a loop segment. But apparently this is not enough because I also noticed these unresolved loops after the first run. Complex loop

I tried to avoid using Arcpy because I want to do it across platform. My first implementation is using gdal python api.

1
  • 3
    Assuming that you are using ArcPy to try and re-construct the topology, can you provide a code snippet that illustrates what you have tried and where you are stuck, please?
    – PolyGeo
    Jan 6, 2021 at 19:15

2 Answers 2

1

Your arrow is pointing to a red loop that is part of a larger blue loop, so it's not clear what you are indicating! Do you want to break the larger blue loop or the smaller red loop?

Your suggestion of identifying segments where they share same from and to node ID's is valid for only the simplest situation, a loop around a small island. Even that could fail if your network contains pseudo nodes and most river network tend to have these.

The loops in your example have tributary junctions so that test would fail.

RivEX has a variety of quality control tools you can run on a river network, suggest you explore the help file and get a better idea of the complexity of the problem you are facing. RivEX can identify bifurcating nodes but it ultimate relies on you making the call on what to delete out as in my experience any automated approach will basically get it wrong.

If you are willing to accept errors, i.e you drop the main channel instead of the side branch but maintain overall network topology then an approach I use with RivEX is to attribute your network with Shreve order and then drop zero order polylines. This creates a single threaded network but as I say it can (and will) follow the first channel it comes across which may be the one you don't want. This can have knock on affects if you are doing barrier analysis.

RivEX can also tag polylines as being part of a loop to aid in identification.

1
  • Indeed, to tackle more complex loops, we need to some advance graphic algorithms. For example, a loop can be formed by more than 2 segments, that is what occurs in my examples as well.
    – Chang
    Jan 7, 2021 at 18:25
0

I was able to fix the simple ones using the following python code. It only uses gdal and numpy apis.

    import sys, os
    import numpy as np
    from osgeo import gdal, osr,ogr, gdalconst
    def search_duplicate_pair(aDic, aPt_pair):        
        nPair = len(aDic)
        iFlag = 0 # not in the dictionary
        for i in np.arange(0, nPair):
            start = aDic[i][0]
            end = aDic[i][1]
            start_x = start[0]
            start_y = start[1]
            end_x = end[0]
            end_y = end[1]
            x1 = aPt_pair[0][0]
            y1 = aPt_pair[0][1]
            x2 = aPt_pair[1][0]
            y2 = aPt_pair[1][1]

            a = np.power(  (start_x - x1 ) ,2)  + np.power(  (start_y - y1 ) ,2) + np.      power(  (end_x - x2 ) ,2)  + np.power(  (end_y - y2 ) ,2)
            #reverse
            b = np.power(  (start_x - x2 ) ,2)  + np.power(  (start_y - y2 ) ,2) + np.      power(  (end_x - x1 ) ,2)  + np.power(  (end_y - y1 ) ,2)
            if a < 0.0001 or b < 0.0001:
                #we found one repeating segment
                iFlag = 1
                break
            else:
                pass
            
        if iFlag == 0:
            #add it into the dic
            aDic.append(aPt_pair)

        return iFlag, aDic

    def remove_stream_network_loops(sFilename_in, sFilename_out):        
        sDriverName = "ESRI Shapefile"
        pDriver = ogr.GetDriverByName( sDriverName )
        if pDriver is None:
            print ("%s pDriver not available.\n" % sDriverName)
        else:
            print  ("%s pDriver IS available.\n" % sDriverName)

        pDataSource = pDriver.Open(sFilename_in, 0) # 0 means      read-only. 1 means       writeable.
        # Check to see if shapefile is found.
        if pDataSource is None:
            print ('Could not open %s' % (sFilename_in))
        else:
            print ('Opened %s' % (sFilename_in))
            pLayer = pDataSource.GetLayer()
            pSpatailRef = pLayer.GetSpatialRef() 
            #output
            pSrs = pSpatailRef
            pDatasetOut = pDriver.CreateDataSource(sFilename_out)
            pLayerOut = pDatasetOut.CreateLayer('flowline', pSrs, ogr.wkbLineString)
            # Add one attribute
            pLayerOut.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))
            pLayerDefn = pLayerOut.GetLayerDefn()
            pFeatureOut = ogr.Feature(pLayerDefn)
            lFeatureCount = pLayer.GetFeatureCount()
            print ( "Number of features in %s: %d" %  (os.path.basename(sFilename_in),      lFeatureCount) )        
            pSpatailRef = pLayer.GetSpatialRef() 
            #aFeature=[]
            j = 0
            aDic =[]
            for pFeature in pLayer:            
                pGeometry = pFeature.GetGeometryRef()                    
                npt = pGeometry.GetPointCount()    
                aPt_pair = [ pGeometry.GetPoint(0), pGeometry.GetPoint(npt-1)]
                iFlag, aDic = search_duplicate_pair(aDic, aPt_pair)
                j = j +1
                if(iFlag ==1):
                    pass    
                else:
                    #save this geomery
                    pFeatureOut.SetGeometry(pGeometry)
                    pFeatureOut.SetField("id", j)
                    pLayerOut.CreateFeature(pFeatureOut)
                    pass
                #aFeature.append(pGeometry)
                pass                
            pDatasetOut = pLayerOut = pFeatureOut  = None      
        return

It is possible to remove complex ones by updating the search_duplicate_pair function by allowing multiple pairs comparisons.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.