Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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 Esriland? 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.

share|improve this question

4 Answers 4

up vote 1 down vote accepted

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

share|improve this answer
    
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

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()
share|improve this answer

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

share|improve this answer

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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.