Spatial Thoughts on LinkedIn posted the following challenge…Can you take a dataset of a river centerline and varying lengths of line transects and make them exactly the same length? using any tool or platform of choice. This post uses ArcPy with ArcGIS Pro and can be used with a Basic license. Here is a link to the posting in LinkedIn.
Interested in learning ArcPy? Check out this course!
we need to import two modules arcpy and math.
import arcpy
import math
Our tool will require three user inputs; transects (linear Feature Layer), centerline (linear Feature Layer), and line_length (Double).
## input transect linear features
transects = arcpy.GetParameter(0)
## input centerline
centerline = arcpy.GetParameter(1)
## the length to make each transect
line_length = arcpy.GetParameter(2)
We need two function to help us with calculations. The first is to get the angle between two points…
## http://wikicode.wikidot.com/get-angle-of-line-between-two-points
def getAngle(pt1, pt2):
"""
Get the angle between two points
Args:
pt1 an arcpy Point object
pt2 an arcpy Point object
Returns:
A float representing the angle between pt1 and pt2
"""
## pt2 X value - pt1 X value
x_diff = pt2.X - pt1.X
## pt2 Y value - pt1 Y value
y_diff = pt2.Y - pt1.Y
return math.degrees(math.atan2(y_diff, x_diff))
… and the second is to get the start and end points of the new line to create based on the new line length.
## get a point based on bearing and distance from another point
## in our case the other point is the intersection of the transect and centerline
def getLineEnds(pt, bearing, dist):
"""
Get the start/end point of the new transect line
Args:
pt1 an arpy Point object
bearing float representing an angle between two points on the line
dist distance along the bearing where we want to locate the point
Returns:
A tuple representing x,y coords of start/end of new transect line
"""
## convert to radians
bearing = math.radians(bearing)
## some mathsy stuff that I dont have a clue about
x = pt.X + dist * math.cos(bearing)
y = pt.Y + dist * math.sin(bearing)
return (x, y)
And we have a couple of required objects to aid with the workflow. We get half of the line input to create the same line length either side of the transect-centerline intersection. ArcPy doesnt seem to like using an Update Cursor on a Feature Layer from a GeoPackage so we need to get the full filepath to the layer in the GeoPackage.
## get half the intended full length of the transect as we will construct a
## line from the intersection of the orginal transect the same length either
## side of the centerline
half_length = line_length / 2
## arcpy doesnt like taking in a GeoPackage layer from the Content and
## using the Update Cursor so we get the full filepath.
transects = arcpy.da.Describe(transects)["catalogPath"]
Next, we get the geometry of the centerline as an ArcPy Polyline object. We need this to find where each transect intersects.
centerline = [row[0] for row in arcpy.da.SearchCursor(centerline,"SHAPE@")][0]
And now for the main event, we iterate through each transect and get the current geometry, we get the intersection point between the transect and the centerline as this should always be the midpoint of the line. Armed with the first, end, and intersect points of the transect we can construct a new line using the bearings of the original and the new line length from the user input. We use the Update Cursor to update each transects geometry as we iterate.
with arcpy.da.UpdateCursor(transects, "SHAPE@") as cursor:
## for each transect
for row in cursor:
## check if straight line, i.e. only contains two points
if row[0].pointCount == 2:
## get start, end, and cnterline intersction point
start_point = row[0].firstPoint
end_point = row[0].lastPoint
intersect_point = row[0].intersect(centerline, 1).firstPoint
####################################################################
## GET NEW LINE START POINT ########################################
intersect_start_angle = getAngle(intersect_point, start_point)
new_start = getLineEnds(intersect_point, intersect_start_angle, half_length)
####################################################################
## GET NEW LINE END POINT ##########################################
intersect_end_angle = getAngle(intersect_point, end_point)
new_end = getLineEnds(intersect_point, intersect_end_angle, half_length)
####################################################################
## CREATE NEW LINE GEOMETRY ########################################
new_line = arcpy.Polyline(arcpy.Array(arcpy.Point(*coords) for coords in [new_start, new_end]))
####################################################################
## UPDATE TRANSECT GEOMETRY ########################################
row[0] = new_line
cursor.updateRow(row)
Save your script and open up ArcGIS Pro. Right-click on your toolbox/toolset of choice and select New > Script. The New Script window will appear. In the General tab set Name to resizeTransects, Label to Resize Centerline Transects , and the Description as per below or any description you feel is apt.
Set the Input Transect Linear Features with a Data Type of Feature Layer and set the Filter to Polyline geometry only and do the same for the Input Centerline Feature parameter. Set the Line Length with a Data Type of Double.
In the Execution tab, click the folder icon in the top-right corner and add your saved Python script.
You can download the tool and other custom tools over on this page. This tool is in the Custom Tools on a Basic License with ArcPy section.
Here is a before…
…and after with shorter transects to 200m in length.
Here is the entire code with links to Esri documents for any ArcPy geoprocessing tool, function, or object used.
import arcpy
import math
################################################################################
## Esri Documentation
## https://pro.arcgis.com/en/pro-app/3.2/arcpy/functions/getparameter.htm
## https://pro.arcgis.com/en/pro-app/3.2/arcpy/functions/describe.htm
## https://pro.arcgis.com/en/pro-app/3.2/arcpy/classes/point.htm
## https://pro.arcgis.com/en/pro-app/3.2/arcpy/classes/multipoint.htm
## https://pro.arcgis.com/en/pro-app/3.2/arcpy/classes/polyline.htm
## https://pro.arcgis.com/en/pro-app/3.2/arcpy/classes/array.htm
## https://pro.arcgis.com/en/pro-app/3.2/arcpy/data-access/searchcursor-class.htm
## https://pro.arcgis.com/en/pro-app/3.2/arcpy/data-access/updatecursor-class.htm
##
## ArcGIS Pro Version 3.2.0
##
#################################################################################
#################################################################################
## USER INPUTS
## input transect linear features
transects = arcpy.GetParameter(0)
## input centerline
centerline = arcpy.GetParameter(1)
## the length to make each transect
line_length = arcpy.GetParameter(2)
################################################################################
## FUNCTIONS ###################################################################
## http://wikicode.wikidot.com/get-angle-of-line-between-two-points
def getAngle(pt1, pt2):
"""
Get the angle between two points
Args:
pt1 an arcpy Point object
pt2 an arcpy Point object
Returns:
A float representing the angle between pt1 and pt2
"""
## pt2 X value - pt1 X value
x_diff = pt2.X - pt1.X
## pt2 Y value - pt1 Y value
y_diff = pt2.Y - pt1.Y
return math.degrees(math.atan2(y_diff, x_diff))
## get a point based on bearing and distance from another point
## in our case the other point is the intersection of the transect and centerline
def getLineEnds(pt, bearing, dist):
"""
Get the start/end point of the new transect line
Args:
pt1 an arpy Point object
bearing float representing an angle between two points on the line
dist distance along the bearing where we want to locate the point
Returns:
A tuple representing x,y coords of start/end of new transect line
"""
## convert to radians
bearing = math.radians(bearing)
## some mathsy stuff that I dont have a clue about
x = pt.X + dist * math.cos(bearing)
y = pt.Y + dist * math.sin(bearing)
return (x, y)
################################################################################
## REQUIRED OBJECTS ############################################################
## get half the intended full length of the transect as we will construct a
## line from the intersection of the orginal transect the same length either
## side of the centerline
half_length = line_length / 2
## arcpy doesnt like taking in a GeoPackage layer from the Content and
## using the Update Cursor so we get the full filepath.
transects = arcpy.da.Describe(transects)["catalogPath"]
################################################################################
## GET THE CENTERLINE AS AN ARCPY POLYLINE OBJECT ##############################
centerline = [row[0] for row in arcpy.da.SearchCursor(centerline,"SHAPE@")][0]
################################################################################
## CREATE NEW GEOMETRY FOR EACH TRANSECT ######################################
with arcpy.da.UpdateCursor(transects, "SHAPE@") as cursor:
## for each transect
for row in cursor:
## check if straight line, i.e. only contains two points
if row[0].pointCount == 2:
## get start, end, and cnterline intersction point
start_point = row[0].firstPoint
end_point = row[0].lastPoint
intersect_point = row[0].intersect(centerline, 1).firstPoint
####################################################################
## GET NEW LINE START POINT ########################################
intersect_start_angle = getAngle(intersect_point, start_point)
new_start = getLineEnds(intersect_point, intersect_start_angle, half_length)
####################################################################
## GET NEW LINE END POINT ##########################################
intersect_end_angle = getAngle(intersect_point, end_point)
new_end = getLineEnds(intersect_point, intersect_end_angle, half_length)
####################################################################
## CREATE NEW LINE GEOMETRY ########################################
new_line = arcpy.Polyline(arcpy.Array(arcpy.Point(*coords) for coords in [new_start, new_end]))
####################################################################
## UPDATE TRANSECT GEOMETRY ########################################
row[0] = new_line
cursor.updateRow(row)
################################################################################