¿Cómo interpolar puntos a lo largo de una línea en Python para ArcGIS 9.3?

Question

Estoy construyendo un Python secuencia de comandos de geoprocesamiento de ArcGIS 9.3. En el script, tengo una simple clase de entidad de LÍNEA con una línea, una línea. ¿Cómo puedo utilizar un sistema de referencia lineal para interpolar un punto a lo largo de la ruta? Puedo utilizar un cursor para acceder a la función de la geometría, y obtener el vértice de geometrías. Pero no puedo encontrar nada que me pudiera ayudar interpolar puntos a lo largo de ella.

Nota: este es ridículamente simple con bien formada:

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

¿Hay algún equivalente en Esriland? Tengo todas las extensiones de ArcGIS a mi disposición. O debo llevar la geometría a bien formada para hacer el trabajo? La geometría de las necesidades de procesamiento para, finalmente, volver a ArcGIS.

Answer 1

En caso de que esto sea útil para otros, pude crear el siguiente código de python usando arcpy que colocará puntos en un intervalo especificado basado en una capa de característica de línea de entrada.

 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])
 

Answer 2

No estoy seguro de cuál es su flujo de trabajo, pero para la interpolación de Z entre los valores conocidos (en todos los vértices existentes) he utilizado ArcObjects IZ.InterpolateZsBetween. He estado tratando de interpolar con la herramienta de Calibración previamente, sin embargo, esta herramienta tiene un error. No estoy seguro de si se ajusta a su propósito, pero véase el código de abajo para 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()

Answer 3

Parece que esto se ha convertido en fácil usando Python / ArcPy en ArcGIS para Desktop en 10.1 - vea Cómo encontrar el punto medio de una línea en arcpy que describe un nuevo método llamado "positionAlongLine" que se ha agregado a la clase PolyLine.

Answer 4

He hecho esto en ArcObjects con la interfaz ICurve (ver el método QueryPoint), pero este viejo hilo esri sugiere que ICurve no está expuesto a través de Python. Si eso es cierto, sólo manténgalo en Shapely. Además, su solución Shapely parece más divertida de todos modos. ::mueca::

Answer 5

Si conoces .NET y quieres quedarte en ESRI, puedes implementar el código de interpolación en .NET (usando ICurve) y luego pegar .NET con Python usando Pythonnet (hice algo similar en el siguiente blogpost http: // gissolved. Blogspot.com/2009/06/python-toolbox-3-pythonnet.html ). Pero si yo tuviera que intentar convertir tu geometría en algo bien comprendido.