Saltar al contenido

¿Encontrar líneas cercanas y paralelas usando ArcPy?

Luego de tanto luchar hemos hallado la solución de este inconveniente que muchos usuarios de esta web han tenido. Si deseas aportar algo no dudes en dejar tu comentario.

Solución:

Aqui esta la respuesta

El rumbo no se aplica a una línea completa, sino a un segmento (línea de 2 puntos), por lo que el primer paso es dividir las polilíneas en dos líneas de puntos. Con los segmentos aplique el ángulo, pero simplificado a solo los dos primeros cuadrantes del círculo (0-180 grados) volteando las líneas si la coordenada Y es más baja al final – las líneas van ambos formas. Luego genere una tabla cercana, compare los ángulos, únase y exporte.

Comencé con algunos datos que tenía por ahí:
Solo canalización

Los resultados:
Tuberías y demasiado cerca

los líneas demasiado cercanas y paralelas son rojos. Las líneas paralelas pero sin marcar están fuera de la distancia (usé 2 metros), están a unos diez metros de distancia.

También encuentra “patas de perro”:
ingrese la descripción de la imagen aquí

No hay ninguna razón para que esos dos vértices estén allí.

¡Ahora, el código!

import sys, os, arcpy, math, time

InFeatureClass = sys.argv[1] # Input, feature class
SearchDistance = sys.argv[2] # Input, distance unit or number
OutFeatureClass = sys.argv[3] # Output, feature class

TempDir = os.environ.get("Temp")
# assumed values
DBname = "Temp_" + time.strftime("%Y-%m-%d_%H%M%S") + ".gdb"
AngleTolerance = 10
FlipAngle = 180 - AngleTolerance

TempDB = TempDir + "\" + DBname
SegmentsFC = TempDB + "\Segments"
NearTable = TempDB + "\NearTable"
NearDist = TempDB + "\NearDist"
NearAngle = TempDB + "\NearDistAngle"

def GetAngle(FromX,FromY,ToX,ToY):
    # Calculate the angle of the segment
    # simplified in the 0 to pi range
    if FromY < ToY:
        fX = FromX
        fY = FromY
        tX = ToX
        tY = ToY
    else:
        fX = ToX
        fY = ToY
        tX = FromX
        tY = FromY
    dX = tX - fX
    dY = tY - fY
    RadAngle = math.atan2(dY,dX)
    DegAngle = (RadAngle * 180) / math.pi
    return DegAngle # In degrees, would work in radians but the numbers are easier

# Get some important information about the input
desc = arcpy.Describe(InFeatureClass)
SR = desc.spatialReference
ShapeName = desc.shapeFieldName

# create temp database and feature class
arcpy.AddMessage("Creating temp database and feature class")
arcpy.CreateFileGDB_management(TempDir,DBname)
arcpy.CreateFeatureclass_management(TempDB,"Segments","POLYLINE",None,"DISABLED","DISABLED",SR)
arcpy.AddField_management(SegmentsFC,"Angle","DOUBLE")
# break the input lines into segments
sCur = arcpy.SearchCursor(InFeatureClass)
iCur = arcpy.InsertCursor(SegmentsFC,SR)
lArray = arcpy.Array()

arcpy.AddMessage("Breaking lines into segments")

for fromRow in sCur:
    geom = fromRow.getValue(ShapeName)
    partNum = 0
    for part in geom:
        thisPart = geom.getPart(partNum)
        firstPoint = True
        for pnt in thisPart:
            if pnt != None:
                if firstPoint:
                    # Skipping the first point, a segment needs two
                    # points therefore segments = points - 1
                    firstPoint = False
                    prePoint = pnt
                else:
                    # use the previous point and create a segment
                    lArray.add(prePoint)
                    lArray.add(pnt)

                    # Create a new row buffer and set shape
                    feat = iCur.newRow()
                    feat.shape = lArray
                    
                    # Set the angle of the segment
                    feat.setValue("Angle",GetAngle(prePoint.X,prePoint.Y,pnt.X,pnt.Y))

                    # insert the new feature and clear the array
                    iCur.insertRow(feat)
                    lArray.removeAll()
                    
                    prePoint = pnt # save the pnt for the next segment
        partNum += 1
del sCur
del iCur

# Generate near table of ALL features within search distance
arcpy.AddMessage("Generating near table")
arcpy.GenerateNearTable_analysis(SegmentsFC,SegmentsFC,NearTable,SearchDistance,"NO_LOCATION","NO_ANGLE","ALL")
# reduce the near table to just the non-touching features
arcpy.TableSelect_analysis(NearTable,NearDist,"NEAR_DIST > 0")
# add fields for from feature angle, to feature angle
arcpy.AddField_management(NearDist,"FAngle","DOUBLE")
arcpy.AddField_management(NearDist,"TAngle","DOUBLE")
arcpy.AddField_management(NearDist,"AngDiff","DOUBLE")

# create a join to copy the angles to the from and to angle fields
arcpy.AddMessage("Copying angles")
arcpy.MakeTableView_management(NearDist,"ND")
arcpy.AddJoin_management("ND","IN_FID",SegmentsFC,"OBJECTID")
arcpy.CalculateField_management("ND","NearDist.FAngle","!Segments.Angle!","PYTHON")
arcpy.RemoveJoin_management("ND")

arcpy.AddJoin_management("ND","NEAR_FID",SegmentsFC,"OBJECTID")
arcpy.CalculateField_management("ND","NearDist.TAngle","!Segments.Angle!","PYTHON")
arcpy.RemoveJoin_management("ND")

# calculate the difference in angle
arcpy.AddMessage("Resolving differences in angle")
arcpy.CalculateField_management(NearDist,"AngDiff","abs(!FAngle! - !TAngle!)","PYTHON")
# flip where one is near 180 and the other is near 0
arcpy.MakeTableView_management(NearDist,"NDA","AngDiff > %s" % FlipAngle) 
arcpy.CalculateField_management("NDA","AngDiff","180 - !TAngle!","PYTHON")

# Reduce the near table to similar angles
arcpy.TableSelect_analysis(NearDist,NearAngle,"AngDiff < %s" % str(AngleTolerance))

# join to the table and export to OutFeatureClass
arcpy.AddMessage("Exporting records")
arcpy.MakeFeatureLayer_management(SegmentsFC,"SegFC")
arcpy.AddJoin_management("SegFC","OBJECTID",NearAngle,"IN_FID","KEEP_COMMON")
arcpy.CopyFeatures_management("SegFC",OutFeatureClass)

arcpy.AddMessage("Cleaning up")
arcpy.Delete(TempDB)

Este código se puede copiar y poner directamente en una herramienta de Python y ejecutar desde una caja de herramientas.

Me doy cuenta de que no está pidiendo una solución SQL para su problema, pero considere por un segundo la rapidez con la que este problema podría resolverse utilizando SQL. ArcGIS tiene la capacidad de integrarse con el servidor Microsoft SQL, y no hay nada que diga que no puede mover todo el procedimiento a PostGIS. El siguiente código resolverá absolutamente su problema, dado que los enlaces tienen ID de nodo de inicio y finalización únicos, y que ha establecido índices espaciales y de identidad adecuados. Esta solución se basa en una función llamada Bearing, que se detalla aquí:

/**
* @function   : Bearing
* @precis     : Returns a bearing between two point coordinates
* @version    : 1.0
* @usage      : FUNCTION Bearing(@p_dE1  float,
*                                @p_dN1 float,
*                                @p_dE2 float,
*                                @p_dN2 float )
*                RETURNS GEOMETRY
*               eg select cogo.Bearing(0,0,45,45) * (180/PI()) as Bearing;
* @param      : p_dE1     : X Ordinate of start point of bearing
* @paramtype  : p_dE1     : FLOAT
* @param      : p_dN1     : Y Ordinate of start point of bearing
* @paramtype  : p_dN1     : FLOAT
* @param      : p_dE2     : X Ordinate of end point of bearing
* @paramtype  : p_dE2     : FLOAT
* @param      : p_dN2     : Y Ordinate of end point of bearing
* @paramtype  : p_dN2     : FLOAT
* @return     : bearing   : Bearing between point 1 and 2 from 0-360 (in radians)
* @rtnType    : bearing   : Float
* @note       : Does not throw exceptions
* @note       : Assumes planar projection eg UTM.
* @history    : Simon Greener  - Feb 2005 - Original coding.
* @history    : Simon Greener  - May 2011 - Converted to SQL Server
  * @copyright  : Licensed under a Creative Commons Attribution-Share Alike 2.5 Australia License. (http://creativecommons.org/licenses/by-sa/2.5/au/)
*/
Create Function [GIS].[Bearing](@p_dE1 Float, @p_dN1 Float,
                                 @p_dE2 Float, @p_dN2 Float)
Returns Float
AS
Begin
    Declare
        @dBearing Float,
        @dEast    Float,
        @dNorth   Float;
    BEGIN
        If (@p_dE1 IS NULL OR
            @p_dN1 IS NULL OR
            @p_dE2 IS NULL OR
            @p_dE1 IS NULL ) 
           Return NULL;

        If ( (@p_dE1 = @p_dE2) AND 
             (@p_dN1 = @p_dN2) ) 
           Return NULL;

        SET @dEast  = @p_dE2 - @p_dE1;
        SET @dNorth = @p_dN2 - @p_dN1;
        If ( @dEast = 0 ) 
        Begin
            If ( @dNorth < 0 ) 
                SET @dBearing = PI();
            Else
                SET @dBearing = 0;
        End
        Else
            SET @dBearing = -aTan(@dNorth / @dEast) + PI() / 2.0;

        If ( @dEast < 0 ) 
            SET @dBearing = @dBearing + PI();

        Return @dBearing;
    End
End;

GO

Luego recurrirá a este procedimiento de selección (suponiendo que desee objetos dentro de 2 metros / pies y considere que menos de 15 grados son casi paralelos):

      SELECT AKey, BKey
      FROM
      ( 
      SELECT  *,CASE  
              WHEN ABS(ABearing - BBearing) <= 90  THEN ABS(ABearing - BBearing)
              WHEN ABS(ABearing - BBearing) <= 180 THEN 180-ABS(ABearing - BBearing)
              WHEN ABS(ABearing - BBearing) <= 270 THEN ABS(ABearing - BBearing)-180
              WHEN ABS(ABearing - BBearing) <= 360 THEN 360-ABS(ABearing - BBearing)
          END AS AngleDiff
      FROM
      (         
      SELECT  A.ID AS AKey,
          B.ID AS BKey,
          A.Bearing AS ABearing,
          B.Bearing AS BBearing
      FROM    (
            SELECT  ID,
                us_node_id,
                ds_node_id,
                GIS.Bearing(
                       Shape.STStartPoint().STX
                       ,Shape.STStartPoint().STY
                       ,Shape.STEndPoint().STX
                       ,Shape.STEndPoint().STY
                      )* 180/PI()  as Bearing,
                Shape
            FROM    YourTableOfPipes
          ) AS A
          INNER JOIN
          (
            SELECT  ID,
                us_node_id,
                ds_node_id,
                GIS.Bearing(
                       Shape.STStartPoint().STX
                       ,Shape.STStartPoint().STY
                       ,Shape.STEndPoint().STX
                       ,Shape.STEndPoint().STY
                      )* 180/PI() as Bearing,
                Shape
            FROM    YourTableOfPipes
          ) AS B
          ON  A.us_node_id <> B.us_node_id
            AND
            A.ds_node_id <> B.ds_node_id
            AND
            A.us_node_id <> B.ds_node_id
            AND
            A.ds_node_id <> B.us_node_id
            AND
            A.Shape.STDistance(B.Shape) <=2        
      ) AS Angle
      ) AS XX
      WHERE AngleDiff <= 15

Para unos pocos miles de objetos, esto debería llevar unos segundos. Un par de millones de objetos deberían tardar unos minutos. Si no tiene un identificador us / ds, simplemente haga coincidir en A.ID <> B.ID. Necesitaría los identificadores us / ds para evitar la coincidencia de conductos que están aguas arriba / aguas abajo entre sí (conectados espacialmente), pero tengo la sensación de que este tipo de coincidencia es algo que no sería una actividad objetivo.

Su objetivo final parece ser la coincidencia de segmentos de línea. Descubrimos que el salto abierto con roadmatcher fue muy útil para este tipo de proceso, y es de código abierto.

Dentro de ArcGIS, la herramienta "contraer línea dual a línea central" también podría ayudarlo a identificar esas líneas (hay un parámetro de ancho mínimo y máximo).

Reseñas y puntuaciones

¡Haz clic para puntuar esta entrada!
(Votos: 0 Promedio: 0)



Utiliza Nuestro Buscador

Deja una respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *