태그 보관물: arcgis-10.0

arcgis-10.0

ArcGIS Desktop과 Python을 사용하여 두 피쳐 클래스에서 교차 피쳐 사이의 각도를 찾으십니까? [닫은]

두 개의 교차 선 피쳐 클래스가 있습니다. ArcGIS 10과 Python을 사용하여 각 교차점에서 각도를 찾고 싶습니다.

누구든지 도울 수 있습니까?



답변

비교적 간단한 워크 플로우가 있습니다. 두 기능이 하나 이상의 지점에서 교차 할 수있는 잠재적 인 문제를 극복합니다. 스크립팅이 필요하지는 않지만 스크립트로 쉽게 전환 할 수 있습니다. 주로 ArcGIS 메뉴에서 수행 할 수 있습니다.

아이디어는 교차하는 폴리 라인의 각 개별 쌍에 대해 하나씩, 교차점 레이어를 활용하는 것입니다. 이 교차점에서 교차하는 각 폴리 라인 의 작은 조각 을 가져와야합니다 . 이 조각들의 방향을 사용하여 교차 각을 계산하십시오.

단계는 다음과 같습니다.

  1. 각 폴리 라인 피처 의 속성 테이블 내에 고유 식별자 가 있는지 확인하십시오 . 나중에 폴리 라인의 일부 기하학적 속성을 교차점 테이블에 결합하는 데 사용됩니다.

  2. Geoprocessing | Intersect 는 점을 얻습니다 ( 출력에 대한 을 지정해야 함 ).

  3. Geoprocessing | Buffer를 사용하면 포인트를 소량 버퍼링 할 수 있습니다. 버퍼 내의 각 라인 부분이 구부러지지 않도록 정말 작게 만드십시오 .

  4. Geoprocessing | Clip (두 번 적용)은 원래 폴리 라인 레이어를 버퍼로만 제한합니다. 이렇게하면 출력에 대한 새 데이터 세트가 생성되므로 후속 작업에서는 원본 데이터를 수정하지 않습니다 (좋은 일).

    그림

    밝은 파란색과 밝은 빨간색으로 표시된 두 개의 폴리 라인 레이어가 어두운 교차점을 생성했습니다. 이 지점 주위에 작은 버퍼가 노란색으로 표시됩니다. 진한 파란색과 빨간색 선분은 원래 기능을이 버퍼로 클리핑 한 결과를 보여줍니다. 나머지 알고리즘은 어두운 세그먼트에서 작동합니다. (여기서는 볼 수 없지만 작은 빨간색 폴리 라인은 공통점에서 두 개의 파란색 선을 교차하여 두 개의 파란색 폴리 라인 주위에 버퍼로 보이는 것을 생성합니다. 이는 실제로 겹치는 두 개의 빨간색-파란색 교차점 주위에 두 개의 버퍼입니다. 이 다이어그램은 5 개의 버퍼를 모두 표시합니다.)

  5. AddField 도구를 사용하여 클리핑 된 각 레이어에 [X0], [Y0], [X1] 및 [Y1]의 새 필드 4 개를 만듭니다. 점 좌표를 유지하므로 두 배로 만들고 많은 정밀도를 부여합니다.

  6. 형상 계산 (새 필드 헤더 각각을 마우스 오른쪽 버튼으로 클릭하여 호출)을 사용하면 클리핑 된 각 폴리 라인의 시작점과 끝점의 x 및 y 좌표를 계산할 수 있습니다. [X0], [Y0], [X1] 및 [Y1]. 이것은 클리핑 된 각 레이어에 대해 수행되므로 8 개의 계산이 필요합니다.

  7. 사용 AddField의 교차점 층에 새로운 [각도] 필드를 생성하는 도구를.

  8. 공통 객체 식별자를 기반으로 클리핑 된 테이블을 교차점 테이블에 조인 합니다. 레이어 이름을 마우스 오른쪽 버튼으로 클릭하고 “조인 및 관련”을 선택하여 조인을 수행합니다.

    이 시점에서 포인트 교차점 테이블에는 9 개의 새로운 필드가 있습니다. 두 개는 [X0] 등이고 다른 하나는 [Angle]입니다. 결합 된 테이블 중 하나에 속하는 [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 (피타고라스 정리를 통해 계산 된)는 단위 벡터로 정규화하는 데 사용됩니다. (클리핑은 양의 길이의 특징을 생성해야하므로 길이가 0 인 경우 문제가되지 않습니다.) 웨지 제품의 크기 dx dya-dydxa (r과 ra로 나눈 후)는 각도의 사인입니다. (일반적인 내부 제품이 아닌 웨지 제품을 사용하면 거의 0에 가까운 각도에 대해 더 나은 수치 정밀도를 제공해야합니다.) 마지막으로 각도는 라디안에서 각도로 변환됩니다. 결과는 0에서 90 사이입니다. 끝날 때까지 삼각법을 피하십시오.이 방법은 신뢰할 수 있고 쉽게 계산 된 결과를 생성하는 경향이 있습니다.

교차점에 일부 점이 여러 번 나타날 수 있습니다. 그렇다면 여러 각도를 얻게됩니다.

이 솔루션의 버퍼링 및 클리핑은 상대적으로 비용이 많이 듭니다 (3 단계 및 4 단계). (a) 교차점 부근의 각 폴리 라인을 따라 두 개의 연속 점을 찾는 프로세스를 단순화하고 (b) 버퍼링이 기본이므로 모든 GIS에서 쉽게 수행 할 수 있으므로 추가 라이센스가 필요하지 않기 때문에 권장했습니다. 기본 ArcMap 수준 이상이며 일반적으로 올바른 결과를 생성합니다. 다른 “지오 프로세싱”작업은 그다지 신뢰할 수 없습니다.


답변

파이썬 스크립트를 만들어야한다고 생각합니다.

지오 프로세싱 도구와 arcpy를 사용하여 수행 할 수 있습니다.

유용한 도구와 아이디어는 다음과 같습니다.

  1. 도구 교차를 사용하여 두 폴리 라인 (PLINE_FC1, PLINE_FC2라고 함) 피처 클래스 (포인트 피처가 필요함-POINT_FC)를 교차 시킵니다. 포인트 POINT_FC에 PLINE_FC1, PLINE_FC2의 ID가 있습니다.
  2. 스플릿 라인 포인트를 사용하여 PLINE_FC1을 POINT_FC로 분할. 결과적으로 폴리 라인이 분할됩니다. 폴리 라인의 주요 이점은 해당 라인의 첫 번째 / 마지막 정점을 다음 / 이전 버텍스 (좌표 차이)와 비교하고 각도를 계산할 수 있다는 것입니다. 따라서 교차점에 선 각도가 생깁니다. 여기에는 한 가지 문제가 있습니다. 출력을 작성하는 방법을 이해하려면이 도구를 수동으로 여러 번 실행해야합니다. 폴리 라인이 필요하다면, 분리하고, 결과에 2 개의 결과 폴리 라인을 쓴 후 다음 폴리 라인으로 진행하여 반복하십시오. 또는이 부분 (분할 결과)이 다른 메모리 피처 클래스에 기록 된 다음 출력에 추가 될 수 있습니다. 이것은 주요 문제입니다. 분할 후 각 폴리 라인의 첫 부분 만 필터링 할 수 있도록 출력을 작성하는 방법을 깨닫는 것입니다. 또 다른 가능한 해결책은 모든 결과 분할 폴리 라인을SearchCursor를 처음 발견합니다 (소스 폴리 라인 PLINE_FC1의 ID 별).
  3. 각도를 얻기 위해 당신이 사용하는 액세스 결과 폴리 라인의 verteces 필요합니다 arcpy을 . 결과 각도를 점에 씁니다 (POINT_FC).
  4. PLINE_FC2에 대해 2-3 단계를 반복하십시오.
  5. POINT_FC의 각도 속성을 빼고 결과를 얻습니다.

2 단계를 코딩하는 것은 매우 어려울 수 있습니다 (일부 도구에는 ArcInfo 라이센스가 필요함). 그런 다음 모든 폴리선의 정점을 분석 할 수 있습니다 (교차 후 ID별로 그룹화).

방법은 다음과 같습니다.

  1. 첫 번째 교차점 POINT_FC를 가져갑니다. 좌표 구하기 ( point_x, point_y)
  2. ID별로 PLINE_FC1에서 각 소스 폴리 라인을 가져옵니다.
  3. (첫번째 가지고 vert0_x, vert0_y)와 초 ( vert1_x, vert1_y그것의) verteces.
  4. 첫 번째 정점의 경우이 정점과 교차점 사이의 선의 접선을 계산하십시오. tan0 = (point_y - vert0_y) / (point_x - vert0_x)
  5. 두 번째 정점에 대해 동일한 것을 계산하십시오. tan1 = (vert1_y - point_y) / (vert1_x - point_x)
  6. 경우 tan1IS가 동일 tan2, 당신은 사이에 교차점이 당신의 라인이 verteces을 발견하고이 라인의 교차 각도를 계산할 수 있습니다. 그렇지 않으면 다음 정점 쌍 (두 번째, 세 번째) 등으로 진행해야합니다.
  7. 모든 교차점에 대해 1-6 단계를 반복하십시오.
  8. 두 번째 폴리 라인 피쳐 클래스 PLINE_FC2에 대해 1 ~ 7 단계를 반복합니다.
  9. PLINE_FC1 및 PLINE_FC2에서 각도 속성을 빼고 결과를 얻습니다.

답변

최근에 나는 그것을 스스로하려고했습니다.

내 단서 기능은 선의 교차점 주위의 원형 점과 교차점에서 1 미터 거리에있는 점을 기반으로합니다. 출력은 교차점 및 각도에 각도 수 속성이있는 폴리 라인 피쳐 클래스입니다.

교차점을 찾으려면 선을 평면화해야하며 공간 참조가 올바른 선 길이 표시로 설정되어야합니다 (광산은 WGS_1984_Web_Mercator_Auxiliary_Sphere 임).

ArcMap 콘솔에서 실행되지만 툴박스에서 스크립트로 쉽게 전환 할 수 있습니다. 이 스크립트는 TOC에서 라인 레이어 만 사용합니다.

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)


답변