Secara efisien menemukan tetangga urutan pertama dari 200 ribu poligon

Untuk setiap kelompok blok Sensus 208.781, saya ingin mengambil ID FIPS dari semua tetangga urutan pertama. Saya memiliki semua batas TIGER yang diunduh dan digabung menjadi 1GB shapefile tunggal.

Saya mencoba skrip ArcPython yang menggunakan SelectLayerByLocation untuk BOUNDARY_TOUCHES pada intinya, tetapi dibutuhkan lebih dari 1 detik untuk setiap grup blok yang lebih lambat dari yang saya inginkan. Ini bahkan setelah saya membatasi pencarian SelectLayerByLocation untuk memblokir grup dalam kondisi yang sama. Saya menemukan skrip ini , tetapi juga menggunakan SelectLayerByLocation secara internal sehingga tidak lebih cepat.

Solusinya tidak harus berbasis Arc - saya terbuka untuk paket lain, meskipun saya paling nyaman coding dengan Python.

Jika Anda memiliki akses ke ArcGIS 10.2 untuk Desktop, atau mungkin sebelumnya, maka saya pikir alat Polygon Neighbors (Analisis) yang:

Membuat tabel dengan statistik berdasarkan pada persentuhan poligon (tumpang tindih, tepi bertepatan, atau node).

mungkin membuat tugas ini lebih mudah sekarang.

Untuk solusi menghindari ArcGIS, gunakan pysal . Anda bisa mendapatkan bobot langsung dari shapefile menggunakan:

w = pysal.rook_from_shapefile("../pysal/examples/columbus.shp")

atau

w = pysal.queen_from_shapefile("../pysal/examples/columbus.shp")

Buka dokumen untuk info lebih lanjut.

Hanya pembaruan. Setelah mengikuti saran Whuber, saya menemukan bahwa Generate Spatial Weights Matrix hanya menggunakan loop Python dan kamus untuk menentukan tetangga. Saya mereproduksi proses di bawah ini.

Bagian pertama loop melalui setiap simpul dari setiap grup blok. Ini membuat kamus dengan koordinat titik sebagai kunci dan daftar ID grup blok yang memiliki titik pada koordinat itu sebagai nilai. Perhatikan bahwa ini membutuhkan dataset yang rapi secara topologi, karena hanya vertex / vertex yang tumpang tindih yang akan terdaftar sebagai hubungan tetangga. Untungnya, pembentukan blok grup TIGER Biro Sensus tidak ada masalah dalam hal ini.

Bagian kedua loop melalui setiap simpul dari setiap kelompok blok lagi. Ini membuat kamus dengan ID grup blok sebagai kunci dan ID tetangga kelompok itu sebagai nilainya.

# Create dictionary of vertex coordinate : [...,IDs,...]
BlockGroupVertexDictionary = {}
BlockGroupCursor = arcpy.SearchCursor(BlockGroups.shp)
BlockGroupDescription = arcpy.Describe(BlockGroups.shp)
BlockGroupShapeFieldName = BlockGroupsDescription.ShapeFieldName
#For every block group...
for BlockGroupItem in BlockGroupCursor :
    BlockGroupID = BlockGroupItem.getValue("BKGPIDFP00")
    BlockGroupFeature = BlockGroupItem.getValue(BlockGroupShapeFieldName)
    for BlockGroupPart in BlockGroupFeature:
        #For every vertex...
        for BlockGroupPoint in BlockGroupPart:
            #If it exists (and isnt empty interior hole signifier)...
            if BlockGroupPoint:
                #Create string version of coordinate
                PointText = str(BlockGroupPoint.X)+str(BlockGroupPoint.Y)
                #If coordinate is already in dictionary, append this BG's ID
                if PointText in BlockGroupVertexDictionary:
                    BlockGroupVertexDictionary[PointText].append(BlockGroupID)
                #If coordinate is not already in dictionary, create new list with this BG's ID
                else:
                    BlockGroupVertexDictionary[PointText] = [BlockGroupID]
del BlockGroupItem
del BlockGroupCursor


#Create dictionary of ID : [...,neighbors,...]
BlockGroupNeighborDictionary = {}
BlockGroupCursor = arcpy.SearchCursor(BlockGroups.shp)
BlockGroupDescription = arcpy.Describe(BlockGroups.shp)
BlockGroupShapeFieldName = BlockGroupDescription.ShapeFieldName
#For every block group
for BlockGroupItem in BlockGroupCursor:
    ListOfBlockGroupNeighbors = []
    BlockGroupID = BlockGroupItem.getValue("BKGPIDFP00")
    BlockGroupFeature = BlockGroupItem.getValue(BlockGroupShapeFieldName)
    for BlockGroupPart in BlockGroupFeature:
        #For every vertex
        for BlockGroupPoint in BlockGroupPart:
            #If it exists (and isnt interior hole signifier)...
            if BlockGroupPoint:
                #Create string version of coordinate
                PointText = str(BlockGroupPoint.X)+str(BlockGroupPoint.Y)
                if PointText in BlockGroupVertexDictionary:
                    #Get list of block groups that have this point as a vertex
                    NeighborIDList = BlockGroupVertexDictionary[PointText]
                    for NeighborID in NeighborIDList:
                        #Don't add if this BG already in list of neighbors
                        if NeighborID in ListOfBGNeighbors:
                            pass
                        #Add to list of neighbors (as long as its not itself)
                        elif NeighborID != BlockGroupID:
                            ListOfBGNeighbors.append(NeighborID)
    #Store list of neighbors in blockgroup object in dictionary
    BlockGroupNeighborDictionary[BlockGroupID] = ListOfBGNeighbors

del BlockGroupItem
del BlockGroupCursor
del BlockGroupVertexDictionary

Di belakang saya menyadari saya bisa menggunakan metode yang berbeda untuk bagian kedua yang tidak memerlukan perulangan melalui shapefile lagi. Tapi ini yang saya gunakan, dan berfungsi cukup baik bahkan untuk 1000-an grup blok sekaligus. Saya belum mencoba melakukannya dengan seluruh AS, tetapi dapat mengeksekusi untuk seluruh negara.

Alternatifnya mungkin menggunakan PostgreSQL dan PostGIS . Saya telah mengajukan beberapa pertanyaan tentang cara melakukan perhitungan serupa di situs ini:

• Perhitungan Tetangga Terdekat di PostGIS?
• Identifikasi hubungan topologi menggunakan PostGIS

Saya menemukan ada kurva belajar yang curam untuk mencari tahu bagaimana berbagai bagian dari perangkat lunak cocok bersama, tetapi saya menemukan itu bagus untuk melakukan perhitungan pada lapisan vektor besar. Saya sudah menjalankan beberapa perhitungan tetangga terdekat pada jutaan poligon dan itu cepat dibandingkan dengan ArcGIS.

Hanya beberapa komentar ... metode esri / ArcGIS saat ini menggunakan kamus untuk menyimpan informasi tetapi perhitungan inti dilakukan dalam C ++ menggunakan Polygon Neighbors Tool. Alat ini menghasilkan Tabel yang berisi informasi kedekatan serta attr opsional seperti panjang batas bersama. Anda dapat menggunakan Generate Spatial Weights Matrix Tool jika Anda ingin menyimpan dan kemudian menggunakan kembali info berulang-ulang. Anda juga dapat menggunakan fungsi ini di WeightsUtilities untuk menghasilkan kamus [akses acak] dengan info kedekatan:

contDict = polygonNeighborDict(inputFC, masterField, contiguityType = "ROOK")

di mana inputFC = semua jenis kelas fitur poligon, masterField adalah bidang "ID unik" dari bilangan bulat dan contiguityType di {"ROOK", "QUEEN"}.

Ada upaya di esri untuk melewati aspek tabular untuk pengguna Python dan langsung ke iterator yang akan membuat banyak kasus penggunaan jauh lebih cepat. Paket PySAL dan spdep dalam R adalah alternatif yang fantastis [lihat jawaban radek] . Saya pikir Anda diharuskan untuk menggunakan shapefile sebagai format data dalam paket ini yang selaras dengan format input utas ini. Tidak yakin bagaimana mereka menangani poligon yang tumpang tindih serta poligon dalam poligon. Hasilkan SWM serta fungsi yang saya jelaskan akan menghitung hubungan spasial tersebut sebagai Tetangga "ROOK" DAN "QUEEN".