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.

Advertisements

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.
  2. # 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]
    
  3. 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
    
  4. 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."