Mencari sudut antara fitur berpotongan di dua kelas fitur menggunakan ArcGIS Desktop dan Python? [Tutup]

Saya memiliki dua kacamata fitur garis berpotongan. Saya ingin menemukan sudut pada setiap titik persimpangan menggunakan ArcGIS 10 dan Python.

Adakah yang bisa membantu?

Ada alur kerja yang relatif sederhana. Ini mengatasi masalah potensial bahwa dua fitur dapat berpotongan di lebih dari satu titik. Itu tidak memerlukan skrip (tetapi dapat dengan mudah diubah menjadi skrip). Ini dapat dilakukan terutama dari menu ArcGIS.

Idenya adalah untuk mengeksploitasi lapisan titik persimpangan, satu titik untuk setiap pasangan berpotongan yang berbeda. Anda perlu mendapatkan bagian kecil dari setiap polyline berpotongan di titik-titik persimpangan ini. Gunakan orientasi potongan-potongan ini untuk menghitung sudut persimpangan mereka.

Berikut langkah-langkahnya:

1. Pastikan setiap fitur polyline memiliki pengidentifikasi unik di dalam tabel atributnya. Ini akan digunakan nanti untuk menggabungkan beberapa atribut geometris dari polyline ke tabel titik persimpangan.

2. Geoprocessing | Intersect mendapatkan poin (pastikan untuk menentukan poin yang Anda inginkan untuk output).

3. Geoprocessing | Buffer memungkinkan Anda melakukan buffer poin dengan jumlah kecil. Buat benar - benar kecil agar porsi setiap garis dalam buffer tidak menekuk.

4. Geoprocessing | Klip (diterapkan dua kali) membatasi lapisan polyline asli hanya pada buffer. Karena ini menghasilkan set data baru untuk outputnya, operasi selanjutnya tidak akan mengubah data asli (yang merupakan hal yang baik).

   

   Berikut adalah skema tentang apa yang terjadi: dua lapisan polyline, ditunjukkan dengan warna biru muda dan merah muda, telah menghasilkan titik persimpangan gelap. Di sekitar titik-titik itu buffer kecil ditampilkan dengan warna kuning. Segmen biru dan merah yang lebih gelap menunjukkan hasil kliping fitur asli ke buffer ini. Algoritma lainnya bekerja dengan segmen gelap. (Anda tidak dapat melihatnya di sini, tetapi sebuah polyline kecil merah memotong dua garis biru pada titik yang sama, menghasilkan apa yang tampaknya menjadi penyangga di sekitar dua polyline biru. Ini benar-benar dua buffer di sekitar dua titik yang tumpang tindih dari persimpangan merah-biru. , diagram ini menampilkan lima buffer secara keseluruhan.)

5. Gunakan alat AddField untuk membuat empat bidang baru di masing-masing lapisan terpotong ini: [X0], [Y0], [X1], dan [Y1]. Mereka akan memegang titik koordinat, jadi buatlah mereka berlipat ganda dan berikan mereka banyak ketepatan.

6. Calculate Geometry (dipanggil dengan mengklik kanan pada setiap header bidang baru) memungkinkan Anda untuk menghitung koordinat x- dan y- dari titik awal dan akhir dari setiap polyline yang terpotong: masukkan ini ke dalam [X0], [Y0], [X1] , dan [Y1], masing-masing. Ini dilakukan untuk setiap layer yang terpotong, sehingga diperlukan 8 perhitungan.

7. Gunakan alat AddField untuk membuat bidang [Angle] baru di lapisan titik persimpangan.

8. Gabungkan tabel terpotong ke tabel titik persimpangan berdasarkan pengidentifikasi objek umum. (Bergabung dilakukan dengan mengklik kanan pada nama layer dan memilih "Bergabung dan Terkait".)

   Pada titik ini tabel persimpangan titik memiliki 9 bidang baru: dua dinamai [X0], dll., Dan satu dinamai [Sudut]. Alias bidang [X0], [Y0], [X1], dan [Y1] yang termasuk dalam salah satu tabel bergabung. Sebut ini (katakanlah) "X0a", "Y0a", "X1a", dan "Y1a".

9. Gunakan Kalkulator Bidang untuk menghitung sudut pada tabel persimpangan. Berikut ini adalah blok kode Python untuk perhitungan:

   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

   Ekspresi perhitungan lapangan, tentu saja, hanyalah

   c

Meskipun panjang blok kode ini, matematika sederhana: (dx, dy) adalah vektor arah untuk polyline pertama dan (dxa, dya) adalah vektor arah untuk yang kedua. Panjangnya, r dan ra (dihitung melalui Teorema Pythagoras), digunakan untuk menormalkannya ke satuan vektor. (Seharusnya tidak ada masalah dengan panjang nol, karena kliping harus menghasilkan fitur panjang positif.) Ukuran produk wedge mereka dx dya - dydxa (setelah pembagian dengan r dan ra) adalah sinus sudut. (Menggunakan produk baji daripada produk dalam biasa harus memberikan presisi numerik yang lebih baik untuk sudut mendekati nol.) Akhirnya, sudut diubah dari radian ke derajat. Hasilnya akan berada di antara 0 dan 90. Perhatikan penghindaran trigonometri sampai akhir: pendekatan ini cenderung menghasilkan hasil yang andal dan mudah dihitung.

Beberapa titik mungkin muncul beberapa kali di lapisan persimpangan. Jika demikian, mereka akan mendapatkan banyak sudut yang terkait dengannya.

Menyangga dan memotong dalam solusi ini relatif mahal (langkah 3 dan 4): Anda tidak ingin melakukannya dengan cara ini ketika jutaan titik persimpangan dilibatkan. Saya merekomendasikannya karena (a) menyederhanakan proses menemukan dua titik berturut-turut di sepanjang setiap polyline dalam lingkungan titik persimpangannya dan (b) buffering sangat mendasar sehingga mudah dilakukan dalam GIS apa pun - tidak diperlukan lisensi tambahan di atas level ArcMap dasar - dan biasanya menghasilkan hasil yang benar. (Operasi "geoprosesing" lainnya mungkin tidak dapat diandalkan.)

Saya yakin Anda perlu membuat skrip python.

Anda dapat melakukannya menggunakan alat geoprocessing dan arcpy.

Berikut adalah alat dan ide utama yang dapat berguna bagi Anda:

1. Buat persimpangan dua polyline Anda (sebut saja mereka PLINE_FC1, PLINE_FC2) dengan kacamata (Anda memerlukan fitur titik sebagai hasilnya - POINT_FC) menggunakan tool Intersect . Anda akan memiliki ID dari PLINE_FC1, PLINE_FC2 dalam poin POINT_FC.
2. Split PLINE_FC1 oleh POINT_FC menggunakan alat Split Line At Point. Dengan hasil, Anda akan memiliki polyline yang terpecah - keuntungan utamanya adalah Anda hanya dapat mengambil simpul pertama / terakhir dari garis tersebut membandingkannya dengan simpul berikutnya / sebelumnya (perbedaan koordinat) dan menghitung sudut. Jadi, Anda akan memiliki sudut garis Anda di titik berpotongan. Ada satu masalah di sini - Anda harus menjalankan alat ini secara manual beberapa kali untuk menyadari bagaimana output ditulis. Maksud saya jika membutuhkan polyline, membaginya, menulis dua polyline hasil ke output dan kemudian melanjutkan ke polyline berikutnya dan ulangi. Atau mungkin bagian ini (hasil pemisahan) ditulis ke kelas fitur memori yang berbeda, kemudian ditambahkan ke output. Ini adalah masalah utama - untuk menyadari bagaimana output ditulis untuk dapat memfilter hanya bagian pertama dari setiap polyline setelah pemisahan. Solusi lain yang mungkin adalah loop melalui semua hasil polylines denganCari kursor dan ambil hanya pertama kali ditemukan (berdasarkan ID dari sumber polylines PLINE_FC1)
3. Untuk mendapatkan sudut, Anda harus mengakses verteces hasil polyline menggunakan arcpy . Tulis sudut yang dihasilkan ke poin (POINT_FC).
4. Ulangi langkah 2-3 untuk PLINE_FC2.
5. Kurangi atribut sudut (di POINT_FC) dan dapatkan hasilnya.

Mungkin akan sangat sulit untuk kode langkah 2 (juga beberapa alat memerlukan lisensi ArcInfo). Kemudian Anda juga dapat mencoba menganalisis vertek setiap polyline (mengelompokkannya berdasarkan ID setelah persimpangan).

Inilah cara untuk melakukannya:

1. Ambil titik persimpangan pertama POINT_FC. Dapatkan koordinatnya ( point_x, point_y)
2. Dengan ID-nya, ambil polyline masing-masing sumber dari PLINE_FC1.
3. Ambil vertek pertama ( vert0_x, vert0_y) dan kedua ( vert1_x, vert1_y).
4. Untuk titik pertama, hitung garis singgung garis antara titik dan titik persimpangan ini: tan0 = (point_y - vert0_y) / (point_x - vert0_x)
5. Hitung hal yang sama untuk simpul kedua: tan1 = (vert1_y - point_y) / (vert1_x - point_x)
6. Jika tan1sama dengan tan2, maka Anda telah menemukan dua vertek dari garis Anda yang memiliki titik persimpangan di antaranya dan Anda dapat menghitung sudut persimpangan untuk garis ini. Kalau tidak, Anda harus melanjutkan ke pasangan vertek berikutnya (kedua, ketiga) dan seterusnya.
7. Ulangi langkah 1-6 untuk setiap titik persimpangan.
8. Ulangi langkah 1-7 untuk kelas fitur polyline kedua PLINE_FC2.
9. Kurangi atribut sudut dari PLINE_FC1 dan PLINE_FC2 dan dapatkan hasilnya.

Baru-baru ini saya mencoba melakukannya sendiri.

Fitur petunjuk saya didasarkan pada titik melingkar di sekitar persimpangan garis serta titik yang terletak pada jarak satu meter dari persimpangan. Outputnya adalah kelas fitur polyline yang memiliki atribut nomor sudut pada persimpangan dan sudut.

Perhatikan bahwa garis harus direncanakan untuk menemukan persimpangan dan referensi spasial harus ditetapkan dengan tampilan panjang garis yang benar (tambang adalah WGS_1984_Web_Mercator_Auxiliary_Sphere).

Berjalan di konsol ArcMap tetapi dengan mudah dapat diubah menjadi skrip di kotak alat. Script ini hanya menggunakan layer line di TOC, tidak lebih.

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)