Ermitteln des Winkels zwischen sich überschneidenden Features in zwei Feature-Classes mithilfe von ArcGIS Desktop und Python? [geschlossen]

Ich habe zwei sich überschneidende Linien-Feature-Classes. Ich möchte den Winkel an jedem Schnittpunkt mit ArcGIS 10 und Python ermitteln.

Kann jemand helfen?

Es gibt einen relativ einfachen Workflow. Es überwindet die potenziellen Probleme, die auftreten können, wenn sich zwei Features in mehr als einem Punkt überschneiden. Es erfordert kein Scripting (kann aber leicht in ein Script umgewandelt werden). Dies kann hauptsächlich über das ArcGIS-Menü erfolgen.

Die Idee ist, eine Schicht von Schnittpunkten auszunutzen, einen Punkt für jedes einzelne Paar sich schneidender Polylinien. Sie müssen ein kleines Stück jeder sich schneidenden Polylinie an diesen Schnittpunkten erhalten. Verwenden Sie die Ausrichtungen dieser Teile, um ihre Schnittwinkel zu berechnen.

Hier sind die Schritte:

1. Stellen Sie sicher, dass jedes der Polylinien-Features eine eindeutige Kennung in seiner Attributtabelle hat. Dies wird später verwendet, um einige geometrische Attribute der Polylinien mit der Schnittpunkttabelle zu verbinden.

2. Geoverarbeitung | Schnitt erhält die Punkte (stellen Sie sicher, dass Sie Punkte für die Ausgabe angeben ).

3. Mit Geoverarbeitung | Puffer können Sie die Punkte um einen winzigen Betrag puffern. Machen Sie es wirklich winzig, damit sich der Teil jeder Zeile innerhalb eines Puffers nicht verbiegt.

4. Geoverarbeitung | Clip (zweimal angewendet) beschränkt die ursprünglichen Polylinienebenen nur auf die Puffer. Da hierdurch neue Datensätze für die Ausgabe erstellt werden, werden durch nachfolgende Vorgänge die ursprünglichen Daten nicht geändert (was eine gute Sache ist).

   

   Hier ist eine schematische Darstellung dessen, was passiert: Zwei hellblaue und hellrote Polylinienebenen haben dunkle Schnittpunkte erzeugt. Um diese Punkte sind winzige Puffer in gelb dargestellt. Die dunkleren blauen und roten Segmente zeigen die Ergebnisse des Zuschneidens der ursprünglichen Features zu diesen Puffern. Der Rest des Algorithmus arbeitet mit den dunklen Segmenten. (Sie können es hier nicht sehen, aber eine winzige rote Polylinie schneidet zwei der blauen Linien an einem gemeinsamen Punkt und erzeugt so etwas wie einen Puffer um zwei blaue Polylinien. Es sind wirklich zwei Puffer um zwei überlappende Punkte des Rot-Blau-Schnittpunkts. Also In diesem Diagramm werden insgesamt fünf Puffer angezeigt.)

5. Verwenden Sie das AddField- Werkzeug, um vier neue Felder in jeder dieser abgeschnittenen Ebenen zu erstellen: [X0], [Y0], [X1] und [Y1]. Sie halten Punktkoordinaten, machen sie also doppelt und geben ihnen viel Präzision.

6. Mit Geometrie berechnen (durch Klicken mit der rechten Maustaste auf jeden neuen Feldkopf) können Sie die x- und y-Koordinaten der Start- und Endpunkte jeder abgeschnittenen Polylinie berechnen: Fügen Sie diese in [X0], [Y0], [X1] ein. bzw. [Y1]. Dies erfolgt für jede abgeschnittene Ebene, sodass 8 Berechnungen erforderlich sind.

7. Verwenden Sie das AddField- Werkzeug, um ein neues [Winkel] -Feld in der Schnittpunktebene zu erstellen.

8. Verknüpfen Sie die beschnittenen Tabellen mit der Schnittpunkttabelle auf der Grundlage gemeinsamer Objektkennungen. (Verknüpfungen werden ausgeführt, indem Sie mit der rechten Maustaste auf den Ebenennamen klicken und "Verknüpfungen und Verknüpfungen" auswählen.)

   Zu diesem Zeitpunkt enthält die Punktkreuzungstabelle 9 neue Felder: Zwei werden mit [X0] usw. bezeichnet, und eines wird mit [Winkel] bezeichnet. Alias der Felder [X0], [Y0], [X1] und [Y1], die zu einer der verknüpften Tabellen gehören. Nennen wir diese (sagen wir) "X0a", "Y0a", "X1a" und "Y1a".

9. Verwenden Sie den Feldrechner , um den Winkel in der Schnittpunkttabelle zu berechnen. Hier ist ein Python-Codeblock für die Berechnung:

   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

   Der Feldberechnungsausdruck ist natürlich nur

   c

Trotz der Länge dieses Codeblocks ist die Mathematik einfach: (dx, dy) ist ein Richtungsvektor für die erste Polylinie und (dxa, dya) ist ein Richtungsvektor für die zweite. Ihre Längen r und ra (berechnet über den Satz von Pythagoras) werden verwendet, um sie zu Einheitsvektoren zu normalisieren. (Bei Nulllängen sollte es kein Problem geben, da das Abschneiden Merkmale positiver Länge erzeugen sollte.) Die Größe ihres Keilprodukts dx dyadydxa (nach Division durch r und ra) ist der Sinus des Winkels. (Die Verwendung des Keilprodukts anstelle des üblichen Innenprodukts sollte eine bessere numerische Genauigkeit für Winkel nahe Null bieten.) Schließlich wird der Winkel vom Bogenmaß in Grad umgewandelt. Das Ergebnis liegt zwischen 0 und 90. Beachten Sie, dass die Trigonometrie bis zum Ende vermieden wird: Dieser Ansatz führt tendenziell zu zuverlässigen und einfach zu berechnenden Ergebnissen.

Einige Punkte werden möglicherweise mehrmals in der Schnittmenge angezeigt. In diesem Fall werden ihnen mehrere Winkel zugeordnet.

Das Puffern und Zuschneiden in dieser Lösung ist relativ teuer (Schritte 3 und 4): Sie möchten dies nicht auf diese Weise tun, wenn Millionen von Schnittpunkten betroffen sind. Ich habe es empfohlen, weil (a) es das Finden von zwei aufeinanderfolgenden Punkten entlang jeder Polylinie in der Nähe ihres Schnittpunkts vereinfacht und (b) das Puffern so einfach ist, dass es in jedem GIS einfach ist - es ist keine zusätzliche Lizenzierung erforderlich über dem ArcMap-Basisniveau - und führt normalerweise zu korrekten Ergebnissen. (Andere "Geoverarbeitungs" -Operationen sind möglicherweise nicht so zuverlässig.)

Ich glaube, Sie müssen Python-Skript erstellen.

Sie können dies mit Geoverarbeitungswerkzeugen und Arcpy tun.

Hier sind die wichtigsten Tools und Ideen, die für Sie nützlich sein können:

1. Erstellen Sie eine Schnittmenge Ihrer beiden Polylinien (nennen wir sie PLINE_FC1, PLINE_FC2) -Featureklassen (als Ergebnis benötigen Sie Punkt-Features - POINT_FC) mit dem Werkzeug " Schnittmenge" . Sie haben IDs von PLINE_FC1, PLINE_FC2 in Punkten POINT_FC.
2. Teilen Sie PLINE_FC1 durch POINT_FC mit dem Werkzeug Linie am Punkt teilen. Als Ergebnis erhalten Sie geteilte Polylinien - der Hauptvorteil besteht darin, dass Sie nur den ersten / letzten Scheitelpunkt einer solchen Linie mit dem nächsten / vorherigen Scheitelpunkt (Koordinatendifferenz) vergleichen und den Winkel berechnen können. Sie haben also einen Winkel Ihrer Linie im Schnittpunkt. Hier gibt es ein Problem: Sie müssen dieses Tool mehrmals manuell ausführen, um zu erkennen, wie die Ausgabe geschrieben wird. Ich meine, wenn es eine Polylinie braucht, teile sie, schreibe zwei Ergebnispolylinien in die Ausgabe und fahre dann mit der nächsten Polylinie fort und wiederhole sie. Oder kann dieser Teil (Ergebnis der Aufteilung) in verschiedene Speicher-Feature-Classes geschrieben und dann an die Ausgabe angehängt werden. Dies ist das Hauptproblem - zu erkennen, wie die Ausgabe geschrieben wird, um nach dem Teilen nur den ersten Teil jeder Polylinie zu filtern. Eine andere mögliche Lösung besteht darin, alle resultierenden geteilten Polylinien mit zu durchlaufenSuchen Sie den Cursor und nehmen Sie nur den ersten gefundenen (anhand der ID der Quellpolylinien PLINE_FC1).
3. Um den Winkel zu bestimmen, müssen Sie mit arcpy auf die Scheitelpunkte der Ergebnispolylinie zugreifen . Schreibe resultierende Winkel zu Punkten (POINT_FC).
4. Wiederholen Sie die Schritte 2-3 für PLINE_FC2.
5. Subtrahieren Sie die Winkelattribute (in POINT_FC) und erhalten Sie das Ergebnis.

Möglicherweise ist es sehr schwierig, Schritt 2 zu codieren (für einige Tools ist eine ArcInfo-Lizenz erforderlich). Dann können Sie auch versuchen, die Scheitelpunkte jeder Polylinie zu analysieren (nach ID nach Schnittpunkt gruppieren).

Hier ist der Weg, um es zu tun:

1. Nehmen Sie den ersten Schnittpunkt POINT_FC. Holen Sie sich seine Koordinaten ( point_x, point_y)
2. Entnehmen Sie anhand der ID die jeweilige Quellpolylinie aus PLINE_FC1.
3. Nehmen Sie den ersten ( vert0_x, vert0_y) und zweiten ( vert1_x, vert1_y) Scheitelpunkt davon.
4. Berechnen Sie für den ersten Scheitelpunkt den Tangens der Linie zwischen diesem Scheitelpunkt und dem Schnittpunkt: tan0 = (point_y - vert0_y) / (point_x - vert0_x)
5. Berechnen Sie dasselbe für den zweiten Eckpunkt: tan1 = (vert1_y - point_y) / (vert1_x - point_x)
6. Wenn tan1gleich ist tan2, haben Sie zwei Scheitelpunkte Ihrer Linie gefunden, zwischen denen sich ein Schnittpunkt befindet, und Sie können den Schnittwinkel für diese Linie berechnen. Andernfalls müssen Sie mit dem nächsten Vertecepaar (zweites, drittes) fortfahren und so weiter.
7. Wiederholen Sie die Schritte 1 bis 6 für jeden Schnittpunkt.
8. Wiederholen Sie die Schritte 1 bis 7 für die zweite Polylinien-Feature-Class PLINE_FC2.
9. Subtrahieren Sie Winkelattribute von PLINE_FC1 und PLINE_FC2 und erhalten Sie das Ergebnis.

Kürzlich habe ich versucht, es alleine zu machen.

Meine Hinweisfunktion basiert auf kreisförmigen Punkten um die Schnittpunkte von Linien sowie auf Punkten, die sich in einem Abstand von einem Meter von Schnittpunkten befinden. Bei der Ausgabe handelt es sich um eine Polylinien-Feature-Class, die Attribute für die Anzahl der Winkel in Bezug auf Schnittpunkte und Winkel enthält.

Beachten Sie, dass Linien planarisiert werden sollten, um Schnittpunkte zu finden, und dass der Raumbezug mit der korrekten Anzeige der Linienlänge festgelegt werden muss (meine ist WGS_1984_Web_Mercator_Auxiliary_Sphere).

Wird in der ArcMap-Konsole ausgeführt, kann aber problemlos in ein Skript in der Toolbox umgewandelt werden. Dieses Skript verwendet nur die Linienebene im Inhaltsverzeichnis, nicht mehr.

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)