ArcGIS 10: create perpendicular transects to stream at specified intervals
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))
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)
If anyone wants to do this with r, I found a script written by Barry Rowlingson on R-sig-geo that still works fairly well.
Two functions: one computes evenly spaced points along a set of xy coordinates, the other works out the end point of transects centered on those points:
# How to generate perpendicular transects along a line feature
evenspace <- function(xy, sep, start = 0, size) {
dx <- c(0, diff(xy[, 1]))
dy <- c(0, diff(xy[, 2]))
dseg <- sqrt(dx^2 + dy^2)
dtotal <- cumsum(dseg)
linelength <- sum(dseg)
pos <- seq(start, linelength, by = sep)
whichseg <- unlist(lapply(pos, function(x) {
sum(dtotal <= x)
}))
pos <- data.frame(
pos = pos, whichseg = whichseg,
x0 = xy[whichseg, 1],
y0 = xy[whichseg, 2],
dseg = dseg[whichseg + 1],
dtotal = dtotal[whichseg],
x1 = xy[whichseg + 1, 1],
y1 = xy[whichseg + 1, 2]
)
pos$further <- pos$pos - pos$dtotal
pos$f <- pos$further / pos$dseg
pos$x <- pos$x0 + pos$f * (pos$x1 - pos$x0)
pos$y <- pos$y0 + pos$f * (pos$y1 - pos$y0)
pos$theta <- atan2(pos$y0 - pos$y1, pos$x0 - pos$x1)
return(pos[, c("x", "y", "x0", "y0", "x1", "y1", "theta")])
}
transect <- function(tpts, tlen) {
tpts$thetaT <- tpts$theta + pi / 2
dx <- tlen * cos(tpts$thetaT)
dy <- tlen * sin(tpts$thetaT)
return(
data.frame(
x0 = tpts$x + dx,
y0 = tpts$y + dy,
x1 = tpts$x - dx,
y1 = tpts$y - dy
)
)
}
Example:
# 1. make a meandering stream:
stream <- data.frame(x = seq(0, 20, len = 50))
stream$y <- sin(stream$x / 3)
head(stream)
#> x y
#> 1 0.0000000 0.0000000
#> 2 0.4081633 0.1356351
#> 3 0.8163265 0.2687633
#> 4 1.2244898 0.3969241
#> 5 1.6326531 0.5177490
#> 6 2.0408163 0.6290046
# 2. plot it - note asp = 1 is essential:
plot(stream, type = "l", asp = 1, lwd = 2, col = 'blue')
# 3. compute evenly spaced (1.5 units) transect centres:
tspts <- evenspace(stream, 1.5)
head(tspts)
#> x y x0 y0 x1 y1 theta
#> 1 0.000000 0.0000000 0.000000 0.0000000 0.4081633 0.1356351 -2.820767
#> 2 1.428493 0.4573134 1.224490 0.3969241 1.6326531 0.5177490 -2.853791
#> 3 2.883469 0.8193843 2.857143 0.8147982 3.2653061 0.8859022 -2.969119
#> 4 4.372037 0.9916890 4.081633 0.9779783 4.4897959 0.9972486 -3.094415
#> 5 5.868659 0.9246324 5.714286 0.9447499 6.1224490 0.8915592 3.012006
#> 6 7.340358 0.6401380 6.938776 0.7370314 7.3469388 0.6385503 2.904839
points(tspts, pch = "X")
# 4. compute perpendicular transects of length 2 units:
tslines <- transect(tspts, 2)
head(tslines)
#> x0 y0 x1 y1
#> 1 0.6307003 -1.897951 -0.6307003 1.897951
#> 2 1.9961840 -1.460427 0.8608022 2.375053
#> 3 3.2267096 -1.150942 2.5402285 2.789711
#> 4 4.4663564 -1.006086 4.2777177 2.989464
#> 5 5.6102103 -1.058598 6.1271081 2.907863
#> 6 6.8712621 -1.304071 7.8094546 2.584347
segments(tslines$x0, tslines$y0, tslines$x1, tslines$y1, col = 'red')
Extra: convert to long format so it's easy to import into ArcGIS
library(dplyr)
library(tidyr)
tslines_long <- tslines %>%
tbl_df() %>%
mutate(id = row_number()) %>%
select(id, everything()) %>%
pivot_longer(c(x0:y1),
names_to = c(".value", "set"),
names_pattern = "(.)(.)"
)
tslines_long
#> # A tibble: 28 x 4
#> id set x y
#> <int> <chr> <dbl> <dbl>
#> 1 1 0 0.631 -1.90
#> 2 1 1 -0.631 1.90
#> 3 2 0 2.00 -1.46
#> 4 2 1 0.861 2.38
#> 5 3 0 3.23 -1.15
#> 6 3 1 2.54 2.79
#> 7 4 0 4.47 -1.01
#> 8 4 1 4.28 2.99
#> 9 5 0 5.61 -1.06
#> 10 5 1 6.13 2.91
#> # … with 18 more rows