¿Encontrar el ángulo entre entidades que se cruzan en dos clases de entidad usando ArcGIS Desktop y Python? [cerrado]

Tengo dos clases de entidad de línea de intersección. Quiero encontrar el ángulo en cada punto de intersección usando ArcGIS 10 y Python.

¿Alguien puede ayudar?

Hay un flujo de trabajo relativamente simple. Supera los problemas potenciales de que dos características pueden cruzarse en más de un punto. No requiere secuencias de comandos (pero puede convertirse fácilmente en una secuencia de comandos). Se puede hacer principalmente desde el menú de ArcGIS.

La idea es explotar una capa de puntos de intersección, un punto para cada par distinto de polilíneas de intersección. Debe obtener una pequeña parte de cada polilínea de intersección en estos puntos de intersección. Usa las orientaciones de estas piezas para calcular sus ángulos de intersección.

Aquí están los pasos:

1. Asegúrese de que cada una de las características de polilínea tenga un identificador único dentro de su tabla de atributos. Esto se usará más adelante para unir algunos atributos geométricos de las polilíneas a la tabla de puntos de intersección.

2. Geoprocesamiento | Intersect obtiene los puntos (asegúrese de especificar que desea puntos para la salida).

3. Geoprocesamiento | Buffer le permite almacenar los puntos en una pequeña cantidad. Hazlo realmente pequeño para que la porción de cada línea dentro de un búfer no se doble.

4. Geoprocesamiento | Clip (aplicado dos veces) limita las capas de polilínea originales solo a los buffers. Debido a que esto produce nuevos conjuntos de datos para su salida, las operaciones posteriores no modificarán los datos originales (lo cual es bueno).

   

   Aquí hay un esquema de lo que sucede: dos capas de polilínea, que se muestran en azul claro y rojo claro, han producido puntos de intersección oscuros. Alrededor de esos puntos se muestran pequeños tampones en amarillo. Los segmentos azul y rojo más oscuros muestran los resultados de recortar las características originales en estos búferes. El resto del algoritmo funciona con los segmentos oscuros. (No puede verlo aquí, pero una pequeña polilínea roja interseca dos de las líneas azules en un punto común, produciendo lo que parece ser un amortiguador alrededor de dos polilíneas azules. En realidad, son dos amortiguadores alrededor de dos puntos superpuestos de intersección rojo-azul. Por lo tanto , este diagrama muestra cinco buffers en total).

5. Use la herramienta AddField para crear cuatro nuevos campos en cada una de estas capas recortadas: [X0], [Y0], [X1] e [Y1]. Tendrán coordenadas de puntos, así que hágales dobles y déles mucha precisión.

6. Calcular geometría (invocado al hacer clic con el botón derecho en cada nuevo encabezado de campo) le permite calcular las coordenadas x e y de los puntos inicial y final de cada polilínea recortada: colóquelas en [X0], [Y0], [X1] , y [Y1], respectivamente. Esto se hace para cada capa recortada, por lo que se necesitan 8 cálculos.

7. Use la herramienta AddField para crear un nuevo campo [Ángulo] en la capa del punto de intersección.

8. Unir las tablas recortadas a la tabla de puntos de intersección en función de identificadores de objetos comunes. (Las uniones se realizan haciendo clic con el botón derecho en el nombre de la capa y seleccionando "Uniones y relaciones").

   En este punto, la tabla de intersección de puntos tiene 9 campos nuevos: dos se llaman [X0], etc., y uno se llama [Ángulo]. Alias los campos [X0], [Y0], [X1] e [Y1] que pertenecen a una de las tablas unidas. Llamemos a estos (digamos) "X0a", "Y0a", "X1a" e "Y1a".

9. Use la Calculadora de campo para calcular el ángulo en la tabla de intersección. Aquí hay un bloque de código Python para el cálculo:

   dx = !x1!-!x0!
   dy = !y1!-!y0!
   dxa = !x1a!-!x0a!
   dya = !y1a!-!y0a!
   r = math.sqrt(math.pow(dx,2) + math.pow(dy,2))
   ra = math.sqrt(math.pow(dxa,2) + math.pow(dya,2))
   c = math.asin(abs((dx*dya - dy*dxa))/(r*ra)) / math.pi * 180

   La expresión de cálculo de campo es, por supuesto, simplemente

   c

A pesar de la longitud de este bloque de código, la matemática es simple: (dx, dy) es un vector de dirección para la primera polilínea y (dxa, dya) es un vector de dirección para la segunda. Sus longitudes, r y ra (calculadas mediante el teorema de Pitágoras), se utilizan para normalizarlas a vectores unitarios. (No debería haber ningún problema con longitudes cero, porque el recorte debería producir características de longitud positiva). El tamaño de su producto de cuña dx dya - dydxa (después de la división por r y ra) es el seno del ángulo. (El uso del producto de cuña en lugar del producto interno habitual debería proporcionar una mejor precisión numérica para ángulos cercanos a cero). Finalmente, el ángulo se convierte de radianes a grados. El resultado estará entre 0 y 90. Tenga en cuenta la evitación de la trigonometría hasta el final: este enfoque tiende a producir resultados confiables y fáciles de calcular.

Algunos puntos pueden aparecer varias veces en la capa de intersección. Si es así, obtendrán múltiples ángulos asociados con ellos.

El almacenamiento en búfer y el recorte en esta solución son relativamente caros (pasos 3 y 4): no desea hacerlo de esta manera cuando están involucrados millones de puntos de intersección. Lo recomendé porque (a) simplifica el proceso de encontrar dos puntos sucesivos a lo largo de cada polilínea dentro del vecindario de su punto de intersección y (b) el almacenamiento en búfer es tan básico que es fácil de hacer en cualquier SIG; no se necesitan licencias adicionales por encima del nivel básico de ArcMap, y generalmente produce resultados correctos. (Otras operaciones de "geoprocesamiento" podrían no ser tan confiables).

Creo que necesitas crear un script de Python.

Puede hacerlo utilizando herramientas de geoprocesamiento y arcpy.

Aquí están las principales herramientas e ideas que pueden ser útiles para usted:

1. Haga una intersección de sus dos polilíneas (llamémoslas PLINE_FC1, PLINE_FC2) entidades de clase (como resultado, necesita entidades de punto - POINT_FC) usando la herramienta Intersecar . Tendrá ID de PLINE_FC1, PLINE_FC2 en los puntos POINT_FC.
2. Dividir PLINE_FC1 por POINT_FC usando la herramienta Dividir línea en el punto. Como resultado, tendrá polilíneas divididas, la principal ventaja de esto es que puede tomar el primer / último vértice de dicha línea y compararlo con el vértice siguiente / anterior (diferencia de coordenadas) y calcular el ángulo. Entonces, tendrá un ángulo de su línea en el punto de intersección. Aquí hay un problema: debe ejecutar esta herramienta manualmente varias veces para darse cuenta de cómo se escribe la salida. Quiero decir que si toma polilínea, divídala, escriba dos polilíneas de resultado en la salida y luego proceda a la siguiente polilínea y repita. O puede ser que esta parte (resultado de la división) se escriba en diferentes clases de entidad de memoria y luego se agregue a la salida. Este es el problema principal: darse cuenta de cómo se escribe la salida para poder filtrar solo la primera parte de cada polilínea después de la división. Otra posible solución es recorrer todas las polilíneas divididas conSearchCursor y tome solo el primero encontrado (por ID de las polilíneas de origen PLINE_FC1).
3. Para obtener el ángulo, necesitará acceder a las vértices de la polilínea resultante utilizando arcpy . Escriba los ángulos resultantes en puntos (POINT_FC).
4. Repita los pasos 2-3 para PLINE_FC2.
5. Resta los atributos de ángulo (en POINT_FC) y obtén el resultado.

Puede ser que sea muy difícil codificar el paso 2 (también algunas herramientas requieren licencia ArcInfo). Luego, también puede intentar analizar los vértices de cada polilínea (agrupándolos por ID después de la intersección).

Aquí está la forma de hacerlo:

1. Tome el primer punto de intersección POINT_FC. Obtenga sus coordenadas ( point_x, point_y)
2. Por su ID, tome la polilínea fuente respectiva de PLINE_FC1.
3. Tome la primera ( vert0_x, vert0_y) y la segunda ( vert1_x, vert1_y) vértices de la misma.
4. Para el primer vértice, calcule la tangente de línea entre este vértice y el punto de intersección: tan0 = (point_y - vert0_y) / (point_x - vert0_x)
5. Calcule lo mismo para el segundo vértice: tan1 = (vert1_y - point_y) / (vert1_x - point_x)
6. Si tan1es igual a tan2, entonces ha encontrado dos vértices de su línea que tienen un punto de intersección en el medio y puede calcular el ángulo de intersección para esta línea. De lo contrario, debe proceder al siguiente par de vértices (segundo, tercero) y así sucesivamente.
7. Repita los pasos 1-6 para cada punto de intersección.
8. Repita los pasos 1-7 para la segunda clase de entidad de polilínea PLINE_FC2.
9. Resta los atributos de ángulo de PLINE_FC1 y PLINE_FC2 y obtén el resultado.

Recientemente estaba tratando de hacerlo por mi cuenta.

Mi característica de pista se basa en puntos circulares alrededor de las intersecciones de líneas, así como en puntos ubicados a una distancia de un metro de las intersicciones. La salida es una clase de entidad de polilínea que tiene atributos de número de ángulos en intersecciones y ángulos.

Tenga en cuenta que las líneas deben planificarse para encontrar intersecciones y la referencia espacial debe establecerse con la visualización correcta de la longitud de línea (la mía es WGS_1984_Web_Mercator_Auxiliary_Sphere).

Se ejecuta en la consola ArcMap pero se puede convertir fácilmente en un script en la caja de herramientas. Este script usa solo una capa de línea en TOC, nada más.

import arcpy
import time

mxd = arcpy.mapping.MapDocument("CURRENT")
df = mxd.activeDataFrame


line = ' * YOUR POLYLINE FEATURE LAYER * ' # paste the name of line layer here    

def crossing_cors(line_layer):
    mxd = arcpy.mapping.MapDocument("CURRENT")
    df = mxd.activeDataFrame
    arcpy.env.overwriteOutput = True
    sr = arcpy.Describe(line_layer).spatialReference

    dict_cors = {}
    dang_list = []

    with arcpy.da.UpdateCursor(line_layer, ['SHAPE@', 'OID@']) as uc:
        for row in uc:
            if row[0] is None:
                uc.deleteRow()

    with arcpy.da.UpdateCursor(line_layer, 'SHAPE@', spatial_reference = sr) as uc:
        for row in uc:
            line = row[0].getPart(0)
            for cor in line:
                coord = (cor.X, cor.Y)
                try:
                    dict_cors[coord] += 1
                except:
                    dict_cors[coord] = 1
    cors_only = [f for f in dict_cors if dict_cors[f]!=1]
    cors_layer = arcpy.CreateFeatureclass_management('in_memory', 'cross_pnt', "POINT", spatial_reference = sr)
    arcpy.AddField_management(cors_layer[0], 'ANGLE_NUM', 'LONG')
    with arcpy.da.InsertCursor(cors_layer[0], ['SHAPE@', 'ANGLE_NUM']) as ic:
        for x in cors_only:
            pnt_geom = arcpy.PointGeometry(arcpy.Point(x[0], x[1]), sr)
            ic.insertRow([pnt_geom, dict_cors[x]])
    return cors_layer

def one_meter_dist(line_layer):
    mxd = arcpy.mapping.MapDocument("CURRENT")
    df = mxd.activeDataFrame
    arcpy.env.overwriteOutput = True
    sr = arcpy.Describe(line_layer).spatialReference

    dict_cors = {}
    dang_list = []
    cors_list = []
    with arcpy.da.UpdateCursor(line_layer, 'SHAPE@', spatial_reference = sr) as uc:
        for row in uc:
            line = row[0]
            length_line = line.length 
            if length_line > 2.0:
                pnt1 = line.positionAlongLine(1.0)
                pnt2 = line.positionAlongLine(length_line - 1.0)
                cors_list.append(pnt1)
                cors_list.append(pnt2)
            else:
                pnt = line.positionAlongLine(0.5, True)
    cors_layer = arcpy.CreateFeatureclass_management('in_memory', 'cross_one_meter', "POINT", spatial_reference = sr)
    ic = arcpy.da.InsertCursor(cors_layer[0], 'SHAPE@')
    for x in cors_list:
        ic.insertRow([x])
    return cors_layer

def circles(pnts):

    import math
    mxd = arcpy.mapping.MapDocument("CURRENT")
    df = mxd.activeDataFrame
    arcpy.env.overwriteOutput = True
    sr = df.spatialReference

    circle_layer = arcpy.CreateFeatureclass_management('in_memory', 'circles', "POINT", spatial_reference = sr)


    ic = arcpy.da.InsertCursor(circle_layer[0], 'SHAPE@')
    with arcpy.da.SearchCursor(pnts, 'SHAPE@', spatial_reference = sr) as sc:
        for row in sc:
            fp = row[0].centroid
            list_circle =[]
            for i in xrange(0,36):
                an = math.radians(i * 10)
                np_x = fp.X + (1* math.sin(an))
                np_y = fp.Y + (1* math.cos(an))
                pnt_new = arcpy.PointGeometry(arcpy.Point(np_x,np_y), sr)

                ic.insertRow([pnt_new])
    del ic 
    return circle_layer

def angles(centers, pnts, rnd):
    mxd = arcpy.mapping.MapDocument("CURRENT")
    df = mxd.activeDataFrame
    sr = df.spatialReference

    line_lyr = arcpy.CreateFeatureclass_management('in_memory', 'line_angles', "POLYLINE", spatial_reference = sr)
    arcpy.AddField_management(line_lyr[0], 'ANGLE', "DOUBLE")
    arcpy.AddField_management(line_lyr[0], 'ANGLE_COUNT', "LONG")

    ic = arcpy.da.InsertCursor(line_lyr[0], ['SHAPE@', 'ANGLE', 'ANGLE_COUNT'])

    arcpy.AddField_management(pnts, 'ID_CENT', "LONG")
    arcpy.AddField_management(pnts, 'CENT_X', "DOUBLE")
    arcpy.AddField_management(pnts, 'CENT_Y', "DOUBLE")
    arcpy.Near_analysis(pnts, centers,'',"LOCATION") 

    with arcpy.da.UpdateCursor(line, ['SHAPE@', 'OID@']) as uc:
        for row in uc:
            if row[0] is None:
                uc.deleteRow()

    with arcpy.da.UpdateCursor(pnts, [u'ID_CENT', u'CENT_X', u'CENT_Y', u'NEAR_FID', u'NEAR_DIST', u'NEAR_X', u'NEAR_Y'], spatial_reference = sr) as uc:
        for row in uc:
            row[0] = row[3]
            row[1] = row[5]
            row[2] = row[6]
            uc.updateRow(row)
            if row[4] > 1.1:
                uc.deleteRow()


    arcpy.Near_analysis(pnts, rnd,'',"LOCATION")     

    list_id_cent = []
    with arcpy.da.UpdateCursor(pnts, [u'ID_CENT', u'CENT_X', u'CENT_Y', u'NEAR_FID', u'NEAR_DIST', u'NEAR_X', u'NEAR_Y', 'SHAPE@'], spatial_reference = sr) as uc:
        for row in uc:
            pnt_init = (row[-1].centroid.X, row[-1].centroid.Y)
            list_id_cent.append([(row[1], row[2]), row[3], pnt_init])

    list_id_cent.sort()
    values = set(map(lambda x:x[0], list_id_cent))
    newlist = [[y for y in list_id_cent if y[0]==x] for x in values]

    dict_cent_angle = {}

    for comp in newlist:
        dict_ang = {}
        for i, val in enumerate(comp):

            curr_pnt = comp[i][2]
            prev_p = comp[i-1][2]
            init_p = comp[i][0]


            angle_prev = math.degrees(math.atan2(prev_p[1]-init_p[1], prev_p[0]-init_p[0]))
            angle_next = math.degrees(math.atan2(curr_pnt[1]-init_p[1], curr_pnt[0]-init_p[0]))

            diff = abs(angle_next-angle_prev)%180


            vec1 = [(curr_pnt[0] - init_p[0]), (curr_pnt[1] - init_p[1])]
            vec2 = [(prev_p[0] - init_p[0]), (prev_p[1] - init_p[1])]

            ab = (vec1[0] * vec2[0]) + (vec1[1] * vec2[1]) 
            mod_ab = math.sqrt(math.pow(vec1[0], 2) + math.pow(vec1[1], 2)) * math.sqrt(math.pow(vec2[0], 2) + math.pow(vec2[1], 2))
            cos_a = round(ab/mod_ab, 2)

            diff = math.degrees(math.acos(cos_a))

            pnt1 = arcpy.Point(prev_p[0], prev_p[1])
            pnt2 = arcpy.Point(init_p[0], init_p[1])
            pnt3 = arcpy.Point(curr_pnt[0], curr_pnt[1])


            line_ar = arcpy.Array([pnt1, pnt2, pnt3])
            line_geom = arcpy.Polyline(line_ar, sr)

            ic.insertRow([line_geom , diff, len(comp)])
    del ic

    lyr_lst = [f.name for f in arcpy.mapping.ListLayers(mxd)]
    if 'line_angles' not in lyr_lst:
        arcpy.mapping.AddLayer(df, arcpy.mapping.Layer(line_lyr[0]))


centers = crossing_cors(line)

pnts = one_meter_dist(line)

rnd = circles(centers)

angle_dict = angles(centers, pnts, rnd)