Interpolating points along line in Python for ArcGIS Desktop?

Question

I am building a Python geoprocessing script for ArcGIS 9.3. In the script, I have a simple LINE feature class with one row, one line.

How can I use a linear referencing system to interpolate a point along the path?

I can use a cursor to access the feature geometry, and get the vertex geometries. But I can't find anything that could help me interpolate points along it.

Note: this is ridiculously simple with Shapely:

import shapely.wkt
line = shapely.wkt.loads('LINESTRING (380 60, 360 130, 230 350, 140 410)')
line.interpolate(0.0).xy # (array('d', [380.0]), array('d', [60.0]))
line.interpolate(100.0).xy # (array('d', [346.16312174913622]), array('d', [153.41625550146173]))
# ... etc

line with point at 100 m

Is there any equivalent in ArcGIS?

I have all the common ArcGIS extensions to my disposal.

Or should I just bring the geometry over to Shapely to do the work?

The geometry processing needs to eventually go back to ArcGIS.

score 3 · Answer 1 · answered Apr 23 '15 at 15:15

In case this is helpful to others, i was able to create the following python code using arcpy which will place points at a specified interval based on an input line feature layer.

import arcpy

line_lyr = 'my_line'
pt_lyr =  "my_point"
interval = 200

insertCursor = arcpy.da.InsertCursor(pt_lyr, ["SHAPE@XY"]) # this is the pre-existing pt feature class

with arcpy.da.SearchCursor(line_lyr, ['OID@','SHAPE@','NAME']) as searchCursor: # this is the line feature on which the points will be based
    for row in searchCursor:
        lengthLine = round(row[1].length) # grab the length of the line feature, i'm using round() here to avoid weird rounding errors that prevent the numberOfPositions from being determined
        if int(lengthLine % interval) == 0:
            numberOfPositions = int(lengthLine // interval) - 1
        else:
            numberOfPositions = int(lengthLine // interval)

        print "lengthLine", lengthLine
        print "numberOfPositions", numberOfPositions
        if numberOfPositions > 0: # > 0 b/c we don't want to add a point to a line feature that is less than our interval
            for i in range(numberOfPositions): # using range, allows us to not have to worry about
                distance = (i + 1) * interval
                xPoint = row[1].positionAlongLine(distance).firstPoint.X
                yPoint = row[1].positionAlongLine(distance).firstPoint.Y
                xy = (xPoint, yPoint)
                insertCursor.insertRow([xy])

This was flagged as a duplicate answer (http://gis.stackexchange.com/a/143746/8104). Please see this post on duplicates: http://meta.stackexchange.com/q/104227 — Aaron, Apr 23 '15 at 15:22
@Damon I would link to your identical answer on the other thread as a comment to this question here, and then delete this answer. (and thanks for the code sample too! useful for me today) — matt wilkie, May 15 '15 at 19:41
I've been trying to use this script (as is). I can see it is creating the features (numberOfPositions is 17) and I can see them in the feature class using ftrCount = int(arcpy.GetCount_management(outPath + "\\" + outName).getOutput(0)) ... but the output is empty. What might I be doing wrong? I'm actually going to post a new question about this. — jbchurchill, Jan 26 '17 at 19:56

score 2 · Answer 2 · edited Apr 13 '17 at 12:34

I am not sure what is your workflow, but for interpolation of Z's between known values (all at existing vertices) I used ArcObjects IZ.InterpolateZsBetween. I've been trying to interpolate with Calibration tool previously, however this tool have a bug. I'm not sure if it fits your purpose, but see code below for IZ.InterpolateZsBetween.

# import arcobjects liberaries
esriSystem = GetModule("C:/Program Files (x86)/ArcGIS/Desktop10.0/com/esriSystem.olb") 
esriGeometry = GetModule("C:/Program Files (x86)/ArcGIS/Desktop10.0/com/esriGeometry.olb")
esriDataSourcesGDB = GetModule("C:/Program Files (x86)/ArcGIS/Desktop10.0/com/esriDataSourcesGDB.olb")
esriGeoDatabase = GetModule("C:/Program Files (x86)/ArcGIS/Desktop10.0/com/esriGeoDatabase.olb")

# open geodatabase and featureclass
pWSF = CreateObject(esriDataSourcesGDB.FileGDBWorkspaceFactory, interface=esriGeoDatabase.IWorkspaceFactory)
pWS = pWSF.OpenFromFile(str(DbPath), 0)
pFWS = pWS.QueryInterface(esriGeoDatabase.IFeatureWorkspace)
pFClass = pFWS.OpenFeatureClass(str(fcName))

# set update cursor on the featureclass
pFCursor = pFClass.Update(None, True)
pFeat = pFCursor.NextFeature()

# loop trough features in featureclass
while pFeat:
    pShape = pFeat.ShapeCopy # clone shape of current feature
    pIZ = pShape.QueryInterface(esriGeometry.IZ2) #set IZ interface on the data - allow for interpolation of the Z value
    IPointCollection = pShape.QueryInterface(esriGeometry.IPointCollection) # set IPointCollection interface on the data - allow for points manipulation within the point collection
    IPoint = CreateObject(esriGeometry.Point, interface=esriGeometry.IPoint) # create Point object with IPoint interface
    pStart = 0 # set pStart parameter to index[0]

# loop trough IPointCollection within the polyline, find pStart and pEnd point within the polyline for IZ.InterpolateZsBetween
    for i in range(IPointCollection.PointCount):
        Point = IPointCollection.QueryPoint(i, IPoint) # query for point within the IPointCollection at index i and insert it in to IPoint

# selection of the pStart and pEnd properties based on points Z value and interpolation of the vertexes within the polyline
        if i==0: # skip value at index[0]
            pass
        elif IPoint.Z != 0: # assign pEnd and pStart if Z value of the point (vertex) is larger than 0.01 (0.01 not 0 as 0 in arcgis is returned in python as 4.54747350886e-013)
            pEnd = i
            pIZ.InterpolateZsBetween(0,pStart,0,pEnd) # program assumes that is dealing with single part polylines 
            pFeat.Shape = pIZ
            pFCursor.UpdateFeature(pFeat)
            pStart = pEnd
    pFeat = pFCursor.NextFeature()

score 1 · Accepted Answer · answered Jan 10 '12 at 05:53

1

I've done this in ArcObjects with the ICurve interface (see QueryPoint method), but this old esri thread suggests ICurve isn't exposed through Python. If that's true just keep it in Shapely. ..besides, your Shapely solution seems more fun anyway. ::grin::

answered Jan 10 '12 at 05:53

elrobis

6,456
1
31
53

QueryPoint from ICurve does exactly this, and although it is possible to use these from Python using COM, it would inflate the code and complexity too much. For my solution, I transferred the vertices Shapely to process the geometry, then return my processed geometry back to ArcGIS. – Mike T Jan 12 '12 at 19:29

score 1 · Answer 4 · edited Apr 13 '17 at 12:33

1

It looks like this has become easy using Python/ArcPy in ArcGIS for Desktop at 10.1 - see Finding mid-point of line using ArcPy? which describes a new method called "positionAlongLine" that has been added to the PolyLine class.

edited Apr 13 '17 at 12:33

Community

1

answered Aug 23 '12 at 22:43

PolyGeo

65,136
29
109
338

score 0 · Answer 5 · answered Jan 12 '12 at 07:50

If you know .NET and want to stay in ESRI land you could implement the interpolation code in .NET (using ICurve) and then glue together .NET with Python using Pythonnet (I did something similar in the following blogpost http://gissolved.blogspot.com/2009/06/python-toolbox-3-pythonnet.html). But if I where you I'd try to convert your geometry to something shapely understands.

Interpolating points along line in Python for ArcGIS Desktop?

5 Answers5

Linked