Znaleźć kąt między przecinającymi się obiektami w dwóch klasach obiektów za pomocą ArcGIS Desktop i Python? [Zamknięte]

Mam dwie klasy cech przecinających się linii. Chcę znaleźć kąt w każdym punkcie przecięcia za pomocą ArcGIS 10 i Pythona.

Czy ktoś może pomóc?

Istnieje stosunkowo prosty przepływ pracy. Eliminuje to potencjalne problemy, które dwie cechy mogą przecinać w więcej niż jednym punkcie. Nie wymaga skryptowania (ale można go łatwo przekształcić w skrypt). Można to zrobić przede wszystkim z menu ArcGIS.

Chodzi o to, aby wykorzystać warstwę punktów przecięcia, po jednym punkcie dla każdej odrębnej pary przecinających się polilinii. Musisz uzyskać mały kawałek każdej przecinającej się polilinii w tych punktach przecięcia. Użyj orientacji tych elementów, aby obliczyć ich kąty przecięcia.

Oto kroki:

1. Upewnij się, że każda funkcja polilinii ma unikalny identyfikator w tabeli atrybutów. Zostanie to wykorzystane później do połączenia niektórych atrybutów geometrycznych polilinii z tabelą punktów przecięcia.

2. Geoprocessing | Intersect uzyskuje punkty (upewnij się, że chcesz określić punkty dla wyniku).

3. Geoprocessing | Buffer pozwala buforować punkty o niewielką ilość. Spraw, aby był naprawdę mały, aby część każdej linii w buforze się nie wyginała.

4. Geoprocessing | Clip (zastosowany dwukrotnie) ogranicza oryginalne warstwy polilinii do samych buforów. Ponieważ tworzy to nowe zestawy danych dla jego wyników, kolejne operacje nie modyfikują oryginalnych danych (co jest dobrą rzeczą).

   

   Oto schemat tego, co się dzieje: dwie warstwy polilinii, pokazane w kolorze jasnoniebieskim i jasnoczerwonym, wytworzyły ciemne punkty przecięcia. Wokół tych punktów małe bufory są zaznaczone na żółto. Ciemniejsze niebieskie i czerwone segmenty pokazują wyniki obcinania oryginalnych elementów do tych buforów. Reszta algorytmu działa z ciemnymi segmentami. (Nie widać tego tutaj, ale niewielka czerwona polilinia przecina dwie niebieskie linie we wspólnym punkcie, tworząc coś, co wydaje się być buforem wokół dwóch niebieskich polilinii. To tak naprawdę dwa bufory wokół dwóch nakładających się punktów przecięcia czerwono-niebieskiego. , ten diagram pokazuje w sumie pięć buforów.)

5. Użyj narzędzia AddField , aby utworzyć cztery nowe pola w każdej z tych przyciętych warstw: [X0], [Y0], [X1] i [Y1]. Będą trzymać współrzędne punktowe, dzięki czemu będą się podwajać i dadzą im dużą precyzję.

6. Oblicz geometrię (wywoływaną przez kliknięcie prawym przyciskiem myszy każdego nowego nagłówka pola) umożliwia obliczenie współrzędnych xiy punktów początkowych i końcowych każdej przyciętej polilinii: umieść je w [X0], [Y0], [X1] i [Y1], odpowiednio. Odbywa się to dla każdej przyciętej warstwy, dlatego potrzeba 8 obliczeń.

7. Użyj narzędzia AddField , aby utworzyć nowe pole [Kąt] w warstwie punktu przecięcia.

8. Połącz przycięte tabele z tabelą punktów przecięcia na podstawie wspólnych identyfikatorów obiektów. (Połączenia są wykonywane przez kliknięcie prawym przyciskiem myszy nazwy warstwy i wybranie „Złączenia i relacje”).

   W tym momencie tabela przecięć punktów ma 9 nowych pól: dwa mają nazwy [X0] itd., A jedno ma nazwę [Kąt]. Alias pola [X0], [Y0], [X1] i [Y1], które należą do jednej z połączonych tabel. Nazwijmy je (powiedzmy) „X0a”, „Y0a”, „X1a” i „Y1a”.

9. Użyj kalkulatora pola, aby obliczyć kąt w tabeli skrzyżowań. Oto blok kodu Pythona do obliczeń:

   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

   Wyrażenie pola obliczeniowego jest oczywiście jedynie

   c

Pomimo długości tego bloku kodu matematyka jest prosta: (dx, dy) jest wektorem kierunku pierwszej polilinii, a (dxa, dya) jest wektorem kierunku drugiej. Ich długości r i ra (obliczone za pomocą twierdzenia Pitagorasa) służą do normalizacji ich do wektorów jednostkowych. (Nie powinno być problemu z zerowymi długościami, ponieważ obcinanie powinno dawać cechy długości dodatniej.) Rozmiar ich produktu klinowego dx dya - dydxa (po podzieleniu przez r i ra) jest sinus kąta. (Użycie iloczynu klinowego zamiast zwykłego iloczynu wewnętrznego powinno zapewnić lepszą dokładność numeryczną dla kątów bliskich zeru.) Wreszcie kąt jest przekształcany z radianów na stopnie. Wynik będzie wynosił od 0 do 90. Zwróć uwagę na unikanie trygonometrii aż do samego końca: takie podejście prowadzi do uzyskania wiarygodnych i łatwych do obliczenia wyników.

Niektóre punkty mogą pojawiać się wiele razy w warstwie przecięcia. Jeśli tak, otrzymają wiele kątów z nimi związanych.

Buforowanie i obcinanie w tym rozwiązaniu są stosunkowo drogie (kroki 3 i 4): nie chcesz tego robić w ten sposób, gdy zaangażowane są miliony punktów przecięcia. Poleciłem go, ponieważ (a) upraszcza proces znajdowania dwóch kolejnych punktów wzdłuż każdej polilinii w sąsiedztwie jego punktu przecięcia i (b) buforowanie jest tak podstawowe, że łatwo jest to zrobić w dowolnym systemie informacji geograficznej - nie jest wymagane dodatkowe licencjonowanie powyżej podstawowego poziomu ArcMap - i zwykle daje prawidłowe wyniki. (Inne operacje „geoprzetwarzania” mogą nie być tak niezawodne).

Uważam, że musisz utworzyć skrypt Pythona.

Możesz to zrobić za pomocą narzędzi geoprzetwarzania i arcpy.

Oto główne narzędzia i pomysły, które mogą Ci się przydać:

1. Zrób przecięcie dwóch polilinii (nazwijmy je PLINE_FC1, PLINE_FC2) klasami funkcji (potrzebujesz w rezultacie operacji punktowych - POINT_FC) za pomocą narzędzia Przecięcie . Będziesz mieć identyfikatory z PLINE_FC1, PLINE_FC2 w punktach POINT_FC.
2. Podziel PLINE_FC1 przez POINT_FC za pomocą narzędzia Podziel linię w punkcie. W rezultacie będziesz miał podzielone polilinie - główną zaletą jest to, że możesz po prostu wziąć pierwszy / ostatni wierzchołek takiej linii, porównać ją z następnym / poprzednim wierzchołkiem (różnica współrzędnych) i obliczyć kąt. Będziesz miał kąt linii w punkcie przecięcia. Jest tutaj jeden problem - musisz uruchomić to narzędzie kilka razy ręcznie, aby zrozumieć, jak zapisywane są dane wyjściowe. Mam na myśli, że jeśli pobierze polilinię, podziel ją, zapisz dwie wynikowe polilinie na wyjściu, a następnie przejdź do następnej polilinii i powtórz. Być może ta część (wynik podziału) jest zapisywana w różnych klasach funkcji pamięci, a następnie dołączana do wyniku. Jest to główny problem - uświadomienie sobie, w jaki sposób zapisywane są dane wyjściowe, aby móc filtrować tylko pierwszą część każdej polilinii po podziale. Innym możliwym rozwiązaniem jest przejście przez wszystkie polilinie podzielone na wynikiWyszukaj kursor i weź tylko pierwszy napotkany (według identyfikatora źródłowych polilinii PLINE_FC1).
3. Aby uzyskać kąt, musisz uzyskać dostęp do wierzchołków polilinii wynikowej za pomocą arcpy . Napisz wynikowe kąty do punktów (POINT_FC).
4. Powtórz kroki 2-3 dla PLINE_FC2.
5. Odejmij atrybuty kąta (w POINT_FC) i uzyskaj wynik.

Może być bardzo trudno zakodować krok 2 (również niektóre narzędzia wymagają licencji ArcInfo). Następnie możesz także spróbować przeanalizować wierzchołki każdej polilinii (grupując je według identyfikatora po przecięciu).

Oto sposób na zrobienie tego:

1. Weź pierwszy punkt przecięcia POINT_FC. Uzyskaj jego współrzędne ( point_x, point_y)
2. Za pomocą jego identyfikatora pobierz odpowiednią źródłową polilinię z PLINE_FC1.
3. Weź pierwsze ( vert0_x, vert0_y) i drugie ( vert1_x, vert1_y) jego wierzchołki.
4. Dla pierwszego wierzchołka obliczyć styczną linii między tym wierzchołkiem a punktem przecięcia: tan0 = (point_y - vert0_y) / (point_x - vert0_x)
5. Oblicz to samo dla drugiego wierzchołka: tan1 = (vert1_y - point_y) / (vert1_x - point_x)
6. Jeśli tan1jest równy tan2, to znalazłeś dwa wierzchołki linii, które mają punkt przecięcia pomiędzy nimi i możesz obliczyć kąt przecięcia dla tej linii. W przeciwnym razie musisz przejść do następnej pary wierzchołków (druga, trzecia) i tak dalej.
7. Powtórz kroki 1-6 dla każdego punktu przecięcia.
8. Powtórz kroki 1-7 dla drugiej klasy cech polilinii PLINE_FC2.
9. Odejmij atrybuty kąta od PLINE_FC1 i PLINE_FC2 i uzyskaj wynik.

Ostatnio próbowałem to zrobić sam.

Moja wskazówka oparta jest na okrągłych punktach wokół przecięć linii, a także punktach znajdujących się w odległości jednego metra od wtargnięć. Dane wyjściowe to klasa elementów polilinii, która ma atrybuty liczby kątów na przecięciach i kącie.

Zauważ, że linie powinny być planowane w celu znalezienia skrzyżowań, a odniesienie przestrzenne musi być ustawione z poprawnym wyświetlaniem długości linii (moja to WGS_1984_Web_Mercator_Auxiliary_Sphere).

Działa w konsoli ArcMap, ale łatwo można go zmienić na skrypt w przyborniku. Ten skrypt używa tylko warstwy linii w spisie treści, nic więcej.

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)