Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I am trying to do what my title says. The goal is to input a line, segment it in an evenly distributed way by a specified distance, convert the vertices to points, then draw lines perpendicular to the segment at a specified distance. These lines will then be converted to point in a similar fashion and return the maximum elevation value both east and west of each point. This whole process is meant to extract samples of an approximation of bank height for a river. I have found tools and processes that do the line segmentation (even though annoyingly it seems one cannot access the nifty Editor toolbar processes in python), but I am having trouble finding a script that makes the perpendicular transect lines. Has anyone done anything like this before? Note: I have found the script here, but it has too many malfunctions and I gave up trying to fix it (http://forums.arcgis.com/threads/49206-Perpendicular-transects-at-regular-intervals). Thank you!

share|improve this question
2  
Starting with the DEM, extracting the focal max (with a circular neighborhood of the desired radius) along the stream is a simple, direct way to obtain the desired information at all locations along the stream. –  whuber Feb 25 '12 at 18:51
    
Thanks for the input, Dr. Huber. I have tried using the River Bathymetry Toolkit, which relies on the kernel density method and was developed via advice you gave on the esri forums. How will the results of this focal max technique differ from that of the kernel density one (relating to bank height)? –  ndimhypervol Feb 26 '12 at 19:09
    
I am unsure of the details of the kernel density method, but in general this will compute a weighted (focal) average of elevations, which will not be the full bank height you're looking for. One necessary thing I lack--which you may have--is an operational definition of "bank height" for a stream. Is there a standard one? –  whuber Feb 26 '12 at 19:45
    
Sorry, correction: the toolbox is called "Riparian Topography Toolbox" by Thomas Dilts (forums.esri.com/Thread.asp?c=93&f=995&t=252754&mc=0). His term was "height above river" and was basically producing a detrended DEM equal to height above the centerline of the river (flow accumulation product). This is then used to "flood" the DEM by specifying a threshold height and running cost distance on a binary raster. This technique is valid and runs very quickly, but for my dataset results in a raster product with heights above river greater than 0 on the river thalweg itself. –  ndimhypervol Feb 26 '12 at 20:24
1  
Re: local bank height. You could always go back to the beginning of that ESRI forums discussion, which took a Euclidean Allocation solution as its point of departure: this will be quite close to the bank height you're looking for. –  whuber Feb 26 '12 at 20:35

4 Answers 4

we had to do something similar and found you can manually do it. It is tedious, but by creating a random point on the line and then adding equal interval points from that origin, you have equal interval points. Then snap to each point and draw a perpendicular line from your stream line or thalweg line, which you have to do individually in each direction (right and left). You can then merge those perpendicular lines and use them as transects, or convert to points. There are two toolsets which may be useful, one is called ET GeoWizards which requires a license, but has tools which do pretty much exactly what you need. The other option is a freely available tool set, called "Geospatial Modelling Environment" (used to be Hawths' tools). See the vector tools "generate random points" and "convert lines to points". Hope this helps...it is a tedious task.

share|improve this answer
up vote 1 down vote accepted

I found code courtesy of Gerry Gabrisch posted on an old ArcGIS 9.2 forum (http://arcscripts.esri.com/details.asp?dbid=15756) and modified it so that it works in ArcGIS 10. The original output was a text file with coordinates, and this version writes geometry to an empty shapefile. Assume "inLine" is your input line feature. All that needs to be done once this code is run is to call the function "buildTransects". The code works pretty well, but all perpendicular lines in the output do not have the same length. Perhaps needs a bit more tweaking, but the lengths match up just enough for the least discriminatory user =D

# Create empty polyline for transects
spatialref = arcpy.Describe(inLine).spatialReference
if arcpy.Exists('transects'): 
    arcpy.Delete_management('transects')
transects = arcpy.CreateFeatureclass_management(temp,'transects.shp',"POLYLINE",'','','',spatialref)
arcpy.AddField_management(transects, "slope", "DOUBLE")

# inLine = input line, transects = new shapefile created above,
# transDist = length of transect (approximate), workspace = workspace
def buildTransects(inLine,transects,transDist,workspace):

    # Define orientation (start, mid, end) and draw transect function
    # This function is defined before it is called
    def orientTransect(feat,ix,iy):

            # If the line is horizontal or vertical, the slope and
            # negative reciprocal calculations will fail, so do this instead
            if starty==endy or startx==endx:
                if starty == endy:
                    y1 = iy + transDist
                    y2 = iy - transDist
                    x1 = ix
                    x2 = ix

                if startx == endx:
                    y1 = iy
                    y2 = iy
                    x1 = ix + transDist
                    x2 = ix - transDist


            else:

                # Get slope of line
                m = ((starty - endy)/(startx - endx))

                # Get negative reciprocal
                negativereciprocal = -1*((startx - endx)/(starty - endy))

                # For all values of slope, calculate perpendicular line
                # with length = transDist
                if m > 0:
                    if m >= 1:
                        y1 = negativereciprocal*(transDist)+ iy
                        y2 = negativereciprocal*(-transDist) + iy
                        x1 = ix + transDist
                        x2 = ix - transDist
                    if m < 1:
                        y1 = iy + transDist
                        y2 = iy - transDist
                        x1 = (transDist/negativereciprocal) + ix
                        x2 = (-transDist/negativereciprocal)+ ix

                if m < 0:
                    if m >= -1:
                        y1 = iy + transDist
                        y2 = iy - transDist
                        x1 = (transDist/negativereciprocal) + ix
                        x2 = (-transDist/negativereciprocal)+ ix

                    if m < -1:
                        y1 = negativereciprocal*(transDist)+ iy
                        y2 = negativereciprocal*(-transDist) + iy
                        x1 = ix + transDist
                        x2 = ix - transDist
            point1.X = x1
            point1.Y = y1
            point2.X = x2
            point2.Y = y2
            lineArray.add(point1)
            lineArray.add(point2)

            del x1
            del x2
            del y1
            del y2

    # Create search cursor in inLine
    rows = arcpy.SearchCursor(inLine)

    # Get number of records in inLine
    numRecords = int(arcpy.GetCount_management(inLine).getOutput(0))


    # Create new point files and array to collect values
    point1 = arcpy.Point()
    point2 = arcpy.Point()
    lineArray = arcpy.Array()

    # Define counter
    counter = 0

    # Loop over rows in outLine
    for row in rows:

        # Create the geometry object
        feat = row.Shape

        # Get coordinate values as lists
        firstpoint = feat.firstPoint
        lastpoint = feat.lastPoint
        midpoint = feat.centroid

        # Get X and Y values for each point
        startx = firstpoint.X
        starty = firstpoint.Y
        endx = lastpoint.X
        endy = lastpoint.Y
        midx = midpoint.X
        midy = midpoint.Y

        m = ((starty - endy)/(startx - endx))

        # For all points besides the last one
        if counter < numRecords - 1:
            orientTransect(feat,startx,starty)
        # For the last point
        else:
            orientTransect(feat,endx,endy)

        #Create insert cursor
        cur = arcpy.InsertCursor(transects)

        #Insert new row from array
        feat = cur.newRow()
        feat.slope = m
        feat.shape = lineArray
        cur.insertRow(feat)
        lineArray.removeAll()

        del cur

        # Increase counter by 1 and start again
        counter = counter + 1

    del row
    del rows

    printit('Added %s transects to inLine and saved to temp/transects.shp' % str(counter))
share|improve this answer

I did this awhile back. Here’s what I did. 1. Using the linear referenced river centerlines, added hatches every one-mile (you can choose the interval). Made the hatches 2 miles wide. (you can choose the width) 2. Converted hatches to graphic lines (draw toolbar) 3. Converted graphic lines into a line feature class. 4. Assigned route ID and measure back to the feature class using “Locate features along Routes”. 5. The resultant lines were perpendicular to the river line, not necessary the general floodplain. Rotated several cross sections to conform to the floodplain trends. This was done on more than one occasion. From here I was able to extract elevations from the DEM and then find min, max, etc.

share|improve this answer

Also see Perpendicular Transects ArcGIS 10 toolbox and python script by Mateus Ferreira. It creates perpendicular lines at regular intervals along the line, but also has an option for deriving the split distance from a field name.

The page links back to several related discussions with partial solutions here on GIS Stack Exchange which may be of interest to those building or enhancing their own tools.

Mateus' script:

# Import system modules
import arcpy
from arcpy import env
import math
arcpy.env.overwriteOutput = True

#Set environments
arcpy.env.overwriteOutput = True
arcpy.env.XYResolution = "0.00001 Meters"
arcpy.env.XYTolerance = "0.0001 Meters"

# Set local variables
env.workspace = arcpy.GetParameterAsText(0)
Lines=arcpy.GetParameterAsText(1)
SplitType=arcpy.GetParameterAsText(2)
DistanceSplit=float(arcpy.GetParameterAsText(3))
TransecLength=arcpy.GetParameterAsText(4)
TransecLength_Unit=arcpy.GetParameterAsText(5)
OutputTransect=arcpy.GetParameterAsText(6)

# Def splitline module
###START SPLIT LINE CODE IN A SAME DISTANCE### Source: http://nodedangles.wordpress.com/2011/05/01/quick-dirty-arcpy-batch-splitting-polylines-to-a-specific-length/
def splitline (inFC,FCName,alongDist):

    OutDir = env.workspace
    outFCName = FCName
    outFC = OutDir+"/"+outFCName

    def distPoint(p1, p2):
        calc1 = p1.X - p2.X
        calc2 = p1.Y - p2.Y

        return math.sqrt((calc1**2)+(calc2**2))

    def midpoint(prevpoint,nextpoint,targetDist,totalDist):
        newX = prevpoint.X + ((nextpoint.X - prevpoint.X) * (targetDist/totalDist))
        newY = prevpoint.Y + ((nextpoint.Y - prevpoint.Y) * (targetDist/totalDist))
        return arcpy.Point(newX, newY)

    def splitShape(feat,splitDist):
        # Count the number of points in the current multipart feature
        #
        partcount = feat.partCount
        partnum = 0
        # Enter while loop for each part in the feature (if a singlepart feature
        # this will occur only once)
        #
        lineArray = arcpy.Array()

        while partnum < partcount:
              # Print the part number
              #
              #print "Part " + str(partnum) + ":"
              part = feat.getPart(partnum)
              #print part.count

              totalDist = 0

              pnt = part.next()
              pntcount = 0

              prevpoint = None
              shapelist = []

              # Enter while loop for each vertex
              #
              while pnt:

                    if not (prevpoint is None):
                        thisDist = distPoint(prevpoint,pnt)
                        maxAdditionalDist = splitDist - totalDist

                        print thisDist, totalDist, maxAdditionalDist

                        if (totalDist+thisDist)> splitDist:
                              while(totalDist+thisDist) > splitDist:
                                    maxAdditionalDist = splitDist - totalDist
                                    #print thisDist, totalDist, maxAdditionalDist
                                    newpoint = midpoint(prevpoint,pnt,maxAdditionalDist,thisDist)
                                    lineArray.add(newpoint)
                                    shapelist.append(lineArray)

                                    lineArray = arcpy.Array()
                                    lineArray.add(newpoint)
                                    prevpoint = newpoint
                                    thisDist = distPoint(prevpoint,pnt)
                                    totalDist = 0

                              lineArray.add(pnt)
                              totalDist+=thisDist
                        else:
                              totalDist+=thisDist
                              lineArray.add(pnt)
                              #shapelist.append(lineArray)
                    else:
                        lineArray.add(pnt)
                        totalDist = 0

                    prevpoint = pnt                
                    pntcount += 1

                    pnt = part.next()

                    # If pnt is null, either the part is finished or there is an
                    #   interior ring
                    #
                    if not pnt:
                        pnt = part.next()
                        if pnt:
                              print "Interior Ring:"
              partnum += 1

        if (lineArray.count > 1):
              shapelist.append(lineArray)

        return shapelist

    if arcpy.Exists(outFC):
        arcpy.Delete_management(outFC)

    arcpy.Copy_management(inFC,outFC)

    #origDesc = arcpy.Describe(inFC)
    #sR = origDesc.spatialReference

    #revDesc = arcpy.Describe(outFC)
    #revDesc.ShapeFieldName

    deleterows = arcpy.UpdateCursor(outFC)
    for iDRow in deleterows:       
         deleterows.deleteRow(iDRow)

    try:
        del iDRow
        del deleterows
    except:
        pass

    inputRows = arcpy.SearchCursor(inFC)
    outputRows = arcpy.InsertCursor(outFC)
    fields = arcpy.ListFields(inFC)

    numRecords = int(arcpy.GetCount_management(inFC).getOutput(0))
    OnePercentThreshold = numRecords // 100

    #printit(numRecords)

    iCounter = 0
    iCounter2 = 0

    for iInRow in inputRows:
        inGeom = iInRow.shape
        iCounter+=1
        iCounter2+=1    
        if (iCounter2 > (OnePercentThreshold+0)):
              #printit("Processing Record "+str(iCounter) + " of "+ str(numRecords))
              iCounter2=0

        if (inGeom.length > alongDist):
              shapeList = splitShape(iInRow.shape,alongDist)

              for itmp in shapeList:
                    newRow = outputRows.newRow()
                    for ifield in fields:
                        if (ifield.editable):
                              newRow.setValue(ifield.name,iInRow.getValue(ifield.name))
                    newRow.shape = itmp
                    outputRows.insertRow(newRow)
        else:
              outputRows.insertRow(iInRow)

    del inputRows
    del outputRows

    #printit("Done!")
###END SPLIT LINE CODE IN A SAME DISTANCE###

# Create "General" file geodatabase
WorkFolder=env.workspace
General_GDB=WorkFolder+"\General.gdb"
arcpy.CreateFileGDB_management(WorkFolder, "General", "CURRENT")
env.workspace=General_GDB

#Unsplit Line
LineDissolve="LineDissolve"
arcpy.Dissolve_management (Lines, LineDissolve,"", "", "SINGLE_PART")
LineSplit="LineSplit"

#Split Line
if SplitType=="Split at approximate distance":
    splitline(LineDissolve, LineSplit, DistanceSplit)
else:
    arcpy.SplitLine_management (LineDissolve, LineSplit)

#Add fields to LineSplit
FieldsNames=["LineID", "Direction", "Azimuth", "X_mid", "Y_mid", "AziLine_1", "AziLine_2", "Distance"]
for fn in FieldsNames:
    arcpy.AddField_management (LineSplit, fn, "DOUBLE")

#Calculate Fields
CodeBlock_Direction="""def GetAzimuthPolyline(shape):
 radian = math.atan((shape.lastpoint.x - shape.firstpoint.x)/(shape.lastpoint.y - shape.firstpoint.y))
 degrees = radian * 180 / math.pi
 return degrees"""

CodeBlock_Azimuth="""def Azimuth(direction):
 if direction < 0:
  azimuth = direction + 360
  return azimuth
 else:
  return direction"""
CodeBlock_NULLS="""def findNulls(fieldValue):
    if fieldValue is None:
        return 0
    elif fieldValue is not None:
        return fieldValue"""
arcpy.CalculateField_management (LineSplit, "LineID", "!OBJECTID!", "PYTHON_9.3")
arcpy.CalculateField_management (LineSplit, "Direction", "GetAzimuthPolyline(!Shape!)", "PYTHON_9.3", CodeBlock_Direction)
arcpy.CalculateField_management (LineSplit, "Direction", "findNulls(!Direction!)", "PYTHON_9.3", CodeBlock_NULLS)
arcpy.CalculateField_management (LineSplit, "Azimuth", "Azimuth(!Direction!)", "PYTHON_9.3", CodeBlock_Azimuth)
arcpy.CalculateField_management (LineSplit, "X_mid", "!Shape!.positionAlongLine(0.5,True).firstPoint.X", "PYTHON_9.3")
arcpy.CalculateField_management (LineSplit, "Y_mid", "!Shape!.positionAlongLine(0.5,True).firstPoint.Y", "PYTHON_9.3")
CodeBlock_AziLine1="""def Azline1(azimuth):
 az1 = azimuth + 90
 if az1 > 360:
  az1-=360
  return az1
 else:
  return az1"""
CodeBlock_AziLine2="""def Azline2(azimuth):
 az2 = azimuth - 90
 if az2 < 0:
  az2+=360
  return az2
 else:
  return az2"""
arcpy.CalculateField_management (LineSplit, "AziLine_1", "Azline1(!Azimuth!)", "PYTHON_9.3", CodeBlock_AziLine1)
arcpy.CalculateField_management (LineSplit, "AziLine_2", "Azline2(!Azimuth!)", "PYTHON_9.3", CodeBlock_AziLine2) 
arcpy.CalculateField_management (LineSplit, "Distance", TransecLength, "PYTHON_9.3")

#Generate Azline1 and Azline2
spatial_reference=arcpy.Describe(Lines).spatialReference
Azline1="Azline1"
Azline2="Azline2"
arcpy.BearingDistanceToLine_management (LineSplit, Azline1, "X_mid", "Y_mid", "Distance", TransecLength_Unit, "AziLine_1", "DEGREES", "GEODESIC", "LineID", spatial_reference)
arcpy.BearingDistanceToLine_management (LineSplit, Azline2, "X_mid", "Y_mid", "Distance", TransecLength_Unit, "AziLine_2", "DEGREES", "GEODESIC", "LineID", spatial_reference)

#Create Azline and append Azline1 and Azline2
Azline="Azline"
arcpy.CreateFeatureclass_management(env.workspace, "Azline", "POLYLINE", "", "", "", spatial_reference)
arcpy.AddField_management (Azline, "LineID", "DOUBLE")
arcpy.Append_management ([Azline1, Azline2], Azline, "NO_TEST")

#Dissolve Azline
Azline_Dissolve="Azline_Dissolve"
arcpy.Dissolve_management (Azline, Azline_Dissolve,"LineID", "", "SINGLE_PART")

#Add Fields to Azline_Dissolve
FieldsNames2=["x_start", "y_start", "x_end", "y_end"]
for fn2 in FieldsNames2:
    arcpy.AddField_management (Azline_Dissolve, fn2, "DOUBLE")

#Calculate Azline_Dissolve fields
arcpy.CalculateField_management (Azline_Dissolve, "x_start", "!Shape!.positionAlongLine(0,True).firstPoint.X", "PYTHON_9.3") 
arcpy.CalculateField_management (Azline_Dissolve, "y_start", "!Shape!.positionAlongLine(0,True).firstPoint.Y", "PYTHON_9.3")
arcpy.CalculateField_management (Azline_Dissolve, "x_end", "!Shape!.positionAlongLine(1,True).firstPoint.X", "PYTHON_9.3")
arcpy.CalculateField_management (Azline_Dissolve, "y_end", "!Shape!.positionAlongLine(1,True).firstPoint.Y", "PYTHON_9.3")

#Generate output file
arcpy.XYToLine_management (Azline_Dissolve, OutputTransect,"x_start", "y_start", "x_end","y_end", "", "", spatial_reference)

#Delete General.gdb
arcpy.Delete_management(General_GDB)
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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