Найти угол между пересекающимися объектами в двух классах объектов, используя ArcGIS Desktop и Python? [закрыто]

У меня есть два пересекающихся линии классов объектов. Я хочу найти угол в каждой точке пересечения, используя ArcGIS 10 и Python.

Кто-нибудь может помочь?

Существует относительно простой рабочий процесс. Это преодолевает потенциальные проблемы, которые две функции могут пересекать в более чем одной точке. Это не требует сценариев (но может быть легко превращено в сценарий). Это можно сделать в первую очередь из меню ArcGIS.

Идея состоит в том, чтобы использовать слой точек пересечения, по одной точке для каждой отдельной пары пересекающихся полилиний. Вам нужно получить небольшой кусочек каждой пересекающейся ломаной в этих точках пересечения. Используйте ориентации этих частей, чтобы вычислить их углы пересечения.

Вот шаги:

1. Убедитесь, что каждый из объектов полилинии имеет уникальный идентификатор в своей таблице атрибутов. Это будет использовано позже, чтобы соединить некоторые геометрические атрибуты полилиний с таблицей точек пересечения.

2. Геообработка | Пересечение получает точки (убедитесь, что вы хотите указать точки для вывода).

3. Геообработка | Буфер позволяет вам буферизовать точки на небольшое количество. Сделайте его действительно крошечным, чтобы часть каждой строки в буфере не сгибалась.

4. Геообработка | Клип (применяется дважды) ограничивает исходные слои полилинии только буферами. Поскольку это создает новые наборы данных для вывода, последующие операции не изменят исходные данные (что хорошо).

   

   Вот схема того, что происходит: два полилинейных слоя, показанные голубым и светло-красным, имеют темные точки пересечения. Вокруг этих точек крошечные буферы показаны желтым цветом. Темно-синий и красный сегменты показывают результаты обрезки оригинальных объектов в эти буферы. Остальная часть алгоритма работает с темными сегментами. (Вы не можете видеть это здесь, но крошечная красная ломаная пересекает две синие линии в общей точке, создавая то, что кажется буфером вокруг двух синих полилиний. Это действительно два буфера вокруг двух перекрывающихся точек пересечения красно-синего цвета. Таким образом, эта диаграмма отображает пять буферов всего.)

5. Используйте инструмент AddField , чтобы создать четыре новых поля в каждом из этих обрезанных слоев: [X0], [Y0], [X1] и [Y1]. Они будут содержать точечные координаты, поэтому сделайте их двойными и придайте им большую точность.

6. Вычислить геометрию (вызывается щелчком правой кнопкой мыши на каждом новом заголовке поля) позволяет вычислять координаты x и y начальной и конечной точек каждой обрезанной полилинии: поместите их в [X0], [Y0], [X1] и [Y1] соответственно. Это делается для каждого обрезанного слоя, поэтому необходимо 8 вычислений.

7. Используйте инструмент AddField , чтобы создать новое поле [Angle] в слое точек пересечения.

8. Соедините вырезанные таблицы с таблицей точек пересечения на основе общих идентификаторов объектов. (Соединения выполняются путем щелчка правой кнопкой мыши на имени слоя и выбора «Объединения и взаимосвязи».)

   На данный момент таблица пересечения точек имеет 9 новых полей: два с именем [X0] и т. Д., А одно с именем [Угол]. Псевдоним полей [X0], [Y0], [X1] и [Y1], которые принадлежат одной из соединяемых таблиц. Давайте назовем их (скажем) «X0a», «Y0a», «X1a» и «Y1a».

9. Используйте Калькулятор поля, чтобы вычислить угол в таблице пересечений. Вот блок кода Python для расчета:

   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
   

   Выражение вычисления поля, конечно, просто

   c

Несмотря на длину этого кодового блока, математика проста: (dx, dy) является вектором направления для первой полилинии, а (dxa, dya) является вектором направления для второй. Их длины r и ra (вычисленные по теореме Пифагора) используются для нормализации их к единичным векторам. (Не должно быть никаких проблем с нулевой длиной, потому что отсечение должно давать характеристики положительной длины.) Размер их клинового произведения dx dya - dydxa (после деления на r и ra) - синус угла. (Использование продукта клина, а не обычного внутреннего продукта должно обеспечить лучшую числовую точность для почти нулевых углов.) Наконец, угол преобразуется из радианов в градусы. Результат будет лежать в диапазоне от 0 до 90. Обратите внимание на избегание тригонометрии до самого конца: этот подход имеет тенденцию давать надежные и легко вычисляемые результаты.

Некоторые точки могут появляться несколько раз в слое пересечения. Если это так, они получат несколько углов, связанных с ними.

Буферизация и отсечение в этом решении относительно дороги (шаги 3 и 4): вы не хотите делать это таким образом, когда задействованы миллионы точек пересечения. Я рекомендовал его, потому что (а) он упрощает процесс поиска двух последовательных точек вдоль каждой ломаной линии в окрестности ее точки пересечения и (б) буферизация настолько проста, что ее легко выполнить в любой ГИС - дополнительное лицензирование не требуется выше базового уровня ArcMap - и обычно дает правильные результаты. (Другие операции «геообработки» могут быть не такими надежными.)

Я считаю, что вам нужно создать скрипт Python.

Вы можете сделать это, используя инструменты геообработки и arcpy.

Вот основные инструменты и идеи, которые могут быть полезны для вас:

1. Сделайте пересечение двух ваших полилиний (давайте назовем их PLINE_FC1, PLINE_FC2) классов объектов (в результате вам понадобятся точечные объекты - POINT_FC), используя инструмент Intersect . У вас будут идентификаторы из PLINE_FC1, PLINE_FC2 в точках POINT_FC.
2. Разделить PLINE_FC1 по POINT_FC, используя инструмент Split Line At Point, В результате вы получите разделенные полилинии - главное преимущество в том, что вы можете просто взять первую / последнюю вершину такой линии, сравнить ее со следующей / предыдущей вершиной (разностью координат) и вычислить угол. Таким образом, у вас будет угол вашей линии в точке пересечения. Здесь есть одна проблема - вам нужно запустить этот инструмент вручную несколько раз, чтобы понять, как записывается вывод. Я имею в виду, что если взять полилинию, разделить ее, записать две результирующие полилинии в вывод, а затем перейти к следующей полилинии и повторить. Или, может быть, эта часть (результат разбиения) записывается в различные классы пространственных объектов памяти, а затем добавляется к выводу. Это главная проблема - понять, как пишется вывод, чтобы иметь возможность фильтровать только первую часть каждой полилинии после разбиения. Другое возможное решение состоит в том, чтобы соединить все результирующие расщепленные полилинии сSearchCursor и взять только первый встреченный (по ID исходных полилиний PLINE_FC1).
3. Чтобы получить угол, вам нужно получить доступ к вершинам результирующей полилинии, используя arcpy . Запишите результирующие углы в точки (POINT_FC).
4. Повторите шаги 2-3 для PLINE_FC2.
5. Вычтите атрибуты угла (в POINT_FC) и получите результат.

Возможно, будет очень сложно кодировать шаг 2 (также для некоторых инструментов требуется лицензия ArcInfo). Затем вы также можете попытаться проанализировать вершины каждой ломаной линии (сгруппировав их по ID после пересечения).

Вот способ сделать это:

1. Возьмите первую точку пересечения POINT_FC. Получить его координаты ( point_x, point_y)
2. По его идентификатору возьмите соответствующую исходную полилинию из PLINE_FC1.
3. Возьмите первый ( vert0_x, vert0_y) и второй ( vert1_x, vert1_y) его стихи.
4. Для первой вершины вычислите тангенс линии между этой вершиной и точкой пересечения: tan0 = (point_y - vert0_y) / (point_x - vert0_x)
5. Вычислите то же самое для второй вершины: tan1 = (vert1_y - point_y) / (vert1_x - point_x)
6. Если tan1равно tan2, то вы нашли две вершины вашей линии, между которыми есть точка пересечения, и вы можете рассчитать угол пересечения для этой линии. В противном случае вам нужно перейти к следующей паре стихов (второй, третий) и так далее.
7. Повторите шаги 1-6 для каждой точки пересечения.
8. Повторите шаги 1-7 для второго класса объектов полилинии PLINE_FC2.
9. Вычтите атрибуты угла из PLINE_FC1 и PLINE_FC2 и получите результат.

Недавно я пытался сделать это самостоятельно.

Моя подсказка основана на круглых точках вокруг пересечений линий, а также на точках, расположенных на расстоянии одного метра от пересечений. Выходными данными является класс пространственных объектов ломаной линии, который имеет атрибуты числа углов на пересечениях и углах.

Обратите внимание, что линии должны быть выровнены, чтобы найти пересечения, и пространственная привязка должна быть установлена ​​с отображением правильной длины линии (у меня WGS_1984_Web_Mercator_Auxili_Sphere).

Работает в консоли ArcMap, но легко может быть превращен в скрипт на панели инструментов. Этот скрипт использует только линейный слой в оглавлении, не более того.

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)