Encontrando ângulo entre os recursos que se cruzam em duas classes de recurso usando o ArcGIS Desktop e o Python? [fechadas]

Eu tenho duas classes de linhas de interseção. Quero encontrar o ângulo em cada ponto de interseção usando o ArcGIS 10 e o Python.

Alguém pode ajudar?

Existe um fluxo de trabalho relativamente simples. Supera os possíveis problemas que dois recursos podem cruzar em mais de um ponto. Não requer scripts (mas pode ser facilmente transformado em um script). Isso pode ser feito principalmente no menu ArcGIS.

A idéia é explorar uma camada de pontos de interseção, um ponto para cada par distinto de polilinhas que se cruzam. Você precisa obter um pequeno pedaço de cada polilinha que se cruza nesses pontos de interseção. Use as orientações dessas peças para calcular seus ângulos de interseção.

Aqui estão os passos:

1. Verifique se cada um dos recursos da polilinha possui um identificador exclusivo em sua tabela de atributos. Isso será usado posteriormente para unir alguns atributos geométricos das polilinhas à tabela de pontos de interseção.

2. Geoprocessamento | O Intersect obtém os pontos (certifique-se de especificar os pontos desejados para a saída).

3. Geoprocessamento | Buffer permite que você armazene os pontos em uma pequena quantidade. Torne-o realmente minúsculo para que a parte de cada linha dentro de um buffer não se dobre.

4. Geoprocessamento | Clipe (aplicado duas vezes) limita as camadas de polilinha originais apenas aos buffers. Como isso produz novos conjuntos de dados para sua saída, as operações subseqüentes não modificarão os dados originais (o que é uma coisa boa).

   

   Aqui está um esquema do que acontece: duas camadas de polilinha, mostradas em azul claro e vermelho claro, produziram pontos de interseção escuros. Em torno desses pontos, pequenos buffers são mostrados em amarelo. Os segmentos azul e vermelho mais escuros mostram os resultados do recorte dos recursos originais para esses buffers. O restante do algoritmo trabalha com os segmentos escuros. (Você não pode vê-lo aqui, mas uma minúscula polilinha vermelha cruza duas das linhas azuis em um ponto comum, produzindo o que parece ser um amortecedor em torno de duas polilinhas azuis. Na verdade, são dois amortecedores em torno de dois pontos sobrepostos da interseção vermelho-azul. , este diagrama exibe cinco buffers no total.)

5. Use a ferramenta AddField para criar quatro novos campos em cada uma dessas camadas cortadas: [X0], [Y0], [X1] e [Y1]. Eles manterão as coordenadas dos pontos; portanto, faça-os dobrar e dê muita precisão.

6. Calcular geometria (chamada clicando com o botão direito do mouse em cada novo cabeçalho de campo) permite calcular as coordenadas x e y dos pontos inicial e final de cada polilinha cortada: coloque-as em [X0], [Y0], [X1] e [Y1], respectivamente. Isso é feito para cada camada cortada; portanto, são necessários 8 cálculos.

7. Use a ferramenta AddField para criar um novo campo [Ângulo] na camada de ponto de interseção.

8. Associe as tabelas cortadas à tabela de ponto de interseção com base em identificadores de objetos comuns. (As junções são realizadas clicando com o botão direito do mouse no nome da camada e selecionando "Joins and Relates".)

   Nesse ponto, a tabela de interseção de pontos possui 9 novos campos: dois são nomeados [X0], etc., e um é chamado [Ângulo]. Alias dos campos [X0], [Y0], [X1] e [Y1] que pertencem a uma das tabelas associadas. Vamos chamá-los (digamos) de "X0a", "Y0a", "X1a" e "Y1a".

9. Use a Calculadora de campo para calcular o ângulo na tabela de interseção. Aqui está um bloco de código Python para o 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

   A expressão de cálculo de campo é, obviamente, apenas

   c

Apesar do comprimento desse bloco de código, a matemática é simples: (dx, dy) é um vetor de direção para a primeira polilinha e (dxa, dya) é um vetor de direção para o segundo. Seus comprimentos, r e ra (calculados pelo Teorema de Pitágoras), são usados ​​para normalizá-los em vetores unitários. (Não deve haver nenhum problema com comprimentos zero, porque o recorte deve produzir características com comprimento positivo.) O tamanho do produto em cunha dx dya - dydxa (depois da divisão por re ra) é o seno do ângulo. (Usar o produto da cunha em vez do produto interno usual deve fornecer uma precisão numérica melhor para ângulos próximos de zero.) Finalmente, o ângulo é convertido de radianos em graus. O resultado ficará entre 0 e 90. Observe a evitação da trigonometria até o final: essa abordagem tende a produzir resultados confiáveis ​​e facilmente calculáveis.

Alguns pontos podem aparecer várias vezes na camada de interseção. Nesse caso, eles terão vários ângulos associados a eles.

O buffer e o recorte nesta solução são relativamente caros (etapas 3 e 4): você não deseja fazê-lo dessa maneira quando milhões de pontos de interseção estão envolvidos. Eu o recomendei porque (a) simplifica o processo de encontrar dois pontos sucessivos ao longo de cada polilinha na vizinhança de seu ponto de interseção e (b) o buffer é tão básico que é fácil de fazer em qualquer GIS - nenhum licenciamento adicional é necessário acima do nível básico do ArcMap - e geralmente produz resultados corretos. (Outras operações de "geoprocessamento" podem não ser tão confiáveis.)

Eu acredito que você precisa criar um script python.

Você pode fazer isso usando ferramentas de geoprocessamento e arcpy.

Aqui estão as principais ferramentas e idéias que podem ser úteis para você:

1. Faça a interseção de suas duas polilinhas (vamos chamá-las de PLINE_FC1, PLINE_FC2) com classes (você precisa de recursos de ponto como resultado - POINT_FC) usando a ferramenta Intersect . Você terá IDs de PLINE_FC1, PLINE_FC2 nos pontos POINT_FC.
2. Divida PLINE_FC1 por POINT_FC usando a ferramenta Split Line At Point. Como resultado, você terá polilinhas divididas - a principal vantagem é que você pode pegar o primeiro / o último vértice dessa linha e compará-lo com o vértice próximo / anterior (diferença de coordenadas) e calcular o ângulo. Então, você terá o ângulo da sua linha no ponto de interseção. Há um problema aqui - você precisa executar essa ferramenta manualmente várias vezes para perceber como a saída é gravada. Quero dizer, se pegar polilinha, dividir, escrever duas polilinhas de resultado na saída e prosseguir para a próxima polilinha e repetir. Ou pode ser que essa parte (resultado da divisão) seja gravada em diferentes classes de recursos de memória e depois anexada à saída. Esse é o principal problema - perceber como a saída é gravada para poder filtrar apenas a primeira parte de cada polilinha após a divisão. Outra solução possível é fazer um loop em todas as polilinhas divididas comSearchCursor e pegue apenas o primeiro encontrado (por ID das polilinhas de origem PLINE_FC1).
3. Para obter o ângulo, você precisará acessar as vertentes da polilinha de resultado usando o arco . Escreva os ângulos resultantes em pontos (POINT_FC).
4. Repita as etapas 2-3 para PLINE_FC2.
5. Subtraia atributos de ângulo (em POINT_FC) e obtenha resultado.

Pode ser que seja muito difícil codificar a etapa 2 (também algumas ferramentas requerem licença do ArcInfo). Depois, você também pode tentar analisar vertigens de cada polilinha (agrupando-as por ID após interseção).

Aqui está o caminho para fazê-lo:

1. Pegue o primeiro ponto de interseção POINT_FC. Obtenha suas coordenadas ( point_x, point_y)
2. Por seu ID, obtenha a polilinha de origem respectiva de PLINE_FC1.
3. Pegue a primeira ( vert0_x, vert0_y) e a segunda ( vert1_x, vert1_y) vertigens.
4. Para o primeiro vértice, calcule a tangente da linha entre esse vértice e o ponto de interseção: tan0 = (point_y - vert0_y) / (point_x - vert0_x)
5. Calcule a mesma coisa para o segundo vértice: tan1 = (vert1_y - point_y) / (vert1_x - point_x)
6. Se tan1for igual a tan2, você encontrou duas vertentes da sua linha com pontos de interseção no meio e pode calcular o ângulo de interseção para esta linha. Caso contrário, você deve prosseguir para o próximo par de vertentes (segunda, terceira) e assim por diante.
7. Repita as etapas de 1 a 6 para cada ponto de interseção.
8. Repita as etapas de 1 a 7 para a segunda classe de repetição de polilinha PLINE_FC2.
9. Subtraia atributos de ângulo de PLINE_FC1 e PLINE_FC2 e obtenha resultado.

Recentemente, eu estava tentando fazer isso sozinho.

Meu recurso de pista é baseado em pontos circulares em torno das interseções de linhas, bem como em pontos localizados a um metro de distância das interseções. A saída é uma classe de recurso de polilinha que possui atributos do número de ângulos em interseções e ângulos.

Observe que as linhas devem ser planarizadas para encontrar interseções e a referência espacial deve ser definida com a exibição correta do comprimento da linha (a minha é WGS_1984_Web_Mercator_Auxiliary_Sphere).

Executando no console do ArcMap, mas facilmente pode ser transformado em um script na caixa de ferramentas. Este script usa apenas a camada de linha no sumário, nada mais.

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)