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í:
Los resultados:
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”:
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).