Descubra el punto que se encuentra entre dos líneas paralelas

8

Estoy enfrentando un problema en ArcGIS. Yo trabajo en una base de datos de navegación. En nuestra base de datos, las calles de un solo carril están representadas por una sola Línea, mientras que una calle de varios carriles (calle con divisor en el centro) está representada por dos líneas paralelas (líneas de color rojo en la imagen).

Tengo un archivo de forma de puntos con algunos puntos que caen dentro de una calle de varios carriles y otros afuera.

Quiero crear un script de ArcPy que encuentre los puntos que se encuentran dentro de las calles de varios carriles. es decir, entre estas líneas paralelas (marcadas en la imagen).

No sé cómo lograr esto, ¿alguien puede ayudarme?

Un problema callejero de varios carriles

Hice algo de ejercicio y descubrí que crear un búfer en un lado de la línea puede crear dentro del polígono Multi-Lane (se muestra en la imagen).

ingrese la descripción de la imagen aquí

pero ahora el problema es que el polígono realmente cruza la línea (es decir, se superpone al límite de varios carriles). entonces atrapará puntos innecesarios. ¿Hay alguna forma de alinear este polígono con la línea de la calle?

Nota: integrar no funcionará aquí, porque también mueve la línea de la calle. Solo necesito alinear el polígono a lo largo de la línea de la calle.

Akhil Kumar
fuente
Algo así como Medir el acimut de la calle - Crear cadenas de líneas desde cada punto hacia el ángulo Azimut + 90 grados - Contar cuántas de sus líneas paralelas se cruzan con esta línea. Si cero o dos -> afuera, si uno -> Lo encontraste. Solo pensando, puede funcionar o no. Otra idea es convertir la calle de doble vía en polígono y seleccionar los puntos que se cruzan. Esto último puede ser complicado de hacer con Python. Bueno, el primero también si las calles son curvas. Pero con el búfer de un solo lado, es posible que pueda construir polígonos callejeros bastante agradables.
user30184
1
¿tienes una licencia avanzada? Sería bastante sencillo con la herramienta cercana.
radouxju
Sí, tengo licencia avanzada.
Akhil Kumar
Al principio pensé en tomar el polígono del búfer y luego en cruzar ese polígono. y descubre qué puntos caen en ese polígono intersectado. pero el mayor problema es que la distancia intermedia no es consistente en todas partes en la calle. en algún lugar solo tiene 10 metros en algún lugar alrededor de 20 metros, en ese caso la lógica de intersección del polígono
fallará
1
Haga el búfer del lado derecho de 10 m del lado izquierdo y el búfer del lado izquierdo del lado derecho. De esa manera cubres el rango de 10-20 m. Las superposiciones no causan ningún daño y también puede fusionar los polígonos primero. O haga un polígono de amortiguación de un lado aún más ancho y recórtelo intersecando con el otro lado. Usa la imaginación y juega.
user30184

Respuestas:

4

Intentaría debajo del algoritmo arcpy (¡incluso manual!)

  1. Encuentre el ancho adecuado de las calles de dos carriles: aquí es posible que necesite agrupar calles con el mismo ancho y siga el siguiente procedimiento para cada grupo.
  2. Cree el búfer en ambas líneas hacia ambas direcciones (derecha e izquierda) con ese ancho (o un poco menos que eso, para garantizar el área de la carretera).
  3. Ejecute la herramienta de intersección para obtener la región superpuesta.
  4. Ejecute Seleccionar por ubicación para seleccionar los puntos que se encuentran dentro de este polígono.
SIslam
fuente
Creo que este es el camino a seguir. Encuentre una manera fácil de unir las líneas, ya sea mediante búfer o de alguna manera cierre las líneas para hacer un solo polígono y luego seleccione dentro.
Barrett
2

Yo diría que este es un ejercicio geométrico.

CÓDIGO PSEUDO:

  • Para cada punto (punto negro) encuentre el camino más cercano y encuentre la proyección del punto en este camino (punto rojo).
  • Dibuje una línea corta (discontinua) en dirección opuesta comenzando en el punto negro
  • Encuentre si hay una intersección entre la línea corta y la carretera del mismo nombre, estrella azul. Si hay uno, el punto negro es el que buscamos.

ingrese la descripción de la imagen aquí

Como se puede ver, hay casos especiales: puntos negros dentro de un círculo:

  1. Muy sinuosa carretera de 1 línea. Esto puede eliminarse a) trabajando solo con carreteras de 2 líneas ob) asegurándose de que los FID de las carreteras que cruzan el punto rojo y la estrella son diferentes. Sin embargo, si la carretera curva tiene un cruce con otra carretera de 1 línea, esto podría no funcionar.
  2. El punto negro se encuentra en la extensión de una carretera exactamente perpendicular de 1 línea. En este caso, existe la posibilidad de que se pueda elegir una carretera de 1 carril como vecino más cercano.
  3. Punto negro se sienta en la línea.

Todos los casos anteriores son muy poco probables, sin embargo, parece que la opción más segura es trabajar solo con carreteras de 2 líneas, es decir, exportarlas a una clase de entidad separada. El caso 3 es divertido, lo dejaremos al azar, porque la distancia más corta a la línea nunca es cero verdadero, por lo tanto, se puede encontrar la dirección 'opuesta' del rayo que conecta 2 puntos.

Implementación de Python:

import arcpy, traceback, os, sys
from arcpy import env
env.overwriteoutput=True

# things to change ---------
maxD=30
mxd = arcpy.mapping.MapDocument("CURRENT")
pointLR = arcpy.mapping.ListLayers(mxd,"NODES")[0]
lineLR = arcpy.mapping.ListLayers(mxd,"LINKS")[0]
sjOneToMany=r'D:\scratch\sj2.shp'
RDNAME='street'
# -------------------------
dDest=arcpy.Describe(lineLR)
SR=dDest.spatialReference

try:
    def showPyMessage():
        arcpy.AddMessage(str(time.ctime()) + " - " + message)
    g = arcpy.Geometry()
    geometryList=arcpy.CopyFeatures_management(pointLR,g)
    n=len(geometryList)
    endPoint=arcpy.Point()

    arcpy.SpatialJoin_analysis(pointLR, lineLR,sjOneToMany,"JOIN_ONE_TO_MANY","KEEP_COMMON","","WITHIN_A_DISTANCE",maxD)
    initFidList=(-1,)
    for fid in range(n):
        query='"TARGET_FID" = %s' %str(fid)
        nearTable=arcpy.da.TableToNumPyArray(sjOneToMany,("TARGET_FID","JOIN_FID"),query)
        if len(nearTable)<2:continue
        fidLines=[int(row[1]) for row in nearTable]
        query='"FID" in %s' %str(tuple(fidLines))
        listOfLines={}
        blackPoint=geometryList[fid]
        with arcpy.da.SearchCursor(lineLR,("FID", "Shape@","STREET"),query) as rows:
            dMin=100000
            for row in rows:
                shp=row[1];dCur=blackPoint.distanceTo(shp)
                listOfLines[row[0]]=row[-2:]
                if dCur<dMin:
                    fidNear,lineNear, roadNear=row
                    dMin=dCur
            chainage=lineNear.measureOnLine(blackPoint)
            redPoint=lineNear.positionAlongLine (chainage).firstPoint
            smallD=blackPoint.distanceTo(redPoint)
            fp=blackPoint.firstPoint
            dX=(redPoint.X-fp.X)*(maxD-smallD)/smallD
            dY=(redPoint.Y-fp.Y)*(maxD-smallD)/smallD
            endPoint.X=fp.X-dX;endPoint.Y=fp.Y-dY
            dashLine=arcpy.Polyline(arcpy.Array([fp,endPoint]),SR)

            for n in listOfLines:
                if n==fidNear:continue
                line, road=listOfLines[n]
                if road!=roadNear:continue
                blueStars=dashLine.intersect(line,1)
                if blueStars.partCount==0:continue
                initFidList+=(fid,); break
    query='"FID" in %s' %str(initFidList)
    arcpy.SelectLayerByAttribute_management(pointLR, "NEW_SELECTION", query)
    arcpy.AddMessage ('\n %i point(s) found' %(len(initFidList)-1))
except:
    message = "\n*** PYTHON ERRORS *** "; showPyMessage()
    message = "Python Traceback Info: " + traceback.format_tb(sys.exc_info()[2])[0]; showPyMessage()
    message = "Python Error Info: " +  str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"; showPyMessage()            

Hay otra posible solución quizás más elegante. Implica triangulación. Avíseme si le interesa y actualizaré mi respuesta

FelixIP
fuente
Esto es bastante complejo, wow. Parece que sería mucho más simple crear un polígono a partir de las líneas y luego usar la proyección de rayos . Determinar si un punto está en una línea también debería ser sencillo.
Paul
1
Si puede crear polígonos a partir de las líneas correctas, no es necesario realizar la conversión. Seleccionar por ubicación lo hará. Sin embargo
FelixIP
¿Funcionará bien en las curvas? Solo para aclarar :)
SIslam
1
@SIslam debería funcionar incluso con curvas grandes similares al caso 1 (ver si n == fidNear: continuar) línea. Bueno, si no hay una carretera de 1 carril entrando. Sigo pensando que disolver puede ayudar, pero no siempre
FelixIP
@Islam ¡Uy! No lo hará, porque la condición (si n == fidNear: continuar) elimina los puntos que se sientan fuera de la curva, pero marca el punto dentro como uno sentado afuera. Sin embargo, se necesita un giro brusco, ¿radio menor que el ancho?
FelixIP
0

Como las calles son paralelas, supuse que se crearon con la Copy Parallelherramienta en la barra de herramientas Editar, haciendo que el par de líneas tenga la misma dirección. Luego podemos iterar sobre las coordenadas de la primera línea y agregarlas a un polígono y luego iterar sobre el reverso de la segunda línea. Definitivamente hay una mejor manera de acercarse a los pares de líneas de agarre; El enfoque OID funciona, pero no es muy bonito.

import collections
import arcpy

FC = "fc"
points = "points"
pgons = "pgons"
arcpy.env.overwriteOutput = True

def buildpoly(oid_coords):
    #create ddict of the form OID:<x1y1, x2y2, ..., xn-1yn-1, xnyn>
    ddict = collections.defaultdict(list)    
    for k,v in oid_coords:
        ddict[k].append(v)

    line1,line2 = ddict.keys()    

    #Assume that the parallel lines have same direction, so reverse the second
    arr = arcpy.Array()
    arr.extend(arcpy.Point(*pt) for pt in ddict[line1])    
    arr.extend(arcpy.Point(*pt) for pt in ddict[line2][::-1])

    return arcpy.Polygon(arr)

#id is an integer field that pairs parallel lines together
unique = list(set(t[0] for t in arcpy.da.SearchCursor(FC, "id")))
polygons = []
for uni in unique:
    polygons.append(buildpoly([r for r in row] for row in arcpy.da.SearchCursor(FC,
                                                                                ["OID@", "SHAPE@XY"],
                                                                                "id={}".format(uni),
                                                                                explode_to_points=True)))


arcpy.CopyFeatures_management(polygons, pgons)

A partir de ahí, es una llamada a Intersecar / Seleccionar capa por ubicación / lo que tienes. Tenga en cuenta que el Spolígono con forma no es perfecto ya que lo dibujé a mano alzada y hay algunos arcos que explode_to_pointsno se manejan correctamente. Solo corre Densifyo equivalente.

ingrese la descripción de la imagen aquí

Pablo
fuente
Este es el conjunto de datos de la red de carreteras, por lo tanto, las carreteras de 1 carril están conectadas a 2 carriles a través del nodo, es decir, no existen pares de características paralelas
FelixIP
Es posible que desee ampliar su solución agregando disolver por nombres de carreteras individuales primero (sin partes m) y considerar casos de 1 o 2 líneas como resultado
FelixIP
@FelixIP, no estoy muy familiarizado con los conjuntos de datos de red. Mi solución fue principalmente una prueba de concepto de cómo se puede hacer con líneas simples (OP puede extenderlo para cubrir mresolución, multiparte, etc.). No sé cómo características como esta están realmente representadas en una red.
Paul
@ Paul La carretera del mismo nombre se puede representar con segmentos de 100 en diferentes filas de la tabla. Además, el camino de doble carril podría convertirse en un solo carril en alguna parte. Disolver fallará gravemente si no hay partes que no
estén
1
@AkhilKumar, no importa si son más o menos paralelas. Esto está trazando las líneas existentes.
Paul