Generate chainage along a line feature

This is a request which comes up quite often. Normally, the request is to create points along a pipeline at every 1000m (or some other integer). Sometimes it will be for chainage along a road, but it’s always a line feature which defines a route of some kind. The first time I did this, I used ET GeoWizards, which was an excellent starting point. However, with access to an ArcInfo Advanced licence, I decided to try it using only geoprocessing tools found in ArcToolbox. This did not go well. I then turned to something which I should have turned to at the beginning – Python.

	#
	# @date 27/11/2014
	# @author Cindy Williams
	#
	# Generates chainage at a set distance along line features
	# and saves it to an existing, empty point feature class.
	# The point feature class should have a predefined schema
	# in case other attributes need to be transferred.
	#
	# For use as a standalone script.
	#
	import arcpy
	arcpy.env.workspace = r"C:\Some\Arb\Folder\work.gdb"
	pts = [] # Empty list for chainage points
	route = arcpy.management.MakeFeatureLayer("ftr_route") # Feature class containing route features
	chain = arcpy.management.MakeFeatureLayer("ftr_chain") # Empty points feature class
	i = 1000 # Set the chainage amount
	# Remove previous attempts at generating chainage
	arcpy.management.TruncateTable(chain)
	print("Chainage feature class emptied.")
	# Loop over all the features in the route feature class
	with arcpy.da.SearchCursor(route,["SHAPE@", "Route_Name"]) as rows:
	for row in rows:
	print("\tGenerating chainage along " + row[1])
	leng = row[0].length # Get the length of the current feature
	pts.append((0,(row[0].positionAlongLine(0)), row[1])) # Get the start point
	n = i # Start the count for current feature
	while i < leng:
	# Add a point at every i metres
	pts.append((i,(row[0].positionAlongLine(i)), row[1]))
	n += i
	pts.append((leng,(row[0].positionAlongLine(leng)), row[1])) # Get the end point
	# Open a cursor on the points layer
	cursor = arcpy.da.InsertCursor(chain,("LEN","SHAPE@XY", "Route_Name"))
	# Insert list of points into chainage feature class
	for row in pts:
	cursor.insertRow(row)
	del pts, cursor
	print("Script completed.")

view raw

generateChainageAlongLineFeatures.py

hosted with ❤ by GitHub

Create points at set distances and split line

This week I wrote a script to create points along an existing line based on distances found in a spreadsheet, and create a new line from the original line which is the length of the last distance recorded. That was a bit of a mouthful, so I will post it here in stages. Some helpful tips were found in the ArcGIS Help (as always), as well as the always helpful ArcPy café (here and here). In reference to that last link, creating the new line would have been much easier using 10.3’s segementAlongLine, but since that’s not in general release yet, it was the long way around for me.

'''
Created on 4 Nov 2014

@author: Cindy Williams
'''

import arcpy

# Environment settings
arcpy.env.overwriteOutput = True
arcpy.env.workspace = r"C:\Some\Arb\Folder\work.gdb"

''' Inputs
'''

ftr_baseline = "ftr_baseline" # Dissolved baseline feature class

ftr_distance = "ftr_distance" # New line feature class

ftr_points = r"ftr_points" # Created points feature class

# Table with distance data
print "Converting table..."
tbl_data = arcpy.conversion.ExcelToTable(r"C:\Some\Arb\Folder\distances.xls",
                                         r"tbl_data", 
                                         "data")

sr = arcpy.Describe(ftr_baseline).spatialReference.exporttostring() # Spatial reference

lst_points = [] # List for created points

The workflow is as follows:

1. Working with only the geometry of the line is more efficient, so retrieve the shape of the first record in the existing line feature class. As there is only one record anyway, I’ve included two methods of doing this in the code. The one that’s not commented out is the preferred one.

# Load line geometry
# geom_line = [row[0] for row in arcpy.da.SearchCursor(ftr_baseline, "SHAPE@", "#", sr)][0]
geom_line = arcpy.da.SearchCursor(ftr_baseline, "SHAPE@").next()[0]

2. Create the point geometries by reading in the km values from the table, converting to metres, accumulating the distance and adding the the new point feature into the list. This is accomplished using an update cursor. Next, use an insert cursor to write the points to disk i.e. to the blank point feature class which has a matching schema.

   acc_dist = 0 # Accumulated distance
   
   # Create points along line according to distance
   with arcpy.da.SearchCursor(tbl_data, ["Day", "KM"], "#", sr) as cursor:
       for row in cursor:
           acc_dist += row[1]*1000 # Accumulate the distance
           pnt = geom_baseline.positionAlongLine(acc_dist) # Create point geometry
           lst_points.append((row[0], row[1], pnt)) # Add value to point list
           print "Appended {}".format(row[0]) 
   
   # Write points to disk
   ins_cursor = arcpy.da.InsertCursor(ftr_points, ["Day", "Dist_km", "SHAPE@XY"])
   for pnt in lst_points:
       ins_cursor.insertRow(pnt)
   del ins_cursor
   

3. Finally, get the last point that was added, split the line at this point and write it to disk. I’m not sure if there is a better way to get the last entry – I did not want to retrieve it from the new point feature class as it already exists in the list.

   # Get point geometry from the tuple in the last list entry
   end_pnt = lst_points[len(lst_points) - 1][2] 
   first_line = arcpy.management.SplitLineAtPoint(geom_baseline, end_pnt, arcpy.Geometry())[0]
   arcpy.management.CopyFeatures(first_line, ftr_distance)
   
   print "Script Complete."
   

   Full script is here.

Optimised version of projecting coordinates on the fly

Last week I posted about creating a csv file containing coordinates projected on the fly. Seeing as I had to update the files as the projection information changed, I decided I might as well optimise my code. I managed to halve the amount of code and not write anything temporary to disk as the final output I want is non-spatial anyway.

'''
Created on 22 Sep 2014

@author: Cindy Williams

Project coordinates on the fly and write to csv
'''

import arcpy
import csv
import os

arcpy.env.overwriteOutput = True

fld = r"C:\Some\Arb\Folder"
sr_cape = arcpy.SpatialReference(r"C:\Some\Arb\Folder\cape17.prj")
sr_mines = arcpy.SpatialReference(r"C:\Some\Arb\Folder\Schwarzeck17_mines.prj")
fields = ["Field3", "SHAPE@XY"]
x_constant = 89
y_constant = 2000018

# Specify topdown=True so that it does not process the created csvs
for root, dirnames, filenames in os.walk(fld, topdown=True):
    for f in filenames:
        if f.endswith(".csv"):
            print "Processing " + f
            cur_csv = os.path.join(root, f)
            lyr = "XY_Layer"
            # Create XY Layer in the original CS
            arcpy.management.MakeXYEventLayer(cur_csv, "Field1", "Field2", lyr, sr_cape)
            print "\tCreated XY Event Layer"
            # Output csv filename created from splitting current csv filename
            final_csv = os.path.join(fld, "Schwarzeck17_" + f.rpartition(" ")[2])
            with open(final_csv, 'wb') as csvfile:
                csvwriter = csv.writer(csvfile)
                # Use a cursor to project the coordinates to the new CS
                with arcpy.da.SearchCursor(lyr, fields, "#", sr_mines) as cursor:
                    for row in cursor:
                        # Apply the constant and write to csv
                        csvwriter.writerow([row[1][0] - x_constant, row[1][1] - y_constant, row[0]])
            print "Processed " + final_csv
            # Clean up
            arcpy.management.Delete(lyr)
print "Script complete."