PostGIS: ST_Equals false, wenn ST_Intersection = 100% der Geometrie?

Ich habe 2 Datensätze, die aus Katasterpaketdaten bestehen - jeweils ungefähr 125.000 Zeilen. Die Geometriespalte besteht aus WKB-Polygonen, die Paketgrenzen darstellen. Alle Daten sind geometrisch gültig (die Polygone sind geschlossen usw.).

Einige neuere Daten kamen in einer anderen Projektion an als die Basisdaten, die für einen Vergleichsjob verwendet wurden. Daher habe ich die neuere neu projiziert (Basis war 4326; die andere war WGA94, die als 900914 in PostGIS aufgenommen wurde ... ich habe sie auf 4326 projiziert). .

Die erste Stufe der Analyse bestand darin, nicht übereinstimmende Pakete zu finden und zu speichern. Dazu gehört die Identifizierung und Lagerung von Paketen mit identischer Geometrie.

Also habe ich eine sehr standardmäßige Abfrage ausgeführt (der folgende Codeblock abstrahiert Schemadetails usw.):

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

NULL Ergebnisse.

"Seltsam ..." dachte ich. "Vielleicht gab es winzige Scheitelpunktverschiebungen, die durch die Reprojektion verursacht wurden: Das wäre ärgerlich und sollte wirklich nicht passieren."

Glücklicherweise gibt es reichlich aspatiale Daten (5 Identifizierungsspalten), mit denen ich Pakete erstellen kann, die räumlich identisch sein sollten: solche mit derselben ID, deren Änderungsdatum in der Tabelle 2014 vor dem maximalen Änderungsdatum in den Daten 2013 lag. Das waren 120.086 verschiedene Zeilen.

Ich habe die Bezeichner und Geometrien in einer separaten Tabelle ( match_id) gespeichert und die folgende Abfrage ausgeführt:

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

Die ersten 16 Werte für intaund intbwaren identisch Null, die nächsten 456 waren 0,99999999-ish (min 0,99999999999994, max 0,999999999999999) und die Zeilen ab 473 waren 1 - bis Zeile 120050, als die Fläche des Schnittpunkts größer als jede der Geometrien war (die größte Wert für intaund intbwar 1.00000000000029, aber immer noch).

Hier ist mein Rätsel: Wenn sich zwei Geometrien räumlich zwischen 99,999999999994% und 100,000000000029% ihrer jeweiligen Bereiche schneiden, möchte ich, dass "ST_Equals" "Ja ... ich gebe Ihnen diese. Nah genug" sagt.

Immerhin entspricht dies einem Ausfall von etwa 1 Teil von 16 Billionen ... dh als ob die US-Staatsverschuldung um weniger als 93 Cent gesunken wäre.

Im Zusammenhang mit dem Erdumfang (bei ~ 40.000 km) ist es so, als wäre man um 0,0000000025 km oben (da dies zu einem so kleinen Flächendifferenz führen muss, muss jede Scheitelpunktverschiebung noch kleiner sein).

Laut TFD (für das ich R'd habe) beträgt die Toleranz ST_Intersects()fiktiv 0,00001 m (1 mm), sodass die impliziten Änderungen in den Eckpunkten (von denen ich gestehe, dass ich sie nicht überprüft habe: Ich werde ST_Dump()sie tun und dies tun) kleiner zu sein scheinen als die Toleranz. (Mir ist klar ST_Intersects !== ST_Intersection(), aber es ist die einzige erwähnte Toleranz).

Ich konnte die entsprechende Toleranz für den von ST_Equals()... durchgeführten Scheitelpunktvergleich nicht herausfinden, aber es scheint wirklich seltsam, dass mindestens 120.000 meiner Zeilen eine vernünftige Bewertung der räumlichen Identität bestehen sollten, aber nicht.

(Hinweis: Ich habe dieselbe Übung auch mit durchgeführt ::geography- mit Ergebnissen, die variabler waren, aber immer noch mehr als 110.000 Einträge mit einer schönen sauberen '1').

Gibt es eine Möglichkeit, die Toleranz von ST_Equals zu lockern, ohne in die Zwischenräume des Codes zu graben? Darauf bin ich nicht scharf.

Wenn nicht, gibt es einen Kludge, den jemand kennt?

Hinweis: Es wäre gut, wenn der 'Kludge' keinen bilateralen Vergleich wie machen würde

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).

Ich kann das umgehen, aber das Schreiben der 20 Seiten zur Dokumentation der Problemumgehung - die nur dann wieder auftaucht, wenn wir zweifelhafte Daten erhalten - ist eine PITA, die ich lieber nicht tun müsste, da es sich wahrscheinlich um eine einmalige handelt .

(Versionen: Postgresql 9.3.5; PostGIS 2.1.3)

Ich vermute, dass Sie Koordinatentransformationen zu winzigen Rundungsfehlern geführt haben (siehe Beispiel unten). Da es keine Möglichkeit gibt, die Toleranz in ST_Equals festzulegen, führt dies dazu, dass ST_Equals für einige Geometrien, die sich nur in der n-ten Dezimalstelle unterscheiden, false zurückgibt, da die Geometrien in jeder Hinsicht identisch sein müssen - siehe die Definition der Schnittmatrix in libgeos . Sie können dies anhand eines wirklich extremen Beispiels überprüfen:

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

was false zurückgibt .

Wenn Sie ST_SnapToGrid verwenden , können Sie eine bestimmte Genauigkeit festlegen , z. B. zehn Dezimalstellen.

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

gibt jetzt true zurück .

Wenn du rennen würdest,

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

Wenn Sie eine angemessene Toleranz festlegen, werden Ihre Probleme vermutlich verschwinden.

Hier ist ein Link zu einer Postgis-Entwicklerdiskussion über Toleranz, die darauf hinweist, dass die Implementierung weniger trivial ist.

Ich habe ein paar Konvertierungen zwischen British National Grid (EPSG: 27700) und lat / lon durchgeführt, um den Punkt der Rundungsgenauigkeit zu veranschaulichen, indem ich irgendwo in London einen Punkt genommen habe.

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

kehrt zurück POINT(-0.19680497282746 51.5949871603888)

und dies umzukehren,

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

kehrt zurück POINT(525000.000880007 189999.999516211)

Das ist weniger als ein Millimeter entfernt, aber mehr als genug, um ST_Equals false zurückzugeben.

Haben Sie Ihre Geometrien mit ST_IsValid überprüft? Wenn sie ungültig sind, sind alle Wetten ungültig. ST_Intersects und die andere Familie von GEOS-Funktionen für räumliche Beziehungen geben häufig nur false zurück, da der Bereich aus Sicht der Schnittmatrix nicht genau definiert ist. Der Grund, warum ST_Buffer wahrscheinlich funktioniert, liegt darin, dass Ihre ungültigen Geometrien in gültige konvertiert werden. ST_Buffer (..., tinybit) ist das sogenannte Werkzeug "Der Versuch eines armen Mannes, meine Geometrie gültig zu machen".

Meine Antwort kommt etwas spät, aber vielleicht hilft sie jemandem, der das gleiche Problem hat. Nach meiner Erfahrung können zwei Dinge helfen, wenn zwei Geometrien gleich sind, aber ST_Equals False zurückgibt:

1. Stellen Sie sicher, dass beim Vergleichen von Geometrien einzelne Geometrien vorliegen (kein MultiLinesting, MultiPoin usw.).
2. versuche ST_Equals(st_astext(a.geom), st_astext(b.geom)) stattST_Equals(a.geom , b.geom)

Die erste ist bereits in der genannten Dokumentation . Der zweite scheint irrational, funktioniert aber. Ich weiß es nicht, aber ich denke, es hat mit dem Binärformat der Standard-PostGIS-Geometrie zu tun.