Trouver l'angle entre les entités qui se croisent dans deux classes d'entités à l'aide d'ArcGIS Desktop et de Python? [fermé]

J'ai deux classes de traits d'intersection. Je veux trouver l'angle à chaque point d'intersection à l'aide d'ArcGIS 10 et de Python.

Quelqu'un peut-il aider?

Il existe un flux de travail relativement simple. Il surmonte les problèmes potentiels que deux caractéristiques peuvent recouper en plusieurs points. Il ne nécessite pas de script (mais peut facilement être transformé en script). Cela peut être fait principalement à partir du menu ArcGIS.

L'idée est d'exploiter une couche de points d'intersection, un point pour chaque paire distincte de polylignes qui se croisent. Vous devez obtenir un petit morceau de chaque polyligne intersectée à ces points d'intersection. Utilisez les orientations de ces pièces pour calculer leurs angles d'intersection.

Voici les étapes:

1. Assurez-vous que chacune des entités polylignes possède un identifiant unique dans sa table attributaire. Ceci sera utilisé plus tard pour joindre certains attributs géométriques des polylignes à la table de points d'intersection.

2. Géotraitement | Intersection obtient les points (assurez-vous de spécifier que vous voulez des points pour la sortie).

3. Le géotraitement | Tampon vous permet de tamponner les points d'une petite quantité. Faites-le vraiment minuscule pour que la partie de chaque ligne dans un tampon ne se plie pas.

4. Le géotraitement | Clip (appliqué deux fois) limite les calques de polyligne d'origine aux seuls tampons. Parce que cela produit de nouveaux jeux de données pour sa sortie, les opérations suivantes ne modifieront pas les données d'origine (ce qui est une bonne chose).

   

   Voici un schéma de ce qui se passe: deux couches de polyligne, représentées en bleu clair et en rouge clair, ont produit des points d'intersection sombres. Autour de ces points, de minuscules tampons sont affichés en jaune. Les segments bleu et rouge plus foncés montrent les résultats de la coupure des caractéristiques originales sur ces tampons. Le reste de l'algorithme fonctionne avec les segments sombres. (Vous ne pouvez pas le voir ici, mais une petite polyligne rouge coupe deux des lignes bleues à un point commun, produisant ce qui semble être un tampon autour de deux polylignes bleues. Il s'agit en fait de deux tampons autour de deux points qui se chevauchent d'intersection rouge-bleu. , ce diagramme affiche cinq tampons en tout.)

5. Utilisez l' outil AddField pour créer quatre nouveaux champs dans chacun de ces calques découpés: [X0], [Y0], [X1] et [Y1]. Ils contiendront des coordonnées de points, alors faites-les doubler et donnez-leur beaucoup de précision.

6. Calculer la géométrie (invoqué par un clic droit sur chaque nouvel en-tête de champ) vous permet de calculer les coordonnées x et y des points de début et de fin de chaque polyligne découpée: mettez-les dans [X0], [Y0], [X1] et [Y1], respectivement. Cela se fait pour chaque couche découpée, donc 8 calculs sont nécessaires.

7. Utilisez l' outil AddField pour créer un nouveau champ [Angle] dans la couche de points d'intersection.

8. Joignez les tables tronquées à la table des points d'intersection en fonction des identificateurs d'objet courants. (Les jointures sont effectuées en cliquant avec le bouton droit sur le nom du calque et en sélectionnant "Jointures et relations".)

   À ce stade, la table d'intersection de points a 9 nouveaux champs: deux sont nommés [X0], etc., et un est nommé [Angle]. Alias les champs [X0], [Y0], [X1] et [Y1] qui appartiennent à l'une des tables jointes. Appelons-les (disons) "X0a", "Y0a", "X1a" et "Y1a".

9. Utilisez la calculatrice de champ pour calculer l'angle dans la table d'intersection. Voici un bloc de code Python pour le calcul:

   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

   L'expression de calcul de champ est, bien sûr, simplement

   c

Malgré la longueur de ce bloc de code, le calcul est simple: (dx, dy) est un vecteur de direction pour la première polyligne et (dxa, dya) est un vecteur de direction pour le second. Leurs longueurs, r et ra (calculées via le théorème de Pythagore), sont utilisées pour les normaliser en vecteurs unitaires. (Il ne devrait pas y avoir de problème avec des longueurs nulles, car l'écrêtage devrait produire des caractéristiques de longueur positive.) La taille de leur produit de coin dx dya - dydxa (après division par r et ra) est le sinus de l'angle. (L'utilisation du produit de coin plutôt que du produit interne habituel devrait fournir une meilleure précision numérique pour des angles proches de zéro.) Enfin, l'angle est converti de radians en degrés. Le résultat se situera entre 0 et 90. Notez l'évitement de la trigonométrie jusqu'à la toute fin: cette approche a tendance à produire des résultats fiables et faciles à calculer.

Certains points peuvent apparaître plusieurs fois dans la couche d'intersection. Si c'est le cas, ils obtiendront plusieurs angles qui leur sont associés.

La mise en mémoire tampon et le découpage dans cette solution sont relativement coûteux (étapes 3 et 4): vous ne voulez pas le faire de cette façon lorsque des millions de points d'intersection sont impliqués. Je l'ai recommandé car (a) il simplifie le processus de recherche de deux points successifs le long de chaque polyligne dans le voisinage de son point d'intersection et (b) la mise en mémoire tampon est si basique qu'elle est facile à faire dans n'importe quel SIG - aucune licence supplémentaire n'est nécessaire au-dessus du niveau de base d'ArcMap - et produit généralement des résultats corrects. (D'autres opérations de «géotraitement» peuvent ne pas être aussi fiables.)

Je pense que vous devez créer un script python.

Vous pouvez le faire en utilisant des outils de géotraitement et arcpy.

Voici les principaux outils et idées qui peuvent vous être utiles:

1. Faites l'intersection de vos deux polylignes (appelons-les PLINE_FC1, PLINE_FC2) avec des classes de traits (vous avez donc besoin d'entités ponctuelles - POINT_FC) en utilisant l'outil Intersection . Vous aurez des ID de PLINE_FC1, PLINE_FC2 aux points POINT_FC.
2. Fractionner PLINE_FC1 par POINT_FC à l'aide de l'outil Fractionner la ligne au point. Par résultat, vous aurez des polylignes divisées - le principal avantage de cela est que vous pouvez simplement prendre le premier / dernier sommet de cette ligne le comparer au sommet suivant / précédent (différence de coordonnées) et calculer l'angle. Ainsi, vous aurez l'angle de votre ligne au point d'intersection. Il y a un problème ici - vous devez exécuter cet outil plusieurs fois manuellement pour comprendre comment la sortie est écrite. Je veux dire que s'il prend une polyligne, la diviser, écrire deux polylignes de résultat dans la sortie, puis passer à la polyligne suivante et répéter. Ou peut-être que cette partie (résultat du fractionnement) est écrite dans différentes classes d'entités mémoire, puis ajoutée à la sortie. C'est le problème principal - pour comprendre comment la sortie est écrite pour pouvoir filtrer uniquement la première partie de chaque polyligne après la division. Une autre solution possible consiste à boucler à travers toutes les polylignesSearchCursor et ne prennent que les premiers rencontrés (par ID des polylignes source PLINE_FC1).
3. Afin d'obtenir l'angle, vous devrez accéder aux sommets de la polyligne résultante en utilisant arcpy . Écrivez les angles résultants aux points (POINT_FC).
4. Répétez les étapes 2-3 pour PLINE_FC2.
5. Soustrayez les attributs d'angle (dans POINT_FC) et obtenez le résultat.

Il sera peut-être très difficile de coder l'étape 2 (certains outils nécessitent également une licence ArcInfo). Ensuite, vous pouvez également essayer d'analyser les sommets de chaque polyligne (en les regroupant par ID après intersection).

Voici la façon de procéder:

1. Prenez le premier point d'intersection POINT_FC. Obtenez ses coordonnées ( point_x, point_y)
2. Par son ID, prenez la polyligne source respective de PLINE_FC1.
3. Prenez les premier ( vert0_x, vert0_y) et deuxième ( vert1_x, vert1_y) sommets de celui-ci.
4. Pour le premier sommet, calculer la tangente de la ligne entre ce sommet et le point d'intersection: tan0 = (point_y - vert0_y) / (point_x - vert0_x)
5. Calculez la même chose pour le deuxième sommet: tan1 = (vert1_y - point_y) / (vert1_x - point_x)
6. Si tan1est égal à tan2, alors vous avez trouvé deux sommets de votre ligne qui ont un point d'intersection entre les deux et vous pouvez calculer l'angle d'intersection pour cette ligne. Sinon, vous devez passer à la prochaine paire de sommets (deuxième, troisième) et ainsi de suite.
7. Répétez les étapes 1 à 6 pour chaque point d'intersection.
8. Répétez les étapes 1 à 7 pour la deuxième classe de traits de polyligne PLINE_FC2.
9. Soustrayez les attributs d'angle de PLINE_FC1 et PLINE_FC2 et obtenez le résultat.

Récemment, j'essayais de le faire par moi-même.

Mon indice est basé sur des points circulaires autour des intersections de lignes ainsi que des points situés à un mètre de distance des intersections. La sortie est une classe d'entités polylignes qui a des attributs de nombre d'angles sur les intersections et l'angle.

Notez que les lignes doivent être planarisées afin de trouver les intersections et la référence spatiale doit être définie avec un affichage de longueur de ligne correct (le mien est WGS_1984_Web_Mercator_Auxiliary_Sphere).

Fonctionnant dans la console ArcMap, mais peut facilement être transformé en script dans la boîte à outils. Ce script utilise uniquement la couche de ligne dans la table des matières, rien de plus.

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)