Umkreissuche/SRID Fragen

Für mich Ist PostGIS und allgemein GIS Neuland. Ich habe mir schon ein paar Postings in diesem Forum durchgelesen, aber nichts gefunden, was meine Fragen beantwortet.

Ich möchte in einer Tabelle Punkte (POIs) als Koordinaten, die ich als GPS Daten (Lat/Long dezimal, also SRID 4326) bekomme, ablegen.
Dazu habe ich mir eine Tabelle geschrieben in der ich die Punkte ablege:
Tabelle POI mit den Spalten:

  • poi_name: String
  • poi_lat: Float
  • poi_long: Float

Jetzt möchte ich (vorläufig) mit diesen Punkten folgendes tun. Gib mir alle Punkte die im Umkreis von X km liegen geordnet nach ihrer Entfernung und als 2. Use Case berechne die Entfernung zwischen zwei POIs.

Also hatte ich eine Spalte poi_geom4326 (mittels SELECT AddGeometryColumn…) angelegt, Index drauf gelegt und mit einem Trigger diese Spalte automatisch gefüllt. Leider waren die Ergebnisse die ich mit SQL erzielt habe nicht das was ich erwartet hatte: Distanzen von ca. 20-30 km ergaben Werte von 0.45, etc. Wenn ich es richtig verstanden habe liegt es daran, dass Geomwerte mit SRID 4326 nicht metrisch ‘behandelt’ werden können.

Also habe ich ein wenigweiter gesuch und bin auf ein SR System (Gauß-Krüger) mit der SRID 31467 gestossen, bei welchem man die Entfernungen in Metern ageben können soll. Hierfür habe ich eine weitere Geo-Spalte mit Index angelegt, die ebenfalls vom Trigger mit befüllt wird.

Die Trigger-Funktion sieht wie folgt aus:
CREATE OR REPLACE FUNCTION setGeoData() RETURNS trigger AS ’
BEGIN
NEW.poi_geom4326 := SetSRID(makepoint(NEW.lat, NEW.lon), 4326);
NEW.poi_geom31467 := Translate(NEW.poi_geom4326, 31467);
RETURN new;
END
’ LANGUAGE plpgsql;

Wenn ich nun die Distanz berechnen lasse bekomme ich zu grosse Werte heraus, z.B.
SELECT start.poi_name, ende.poi_name, distance(start.poi_geom31467,ende.poi_geom31467) FROM POI As start, POI As ende where start.poi_id=12 AND ende.poi_id=25
liefert eine Distanz von ca 36000 Metern (?). Die Luftlinie zwischen den beiden Punkten mit Google Earth handgemessen sind ca 27.1 km und
ST_Distance_Sphere liefert 27.3 km.

Jetzt könnte ich ja ST_Distanz_Sphere benutzen um alle Punkte innerhalb eines Radius zu bekommen, aber davon wurde abgeraten, weil es nicht performant genug ist, da es nicht über den Index gehen kann und für jeden Punkt individuell berechnet werden muss. Die Lösung welche ich gefunden haben benutzt die extend(mittelpunkt, radius) Method mit der intersects (&&) operation. Wenn ich dort den Radius in Metern angebe bekomme ich zu viele Ergebnisse.

Wie gesagt, es kann sein, dass ich ein grundlegendes Problem habe (Brett vor dem Kopf), da ich mich zum ersten mal mit GIS befasse. Von daher bin ich für jede Hilfe dankbar!

Hallo!

Also:
Im Prinzip hast du die meisten Dinge ganz richtig gemacht. Man kann die Funktionen natürlich auch verschachteln, was die 2. Geometriespalte überflüssig macht und generell ist es fixer und übersichtlicher, die Geometrie mit GeomFromText zu erzeugen. Du scheint ja aber generell mit der DB kein Problem zu haben, also überlasse ich das dem geneigten Leser zu Übung :wink:

Zur Ungenauigkeit:
Das Gauß-Krüger-System teilt den Globus in 3° breite Meridianstreifen und projiziert jeden einzelnen auf eine Zylinderwandung, in der die (Erd-)Kugel steckt. Dabei liegen die Pole auf der Wand, eine sog. transversale Mercatorprojektion. Die Punkte der Kugeloberfläche werden nun auf die Zylinderwand projiziert, dabei entstehen bei weit von der Schnittlinie zwischen Kugel und Zylinder große Verzerrungen, deshalb werden immer nur 3°-Streifen projiziert und dann die Kugel um 3° ‘weitergedreht’. Die jeweilige Nummers des Streifens(1,2,3,4,5…) wird als Kennung der x-Koordinate vorangestellt (z.b. 3587132). Trotzdem hattest du die richtige Idee (nur vermutlich den falschen Streifen) :wink:

Zur Auswahl hättest du in D:
31466 -> in WGS84 4,5-7,5° (nur ein Zipfelchen im Westen)
31467 -> in WGS84 7,5-10,5°
31468 ->in WGS84 10,5-13,5°
31469 -> in WGS84 13,5-16,5° (nur ein Zipfelchen im Osten)

Es hängt also auch davon ab, wo deine WGS84-Koordinaten liegen, welchen sog. GaußKrügerstreifen du nehmen musst.

Ansonsten gibts zur Entfernung nichts mehr zu sagen:

Die Umkreisabfrage kannst du ähnlich gestalten, in dem du die Distanz im 3146x-System als Where-Kriterium nimmst, dann kannst du einerseits den index benutzen, andererseits in m rechnen, zumindest solange die Entfernungen nicht durch ganz Deutschland gehen. Andernfalls könntest du einen Buffer um den Punkt konstruieren und schauen, welche Punkte dadrin liegen.
Soweit mal

Gruß Uwe

Hi,

vielen Dank für die ausführliche Antwort. Wenn ich es also richtig verstanden habe liegen die Punkte, die ich in meiner Tabellespalte speichere im falschen 3° Streifen. Falls dem so ist, dann ist meine verwirrung allerings noch grösser geworden, denn die Punkte liegen prinzipiell über ganz Deutschland (sogar ganz Europa) verteilt. Da ich beim Anlegen der Geo-Spalte schon den SRID (und damit den 3° Streifen) angeben muss, kann ich die Punkte also gar nicht in einer Spalte ablegen, sondern muss für jeden möglichen 3° Streifen eine eigene Spalte anlegen, oder habe ich da einen Denkfehler? Wie finde ich dann den nächsten Punkt zu einem gegebenen Punkt; vor allem wenn er am Rand eines 3° Streifens liegt?
Gibt es vielleicht eine “universellere” SRID, die besser geeignet ist?

Naja, die universellste ist halt WGS84, also das was du im Orginal hast. Verzerrungsfrei arbeiten über grosse Entfernungen lässt sich halt nur mit Winkeln. Wer den Blödsinn mit der Erd-‘KUGEL’ erfunden hat, gehört geschlagen… :wink:
An deiner Stelle würde ich nur eine Geometriespalte für die Orginaldaten in WGS84 verwenden. Das vermindert den Overhead, den du mit mehreren Geometrien hast. Daraus könntest du z.B. auch einen View erstellen, in dem die jeweiligen ‘Rechengeometrien’ mit einem case-Statement, das die unterschiedlichen Streifen unterscheidet, erzeugt werden.
Wenn du allerdings über Streifengrenzen hinweg rechnen willst, wäre WGS84 das am ehesten geeignete, andernfalls hast du jede Menge Rechnerei.
Ein anderes, etwas weniger aufwändiges System wäre UTM, da sind die Streifen 6° breit, ansonsten funktioniert es analog zu Gauß-Krüger.

Gruß Uwe

´Mir persönlich wäre es mittlerweile auch lieber die Erde wäre eine Scheibe.

Ich möchte aber trotzdem noch einmal Nachharken, weil ich das Ergebnis etwas ernüchternd finde. Ich habe in meiner Tabelle ca 400.000 POIs und möchte zu einem gegebenen Punkt den nähesten in der Tabelle finden, also habe ich zunächst folgendes gemacht:

SELECT poi_name, ST_Distance_Sphere(poi_geom4326, GeomFromText(‘POINT(50.1 8.7)’, 4326)) As dist FROM POI ORDER BY dist

Die Abfrage dauert: 11.3 sec. Das ist zu langsam, also habe ich mit LIMIT gespielt:

SELECT poi_name, ST_Distance_Sphere(poi_geom4326, GeomFromText(‘POINT(50.4 8.7)’, 4326)) As dist FROM POI ORDER BY dist LIMIT 1

Die Abfrage dauert: 600 ms. Das ist schon viel besser.

Dann bin ich aber auf das hier gestossen:
Solving the Nearest Neighbor Problem in PostGIS
Solving the Nearest Neighbor Problem in PostGIS

Dort wird folgender Weg vorgeschlagen. Index auf die the_geom Spalte und dann:
SELECT g1.gid As gref_gid, g1.description As gref_description, g2.gid As gnn_gid,
g2.description As gnn_description
FROM sometable As g1, sometable As g2
WHERE g1.gid = 1 and g1.gid <> g2.gid AND expand(g1.the_geom, 300) && g2.the_geom
ORDER BY distance(g1.the_geom,g2.the_geom)

Allerdings funktioniert das bei meinen Daten nicht, weil the_geom bei mir SRID 4326 sind
und die 300 Units die man im expand übergibt auf jeden Fall keine Meter sind (ja was sind
das eigentlich?).

Meine Frage ist kann ich mein Problem überhaupt performanter Lösen als mit
SELECT poi_name, ST_Distance_Sphere(poi_geom4326, GeomFromText(‘POINT(50.4 8.7)’, 4326)) As dist
FROM POI ORDER BY dist LIMIT 1 ?

Sorry, ich hatte schon mal geantwortet aber anscheinend das Abschicken vergessen :angry:

Sphärische Berechnungen sind halt recht aufwändig, da kommt man nicht drumrum. Der Weg im Link (Danke!) basiert darauf die Daten einzuschränken, da ist dann der Suchraum ebenfalls beschränkt. Wenn man da eine sinnvolle Annahme machenkann ist da eine Möglichkeit das zu beschleunigen. Die Einheiten, die in WGS 84 verwendet werden sind übrigens Grad. D.h. Du gibst den Radius halt nicht in m an sondern als Öffnungswinkel eines Kegel mit Spitze im Erdmittelpunkt. So könntest du übrigens auch rechnen. Da du WGS84 Koordinaten hast, könntest du dir die Winkeldifferenz (die ist proportional zur Distanz auf der Kugeloberfläche) direkt ausrechnen, anstatt über Geometrien zu gehen, das dürfte ebenfalls schneller sein. Du sparst dir auch die ganzen Hinundher-Projektionen.

Gruß Uwe

PS: Wenn du Limit 1 nimmst, dann nimmt er den ersten Besten, logisch, das er dann nicht suchen muß… :wink:

Sorry, dass ich den alten Thread noch mal aufwärme, aber ich bin neu dabei :astonished:

Gismans Vorschlag mit der Errechnung der Entfernung aus der Winkeldifferenz fand ich sehr interessant und musste ihn gleich mal ausprobieren. Ich hab nämlich ein sehr ähnliches Problem wie Vermonster (Umkreissuche) und möchte meine Suche möglichst effizient einschränken, ohne dabei Ergebnisse zu verlieren. Anders gesagt, ich brauche eine Entfernungsberechnung, die möglichst genaue, aber auf keinen Fall zu große Werte liefert. Zu klein ist nicht so schlimm.

Nun sehe ich folgendes:

select st_distance(
st_pointfromtext(‘POINT(8.680833 50.111111)’, 4326),
st_pointfromtext(‘POINT(8.680506 50.111512)’, 4326));

…liefert mir 0.000517426.

select st_distance(
st_transform(st_pointfromtext(‘POINT(8.680833 50.111111)’, 4326), 31467),
st_transform(st_pointfromtext(‘POINT(8.680506 50.111512)’, 4326), 31467));

…liefert mir 50.365.

Wenn ich bisher alles richtig verstanden habe, heißt dass, das ich mit Gauß-Krüger Zone 3 (SRID 31467) einen Abstand von 50,365 Metern habe, und ohne Projektion einen Abstand von 0,000517426 Grad. Wenn ich nun die Grad in Meter umrechnen will, hätte ich erwartet, dass ich sie mit (Erdumfang in m / 360) multiplizieren müsste. Der Erdumfang ist an der kleinsten Stelle knapp 40.000 km, damit komme ich auf einen Multiplikator von ca. 110000 und somit auf eine Entfernung von ca. 57 Metern.
Nach dieser Rechnung dürfte der Erdumfang höchstens 35.000 km betragen, damit die Umrechnung funktioniert.

Bei größeren Abständen sind die Ergebnisse übrigens noch schlimmer (Beispiel Distanz Hamburg - Berlin: 256 km in Gauß-Krüger 3, ohne Projektion 3,5676 Grad -> Erdumfang höchstens 25.000 km!?!).

Habe ich einen Denkfehler? Wenn nein, woher kommt diese deutliche Abweichung, und was für Abweichungen muss ich maximal befürchten? Gerade die letzte Teilfrage ist wahrscheinlich nicht unbedingt trivial :smiley: . Vielleicht kann sie ja trotzdem jemand beantworten - ich würde mich freuen…