# Interpolating points along line in Python for ArcGIS Desktop?

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
`````` 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.

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::

• 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

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.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.positionAlongLine(distance).firstPoint.X
yPoint = row.positionAlongLine(distance).firstPoint.Y
xy = (xPoint, yPoint)
insertCursor.insertRow([xy])
``````
• This was flagged as a duplicate answer (gis.stackexchange.com/a/143746/8104). Please see this post on duplicates: 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

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

# 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
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()
``````

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.

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.