PostGIS: ST_Equals false ketika ST_Intersection = 100% dari geometri?

9

Saya memiliki 2 dataset yang terdiri dari data paket kadaster - masing-masing sekitar 125.000 baris. Kolom geometri adalah poligon WKB yang mewakili batas parsel; semua data secara geometris valid (poligon ditutup dll).

Beberapa data baru-baru ini tiba dalam proyeksi yang berbeda dengan data dasar yang digunakan untuk pekerjaan perbandingan - jadi saya memproyeksikan yang baru (basis adalah 4326; yang lain adalah WGA94 yang dibawa ke PostGIS sebagai 900.914 ... Saya memproyeksikan ulang ke 4326) .

Tahap pertama analisis adalah menemukan dan menyimpan paket yang tidak cocok; bagian dari itu adalah untuk mengidentifikasi dan menyimpan parsel dengan geometri yang identik.

Jadi saya menjalankan query yang sangat standar (blok kode di bawah ini meringkas detail skema dll):

create table matchdata as
  select  a.*
  from gg2014 a, gg2013 b
  where ST_Equals(a.g1,b.g1)

Hasil NOL.

"Aneh ..." pikirku. "Mungkin ada pergeseran vertex kecil yang disebabkan oleh proyeksi ulang: itu akan mengganggu, dan benar-benar tidak boleh terjadi."

Untungnya ada banyak data aspalal (5 kolom pengidentifikasi) yang memungkinkan saya membuat paket yang harus identik secara spasial: mereka yang memiliki pengidentifikasi yang sama, yang tanggal perubahannya di tabel 2014 adalah sebelum tanggal perubahan maksimal pada data 2013. Itu berjumlah 120.086 baris berbeda.

Saya menyimpan pengidentifikasi dan geometri dalam tabel terpisah ( match_id), dan menjalankan kueri berikut:

select apid, 
       bpid, 
       ST_Area(ag::geometry) as aa, 
       ST_Area(bg::geometry) as ab,
       ST_Area(ST_Intersection(ag,bg)::geometry)/ST_Area(ag::geometry) as inta,
       ST_Area(ST_Intersection(ag,bg)::geometry)/ST_Area(ag::geometry) as intb
from match_id
order by inta

16 nilai pertama untuk intadan intbidentik nol, 456 berikutnya adalah 0,99999999-ish (min 0,99999999999994, maks 0,99999999999999999), dan baris 473 dan seterusnya adalah 1 - hingga baris 120050, ketika area persimpangan lebih besar daripada geometri (terbesar) nilai untuk intadan intb1,00000000000029, tetapi masih).

Jadi inilah teka-teki saya: jika dua geometri berpotongan secara spasial antara 99,999999999994% dan 100,000000000029% dari daerah masing-masing, saya ingin "ST_Equals" untuk mengatakan "Yap .... Saya akan memberi Anda yang satu. Cukup dekat".

Bagaimanapun, itu setara dengan keluar sekitar 1 bagian dalam 16 triliun ... yaitu, seolah-olah utang nasional AS turun kurang dari 93 sen.

Dalam konteks keliling Bumi (pada ~ 40.000 km), rasanya seperti berada di ketinggian 0,0000000025 km, (karena menghasilkan perbedaan area yang kecil, setiap pergeseran titik harus lebih kecil).

Menurut TFD (yang saya punya R'd) toleransi untuk ST_Intersects()adalah 0,00001m (1mm), jadi perubahan tersirat dalam simpul (yang saya akui saya belum memeriksa: Saya akan ST_Dump()melakukannya dan melakukannya) akan tampak lebih kecil selain toleransi. (Saya menyadari itu ST_Intersects !== ST_Intersection(), tapi itu satu-satunya toleransi yang disebutkan).

Saya belum dapat menemukan toleransi yang sesuai untuk perbandingan titik yang dilakukan oleh ST_Equals()... tetapi tampaknya benar-benar aneh bahwa setidaknya 120.000 dari baris saya harus melewati penilaian yang masuk akal tentang identitas spasial, tetapi tidak.

(Catatan: Saya juga melakukan latihan yang sama menggunakan ::geography- dengan hasil yang memiliki lebih banyak variabilitas, tetapi masih lebih dari 110.000 entri dengan '1' bersih yang bagus).

Apakah ada cara untuk melonggarkan toleransi ST_Equals, yang tidak perlu menggali celah-celah kode? Saya tidak tertarik melakukan itu.

Jika tidak, adakah kludge yang diketahui orang?

Catatan: akan lebih baik jika 'kludge' tidak melakukan perbandingan bilateral seperti

where ST_within(g1, ST_Buffer(g2, 0.0000001))
  and ST_within(g2, ST_Buffer(g1, 0.0000001))


   - I've done that: sure, it works... but it's a gigantic documentation PITA).

Saya dapat mengatasi ini, tetapi menulis 20 halaman untuk mendokumentasikan penyelesaiannya - yang hanya akan muncul lagi jika kita mendapatkan data yang cerdik - adalah PITA yang saya lebih suka tidak harus melakukannya mengingat bahwa itu kemungkinan hanya sekali saja .

(Versi: Postgresql 9.3.5; PostGIS 2.1.3)

GT
sumber
Pikirkan di sini, tetapi apakah Anda sudah mencoba mengkanonikkan parsel baru ke kisi yang sesuai dengan data yang ada menggunakan st_snaptogrid?
nickves
Saya dapat mengerti tidak ingin melihat kode sumber, tetapi pertanyaan Anda mendorong saya untuk melakukannya (meskipun C ++ saya menyebalkan), jadi saya berterima kasih untuk itu. Jika Anda tertarik, saya dapat memposting bagian yang relevan, yang semuanya ada di github.com/libgeos .
John Powell
ST_Equalshanya mengembalikan trueketika geometri sama - tipe geometri, jumlah simpul, SRID, dan nilai simpul (dalam semua dimensi, dalam urutan yang sama). Jika ada varian, perbandingan berhenti, dan falsedikembalikan.
Vince
@Vince: seperti yang saya mengerti (dari dokumen), ST_Equals()mengabaikan directionality. Saya menganggap itu berarti bahwa untuk poligon 2-D tertutup, tidak ada bedanya jika titik-titiknya dihitung searah jarum jam vs berlawanan arah jarum jam. ST_OrderingEquals()adalah tes yang lebih ketat. Yang mengatakan, setelah memeriksa titik (menggunakan ST_Dump()dan menghitung delta untuk setiap titik) jelas bahwa jawaban luar biasa @John Barca adalah pada uang. ST_equals()dikontraindikasikan, bahkan untuk data identik-identik ex-ante , jika satu geometri diproyeksikan ulang - kecuali jika perbandingan dibuat dengan ST_SnapToGrid ().
GT.
Terlambat kembali ke ini: cara cepat yang bagus untuk mendapatkan tes yang dapat diterima untuk kesetaraan spasial dekat adalah untuk memeriksa berapa proporsi masing-masing geometri adalah bagian dari persimpangan. Ini sedikit memberatkan secara komputasi; menghitung (100*(ST_Area(ST_Intersection(a.g1, b.g1))/ST_Area(a.g1)))::int as int_pcadan (100*(ST_Area(ST_Intersection(a.g1, b.g1))/ST_Area(b.g1)))::int as int_pcb(pastikan Anda JOINmenyertakan ST_Intersects(a.g1,b.g1)). Tes jika (int_pca, int_pcb)=(100,100)(atau beberapa set cutoff lainnya). Kludgy, tapi itu akan melakukan 2,6 juta paket dalam ~ 30 menit (selama g1 diindeks GIST).
GT.

Jawaban:

20

Dugaan saya adalah bahwa Anda mengoordinasikan transformasi telah memperkenalkan kesalahan pembulatan kecil (lihat contoh di bawah). Karena tidak ada cara untuk mengatur toleransi dalam ST_Equals, ini menyebabkan ST_Equals mengembalikan false untuk beberapa geometri yang hanya berbeda di tempat desimal ke-n, karena geometri harus identik dalam segala hal - lihat definisi matriks persimpangan dalam libgeos . Anda dapat memeriksa ini dengan contoh yang sangat ekstrem,

SELECT ST_Equals(
      ST_MakePoint(0,0),
      ST_MakePoint(0,0.000000000000000000000000000000000000000000000000000000000001));

yang mengembalikan false .

Jika Anda menggunakan ST_SnapToGrid, Anda dapat memaksakan presisi yang diberikan, misalnya, ke sepuluh tempat desimal,

SELECT ST_Equals(
      ST_MakePoint(0,0),
      ST_SnapToGrid(
             ST_MakePoint(0,0.00000000000000000000000000000000000000000000001),
      10));

sekarang mengembalikan true .

Jika Anda lari,

CREATE table matchdata AS
SELECT  a.*
FROM gg2014 a, gg2013 b
WHERE ST_Equals(ST_SnapToGrid(a.g1, 5), ST_SnapToGrid(b.g1, 5));

pengaturan toleransi yang sesuai, saya kira masalah Anda akan hilang.

Berikut ini adalah tautan ke diskusi pengembang Postgis tentang toleransi yang menunjukkan bahwa penerapannya kurang dari sepele.

Saya melakukan beberapa konversi antara British National Grid (EPSG: 27700) dan lat / lon untuk mengilustrasikan poin tentang pembulatan presisi, mengambil titik di suatu tempat di London,

SELECT ST_AsText(ST_Transform(ST_SetSrid(ST_MakePoint(525000, 190000),27700),4326));

kembali POINT(-0.19680497282746 51.5949871603888)

dan membalikkan ini,

SELECT ST_AsText(ST_Transform(ST_SetSrid(ST_MakePoint(-0.19680497282746, 51.5949871603888),4326),27700));

kembali POINT(525000.000880007 189999.999516211)

yang mati kurang dari satu milimeter, tetapi lebih dari cukup untuk membuat ST_Equals kembali salah.

John Powell
sumber
Jawaban John Barca benar - bahwa kesalahan pembulatan kecil dapat membuang ST_Equals. Pengalaman saya (tidak menyenangkan) adalah ketika bekerja dengan dua set data - keduanya diproyeksikan dari EPSG 4326 ke EPSG 3857 - satu melalui ArcCatalog (ArcToolbox -> Alat Manajemen Data -> Proyeksi dan Transformasi) , sedangkan yang lain melalui GDAL ogr2ogr.
Ralph Tee
Jawaban ini sangat membantu saya. Tapi saya perhatikan bahwa indeks geografis tidak lagi digunakan dan permintaan terlalu lama. Solusi saya adalah membuat tabel sementara dengan geometri yang diambil dan menambahkan indeks sebelum menjalankan kueri. Apakah ada cara yang lebih baik untuk mempercepat?
hfs
1
@ hfs. Saya yakin Anda bisa membuat indeks fungsional menggunakan ST_SnapToGrid. Anda benar bahwa menggunakan panggilan fungsi di dalam yang sama / berpotongan / berisi operasi spasial dll akan menyebabkan indeks tidak digunakan dan membuat indeks fungsional akan menyelesaikan ini. Atau Anda dapat memperbarui data secara permanen jika menurut Anda ketepatannya palsu dan kemudian tidak harus menggunakan ST_SnapToGrid dalam kueri. Itu tergantung pada data Anda dan kasus penggunaan, tentu saja.
John Powell
2

Apakah Anda menjalankan pemeriksaan ST_IsValid pada geometri Anda? Jika tidak valid, semua taruhan dibatalkan. ST_Intersects dan keluarga lainnya dari fungsi hubungan spasial GEOS sering kali hanya akan kembali salah karena area tidak didefinisikan dengan baik dari sudut pandang matriks persimpangan. Alasan melakukan ST_Buffer mungkin berhasil adalah karena ia mengubah geometri tidak valid Anda menjadi yang valid. ST_Buffer (..., tinybit) adalah apa yang dikenal sebagai alat "orang miskin yang mencoba membuat geometri saya valid".

LR1234567
sumber
Langkah pertama dengan set data baru adalah memilih hanya menggunakan geometri yang valid ST_isValid(g1)- yang disebutkan (miring) "[] kolom geometri adalah WKB poligon yang mewakili batas-batas parsel; semua data secara geometris valid (poligon ditutup dll) ."
GT.
0

Jawaban saya datang agak terlambat, tetapi mungkin itu akan membantu seseorang yang memiliki masalah yang sama. Dari pengalaman saya, ketika dua geometri yang memang sama tetapi ST_Equals mengembalikan False dua hal dapat membantu:

  1. pastikan bahwa membandingkan geometri adalah geometri tunggal (Tidak Ada MultiLinesting, MultiPoin, dll.)
  2. coba ST_Equals(st_astext(a.geom), st_astext(b.geom)) bukanST_Equals(a.geom , b.geom)

Yang pertama sudah disebutkan dalam dokumentasi . Yang kedua tampaknya tidak rasional tetapi berhasil. Saya tidak tahu, tetapi kira itu ada hubungannya dengan format biner dari geometri postGIS default.

ioanna tsak
sumber