Membandingkan dua geometri di ArcPy?

18

Saya mencoba membandingkan dua kelas fitur terpisah untuk mengidentifikasi perbedaan di antara mereka (semacam fungsi diff). Alur kerja dasar saya:

  1. Saya mengekstrak geometri menggunakan SearchCursor
  2. Simpan geometri dari dua kelas fitur sebagai GeoJSON menggunakan modifikasi __geo_interface__(mendapatkannya dari valveLondon return {'type': 'Polygon', 'coordinates': [[((pt.X, pt.Y) if pt else None) for pt in part] for part in self]} ). Ini untuk menghindari objek geometri bersama yang digunakan ESRI dengan kursor dan ketidakmampuan untuk membuat salinan yang dalam (beberapa diskusi di sini di gis.stackexchange membicarakannya).
  3. Periksa geometri dari dua kelas fitur berdasarkan pengidentifikasi unik. Misalnya, bandingkan geometri FC1 OID1 dengan geometri FC2 OID1. Untuk mendapatkan geometri sebagai instance objek ESRI, panggil arcpy.AsShape()(dimodifikasi untuk membaca poligon berlubang (lihat poin 2 di atas) dengan return cls(Array([map(lambda p: Point(*p) if p is not None else Point(), part) for part in coordinates])). Perbandingannya hanya geom1.equals(geom2)seperti yang ditunjukkan dalam Kelas Geometri .

Saya berharap menemukan ~ 140 perubahan dalam geometri, tetapi skrip saya bersikeras ada 430. Saya mencoba memeriksa representasi GeoJSON dan mereka identik, namun Kelas Geometri sama dengan () menolak untuk mengatakan demikian.

Contohnya di bawah ini:

>>> geom1geoJSON 
{'type': 'Polygon', 'coordinates': [[(-122.8423481559999, 47.060497293000083), (-122.84239755599992, 47.059262423000064), (-122.84416913599989, 47.059309693000046), (-122.84416913599989, 47.060497293000083), (-122.8423481559999, 47.060497293000083)]]}
>>> geom2geoJSON 
{'type': 'Polygon', 'coordinates': [[(-122.8423481559999, 47.060497293000083), (-122.84239755599992, 47.059262423000064), (-122.84416913599989, 47.059309693000046), (-122.84416913599989, 47.060497293000083), (-122.8423481559999, 47.060497293000083)]]}
>>> geom1 = arcpy.AsShape(geom1geoJSON)
>>> geom2 = arcpy.AsShape(geom2geoJSON)
>>> geom1.equals(geom2)
False
>>> geom2.equals(geom1)
False

Perilaku yang diharapkan di sini harus Benar (bukan Salah).

Adakah yang punya saran sebelum saya memindahkan semuanya ke geometri? (Saya ragu karena ogr.CreateGeometryFromGeoJSON () mengharapkan string, dan arcpy __geo_interface__mengembalikan kamus dan saya merasa seperti saya menambah kompleksitas tambahan).

Menemukan sumber daya berikut bermanfaat, meskipun mereka tidak menjawab pertanyaan:

  1. arcpy.Geometry pertanyaan di sini di gis.stackexchange.com yang ditautkan di atas dalam teks saya.
  2. Kesalahan di kelas Polygon arcpy dari forum arcgis.com (ternyata ada banyak kesalahan presisi di ArcGIS 10.0 yang secara teoritis diperbaiki pada 10.1 tapi saya tidak dapat memverifikasi bahwa, dalam 10.0 SP5 Anda masih mendapatkan kesalahan).
arcpy geometry Michalis Avraam
sumber

Jawaban:

12

Masalahnya kemungkinan besar adalah salah satu dari presisi floating point . Dalam kasus Anda, Anda telah mengekstraksi geometri menggunakan arcpy, dan Anda telah mencocokkannya dengan RUID Anda.

Untungnya karena Anda memiliki arcpy terpasang, Anda punya numpy, yang membuat membandingkan set array numerik menjadi mudah. Dalam hal ini saya akan menyarankan numpy. Allclose fungsi , yang tersedia di numpy 1.3.0 (diinstal dengan ArcGIS 10).

Dari sampel yang Anda berikan di atas

geom1geoJSON = {'type': 'Polygon', 'coordinates': [[(-122.8423481559999, 47.060497293000083), (-122.84239755599992, 47.059262423000064), (-122.84416913599989, 47.059309693000046), (-122.84416913599989, 47.060497293000083), (-122.8423481559999, 47.060497293000083)]]}
geom2geoJSON = {'type': 'Polygon', 'coordinates': [[(-122.8423481559999, 47.060497293000083), (-122.84239755599992, 47.059262423000064), (-122.84416913599989, 47.059309693000046), (-122.84416913599989, 47.060497293000083), (-122.8423481559999, 47.060497293000083)]]}

import numpy as np

close = np.allclose(np.array(geom1geoJSON["coordinates"]), np.array(geom2geoJSON["coordinates"]), atol=1e-7)
#Returns True

Itu atol kunci menentukan nilai toleransi.

Perhatikan bahwa Anda tidak boleh menggunakannya arcpy.AsShapesama sekali. Pernah. Seperti yang saya catat dalam pertanyaan ini (/ plug shameless) ada bug yang dikenal di ArcGIS yang memotong geometri ketika mereka dibuat tanpa sistem koordinat (bahkan setelah mengatur env.XYTolerancevariabel lingkungan). Dalam arcpy.AsShapetidak ada cara untuk menghindari hal ini. Untungnya geometry.__geo_interface__mengekstrak geometri yang benar dari geometri yang ada (meskipun tidak menangani poligon kompleks tanpa perbaikan dari @JasonScheirer).

om_henners
sumber
Terima kasih. Saya tidak berpikir menggunakan numpy untuk melakukan ini. Solusi lain tampaknya menggunakan modul desimal dan bekerja melalui itu, tetapi membutuhkan lebih banyak pekerjaan.
Michalis Avraam
Saya pikir itu akan menjadi penting untuk mengatur numpy.allclose() rtolparameter ke 0. Secara default itu 1e-05 dan itu dapat menyebabkan toleransi yang besar jika nilai array besar lihat: stackoverflow.com/a/57063678/1914034
Below Radar
11

Ketepatan koordinat akan menjadi pertimbangan penting di sini. Nomor titik apung tidak dapat disimpan dengan tepat.

Jika Anda menggunakan alat Bandingkan Fitur , apakah muncul dengan hasil yang diharapkan menggunakan toleransi XY default?

blah238
sumber
Saya tidak memeriksa alat Bandingkan Fitur karena alat yang saya bangun sebenarnya membandingkan fitur individual yang bergerak di antara kelas fitur yang berbeda. Dengan kata lain, suatu fitur dapat berpindah dari CityRoads ke CountyRoads, jadi saya perlu mencari tahu apakah ada yang berubah dalam geometri dan atribut selain dari kelas fitur yang menahannya. Ada total 24 kelas fitur, dan fitur dapat bergerak di antaranya. Feature Compare hanya akan membandingkan 2 kelas fitur, sehingga dapat memberi tahu saya jika tidak ada lagi di FC. Maka saya masih perlu membandingkan fitur untuk memastikan tidak berubah
Michalis Avraam
Saya memeriksa alat Bandingkan Fitur dengan toleransi default (8.983e-009 yang cukup kecil tapi ini adalah File GDB) dan melaporkan beberapa perubahan, tetapi bukan yang benar. Secara khusus, dikatakan ada 69 perubahan geometri (saya kira lebih baik dari sebelumnya) tetapi tampaknya mengasumsikan bahwa OID adalah cara untuk mengidentifikasi fitur unik (mencari OID1 lama dan OID1 baru) yang belum tentu benar (saya telah mengaturnya untuk menggunakan RUID saya sebagai semacam tetapi tidak suka). Jadi kembali ke papan gambar.
Michalis Avraam
4

di samping @ blah328 jawaban, Anda harus memilih untuk membandingkan dua tabel untuk melaporkan perbedaan dan persamaan dengan nilai tabular dan definisi bidang dengan Tabel Bandingkan .

Contoh:

import arcpy
from arcpy import env
arcpy.TableCompare_management(r'c:\Workspace\wells.dbf', r'c:\Workspace
\wells_new.dbf', 'WELL_ID', 'ALL', 'IGNORE_EXTENSION_PROPERTIES', 'WELL_DEPTH 0.001',
'#','CONTINUE_COMPARE', r'C:\Workspace\well_compare.txt'
Aragon
sumber
Terima kasih, saya akan memeriksanya ketika saya mencoba membandingkan data atribut. Untuk saat ini, sepertinya saya tidak dapat membandingkan geometri, yang lebih penting.
Michalis Avraam
3 
def truncateCoordinates(myGeometry)
    trucated_coords = []
    partnum = 0

    for part in (myGeometry):
        for pnt in myGeometry.getPart(partnum):
            if pnt:
                trucated_coords.append("{:10.4f}".format(pnt.X))
                trucated_coords.append("{:10.4f}".format(pnt.Y))
             else:
                continue
        partnum += 1     
    return truncated_coords

Jika .equals()fungsi tidak berfungsi seperti yang diharapkan dan / atau koordinat sedikit diubah di ArcGIS, Anda dapat memijat koordinat XY, lalu membandingkan setara String dari Geometri. Perhatikan, truncateCoordinates()matikan semua nilai di luar angka desimal ke-4.

geom1 = truncateCoordinates(feature1.Shape)
geom2 = truncateCoordinates(feature2.Shape)

geom1 == geom2
Klewis
sumber
@ klewis- Itu adalah salah satu cara untuk membandingkan geometri, tetapi rasanya seperti geometri. sama dengan (geometri) harus mengembalikan true ketika Anda membandingkan geometri yang sama persis. Memotong koordinat adalah semacam peretasan dalam arti tertentu. Mungkin ESRI perlu mulai menggunakan tipe desimal () alih-alih float jika mereka tidak dapat menangani nilai floating point dengan benar secara internal tetapi dapat mewakili mereka sebagai string yang sama.
Michalis Avraam