1

Currently I am working in project that require from me extracting channel banks. I have been thinking about the best method that I can use to come up with good function that extract the location of channel banks. I don't want to use borrowed software such as Geo-ras or HEC-Ras. My main focus is to come up with a function which I can apply to many DEMs based on the the shape of the cross-section and extract the width and height from it. It is similar to the this picture. At the current stage, my thought is create multiple points at one point from the channel network and then create multiple consecutive points across the channel with the same interval and use their elevation to create slopes and their is a change in the slope (threshold that will create my first bank point until the it finds the second bank then it stops. I believe that what they call slope breaking method. From that I can extract the depth and width in the channel correctly.

Graph

5
  • What are your input data? A single DEM? Do you have a water surface raster?
    – Radar
    Jan 13, 2016 at 22:38
  • I am using DEM constructed from Lidar data
    – Mrmr
    Jan 14, 2016 at 0:44
  • What type of LiDAR? Green LiDAR would give you bathymetry (penetrating the water surface) - others may be giving you something close to the water surface.
    – Radar
    Jan 14, 2016 at 0:50
  • I am not interested in the water surface. I am more focused on the corss sections. Imagine that we have location each cell and elevation and when we draw it, it will give something similar to the schematic above. the challenge is to find the banks of that dummy schematic.
    – Mrmr
    Jan 14, 2016 at 0:55
  • You can try the River Bathymetry Toolkit mentioned in @GBG's answer below. It will do much of what you state above (create a centerline, create cross sections at intervals, extract elevations). Defining the river banks can be much trickier though depending on how variable your water level is.
    – Radar
    Jan 14, 2016 at 17:59

1 Answer 1

1

Here is some code that will take a line, and generate a perpendicular line to that line. You could take a river center line, divide it by vertex, run this tool, proportion the perpendicular lines, populate the sections with z values, and MAYBE glean some river banks from the results assuming that the surface model was collected when the river was at low flow and the banks are steep... But I have found there are no magic bullets for this type of work. You will need to do a lot of manual editing. A link to the tool is here and includes a readme: ftp://lnnr.lummi-nsn.gov/GIS_Scripts/createperpendicularlines.zip. Of course PolyGeo will scold me for just providing a potentially unmaintained link ;) so I have include the Python code below.

Also search the web for ESSA River Bathymetry Toolkit.

import sys, string, os, arcpy, traceback, math
def AddAndSubtractRadians(theta):
    return (theta + 1.57079632679, theta - 1.57079632679)       
def Distance(x1, y1, x2, y2):
    '''Cartesian distance formula'''
    return float(math.pow(((math.pow((x2-x1),2)) + (math.pow((y2 - y1),2))),.5))   
def CartesianToPolar(xy1, xy2):
    '''Given coordinate pairs as two lists or tuples, return the polar
    coordinates with theta in radians. Values are in true radians along the
    unit-circle, for example, 3.14 and not -0 like a regular python
    return.'''
    try:
        x1, y1, x2, y2 = float(xy1[0]), float(xy1[1]), float(xy2[0]), float(xy2[1])
        xdistance, ydistance = x2 - x1, y2 - y1
        distance = math.pow(((math.pow((x2 - x1),2)) + (math.pow((y2 - y1),2))),.5)
        if xdistance == 0:
            if y2 > y1:
                theta = math.pi/2
            else:
                theta = (3*math.pi)/2
        elif ydistance == 0:
            if x2 > x1:
                theta = 0
            else:
                theta = math.pi
        else:
            theta = math.atan(ydistance/xdistance)
            if xdistance > 0 and ydistance < 0:
                theta = 2*math.pi + theta
            if xdistance < 0 and ydistance > 0:
                theta = math.pi + theta
            if xdistance < 0 and ydistance < 0:
                theta = math.pi + theta
        return [distance, theta]
    except:
        print"Error in CartesianToPolar()"  
def PolarToCartesian(polarcoords):
    '''A tuple, or list, of polar values(distance, theta in radians) are
    converted to cartesian coords'''
    r = polarcoords[0]
    theta = polarcoords[1]
    x = r * math.cos(theta)
    y = r * math.sin(theta)
    return [x, y]
def gPrint(string):
    #print string
    arcpy.AddMessage(string) 
def GetMidpointOfAStraightLine( x1, y1, x2, y2):
    midx = (x1 + x2)/2   
    midy = (y1 + y2)/2
    return [midx, midy]
def GetSegmentLengths(thisrecordsgeometry):
    lengthlist = []
    positiontracker = 0
    totaldistance = 0
    try:
        for item in thisrecordsgeometry:
            x1 = thisrecordsgeometry[positiontracker][0]
            y1 = thisrecordsgeometry[positiontracker][1]
            x2 = thisrecordsgeometry[positiontracker +1][0]
            y2 = thisrecordsgeometry[positiontracker +1][1]
            segmentdistance = Distance(x1, y1, x2, y2)
            totaldistance = totaldistance + segmentdistance
            lengthlist.append(totaldistance)
            positiontracker += 1  
    except:
        pass
    return lengthlist
def stringToBoolean(x):
    if x == 'true':
         x = True
    else:
         x = False
    return x
gPrint("\n      Runing...Create Perpendicular Lines At The Midpoint of Line Features")
gPrint("      A Two-Bit Algorithms product.\n      Copyright 2011 Gerry Gabrisch ([email protected])\n")
infc = arcpy.GetParameterAsText(0)
distance = arcpy.GetParameterAsText(1)
fcname = arcpy.GetParameterAsText(2)
leftside = arcpy.GetParameterAsText(3)
rightside = arcpy.GetParameterAsText(4)
bothsides = arcpy.GetParameterAsText(5)

leftside = stringToBoolean(leftside)
rightside = stringToBoolean(rightside)
bothsides = stringToBoolean(bothsides)

gPrint("      Cracking Features...")
#Get the input line files geometry as a python list.
desc = arcpy.Describe(infc)
shapefieldname = desc.ShapeFieldName
rows = arcpy.SearchCursor(infc)
listofpointgeometry = []
gPrint("      Finding Midpoints...")
for row in rows:
    feat = row.getValue(shapefieldname)
    midpointDistance = (feat.length/2)
    partnum = 0
    partcount = feat.partCount
    thisrecordsgeometry = []
    while partnum < partcount:
        part = feat.getPart(partnum)
        pnt = part.next()
        pntcount = 0
        while pnt:
            thetuple = [pnt.X, pnt.Y]
            thisrecordsgeometry.append(thetuple)
            pnt = part.next()
            pntcount += 1
        partnum += 1
    startnode = [thisrecordsgeometry[0][0], thisrecordsgeometry[0][1]]
    endnode = [thisrecordsgeometry[-1][0], thisrecordsgeometry[-1][1]]
    #If the line is straight, then do this simple calculation to get its midpoint....
    if len(thisrecordsgeometry) == 2:
        midpoint = GetMidpointOfAStraightLine( thisrecordsgeometry[0][0], thisrecordsgeometry[0][1],thisrecordsgeometry[1][0], thisrecordsgeometry[1][1])
    #If the line is made of more than two verticies, sum up the line segments until the center point is located...    
    else:
        listofsegmentlengths = GetSegmentLengths(thisrecordsgeometry)
        nodecounter = 1
        for item in listofsegmentlengths:
            if item < midpointDistance:
                nodecounter += 1
            else:
                distanceback = item - midpointDistance 
                break
        polar = CartesianToPolar(thisrecordsgeometry[nodecounter], thisrecordsgeometry[nodecounter -1])            
        polardistancetomidpoint = [distanceback, polar[1]]
        midpointascartesiantemp = PolarToCartesian(polardistancetomidpoint)
        midx = thisrecordsgeometry[nodecounter][0] + midpointascartesiantemp[0]
        midy = thisrecordsgeometry[nodecounter] [1] + midpointascartesiantemp[1]
        midpoint = [midx,midy]
    listofpointgeometry.append([startnode,midpoint,endnode])  
#Create the data models to store the new geometery....    
featureList = []
array = arcpy.Array()
pnt = arcpy.Point()
gPrint("      Defining Perpendicular Lines...")
for pt in listofpointgeometry:
    startx = pt[0][0]
    starty = pt[0][1]
    midx = pt[1][0]
    midy = pt[1][1]
    endx = pt[2][0]
    endy = pt[2][1]
    #use the start point and end point to get a theta
    polarcoor = CartesianToPolar((startx,starty), (endx,endy))
    #Add and subtract the 90 degrees in radians from the line...
    ends = AddAndSubtractRadians(polarcoor[1])
    firstend = PolarToCartesian((float(distance),float(ends[0]))) 
    secondend = PolarToCartesian((float(distance),float(ends[1])))
    firstx2 = midx + firstend[0]
    firsty2 = midy + firstend[1]
    secondx2 = midx + secondend[0]
    secondy2 = midy + secondend[1]
    if bothsides:
        pnt.X, pnt.Y = firstx2 , firsty2
        array.add(pnt)
        pnt.X, pnt.Y = secondx2 , secondy2
        array.add(pnt)
        polyline = arcpy.Polyline(array)
        array.removeAll()
        featureList.append(polyline)
    if leftside:
        pnt.X, pnt.Y = midx , midy
        array.add(pnt)
        pnt.X, pnt.Y = firstx2 , firsty2
        array.add(pnt)
        polyline = arcpy.Polyline(array)
        array.removeAll()
        featureList.append(polyline)
    if rightside:
        pnt.X, pnt.Y = pnt.X, pnt.Y = midx , midy
        array.add(pnt)
        pnt.X, pnt.Y = secondx2 , secondy2
        array.add(pnt)
        polyline = arcpy.Polyline(array)
        array.removeAll()
        featureList.append(polyline)
gPrint("      Creating New Feature Class...")
arcpy.CopyFeatures_management(featureList, fcname)
spatialRef = arcpy.Describe(infc).spatialReference
arcpy.DefineProjection_management(fcname, spatialRef)   
gPrint ("      Done")      
4
  • As far as the River Bathymetry Toolkit goes, you can use it to "flood" an input DEM as a means of estimating the maximum bank height. It will also allow you to detrend the DEM to remove the influence of longitudinal slope (among many other things). Source: was a developer on the RBT.
    – Radar
    Jan 14, 2016 at 17:53
  • where I can modify the input files or directory from the code
    – Mrmr
    Jan 15, 2016 at 4:59
  • Thank you Nick O for your response it is useful but still the issue that I am working on now is how to extract from this created lines the properties of the cross section such as depth and width. Conceptually, I am thinking of extracting the elevation of multiple points and find the slope change. From the change I can create breaking points as my left and right banks.
    – Mrmr
    Jan 15, 2016 at 21:03
  • 1
    Thanks for pointing that out. I removed the error handler for clarity.
    – GBG
    Apr 24, 2017 at 14:46

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.