Schnittalgorithmus, der 180 ° Meridian und Pole korrekt handhabt?

9

Gibt es einen bekannten Schnittalgorithmus, der 180 Meridiane und Pole korrekt verarbeitet?

Nehmen wir zum Beispiel an, wir haben eine Liste von Breiten- und Längengraden, die die Antarktis darstellen. Nehmen wir an, wir haben auch ein einfaches Polygon, das ein Flugzeug darstellt. Wir wollen wissen, ob sich das Flugzeug über der Antarktis befindet.

Der generische 2D-Polygonschnitt funktioniert in diesem Fall nicht, wenn Sie einfach den Breitengrad für y und den Längengrad für x verwenden, da das flache Koordinatensystem Kanten mit einem Längengrad von -180 und 180 sowie einem Breitengrad von -90 und 90 aufweist. Die Antarktis wird an drei Rändern von der Seite verschwinden.

Sean
quelle

Antworten:

14

GIS (als Feld) hat es nicht allzu heiß gemacht, wenn es darum geht, sich wirklich mit der Oberfläche des Globus auseinanderzusetzen.

Zum Beispiel ist Ihr Problem nicht vollständig definiert. Anders als in 2D, wo wir wissen, dass die Kanten eines Polygons aus geraden Linien bestehen, was sind sie auf dem Globus? Bögen mit großen Kreisen, die den Abstand zwischen Scheitelpunkten minimieren, sind eine gute Wahl, aber nicht die einzige. Eine andere Wahl sind beispielsweise gerade Linien (die sich somit unter der Erdoberfläche bewegen). (Rhumb-Linien sind ebenfalls eine Wahl, aber möglicherweise eine dumme, insbesondere in der Nähe der Pole ...) Hier ist ein großartiger Link zu Problemen bei der Interpretation von WRT-Kanten: http://blog.opengeo.org/2010/08/10/shape-of -a-Polygon /

Ein weiteres großes Problem beim Projekt-zu-2d-dann-schneiden-Ansatz ist neben den anderen erwähnten Singularitäten, dass alle neu erstellten Schnittpunkte (an denen sich Polygonkanten kreuzen) fehl am Platz sind und von der verwendeten Projektion abhängen. (Ich glaube, hier kommt die Empfehlung für die Verdichtung her: Durch Hinzufügen einer Tonne Zwischenscheitelpunkte erhalten Sie einen reduzierten Fehler durch die Projektion der Mitte der Polygonkanten.)

Angenommen, Sie möchten nicht alle Kompromisse und Problemumgehungen bei der Projektion auf 2D eingehen und möchten selbst etwas durchdenken und codieren, dann habe ich ein wenig Prototyp-Code (und nur Prototyp-Code!) Für einen Kunden erstellt mit diesem.

Hier ist eine Skizze des Ansatzes. Sie müssen wissen, was ein Vektor ist und welche Bedeutung die Punkt- und Kreuzprodukte haben . (Vorsichtsmaßnahme: Wikipedia ist praktisch für schnelle Links, aber ausreichend, wenn es Ihre erste Einführung in die Themen ist. Ein gutes Tutorial für 3D-Grafiken hilft.)

  • Stellen Sie einen Punkt auf der Kugel durch einen kartesischen 3D-Vektor dar. Der Punkt auf der Erdoberfläche ist der Punkt, an dem der Vektor, wenn Sie ihn zu einem unendlich langen Strahl ausdehnen, die Erdoberfläche schneiden würde.
  • Stellen Sie große Kreise durch eine Ebene durch den Ursprung dar. (In 3D reicht ein einzelner Einheitsvektor aus, um eine Ebene durch den Ursprung zu definieren. Er ist die Normale zur Ebene.) Der Großkreis ist der Schnittpunkt der gesamten Ebene mit der Erdoberfläche.
  • Sie können die Schnittpunkte zweier Großkreise finden, indem Sie ihre beiden Ebenen schneiden.
  • Definieren Sie einen Bogen eines Großkreises durch zwei Punkte. Der Vektor, der die Großkreisnormalen definiert, ist das Kreuzprodukt der Vektoren der Start- und Endpunkte.
  • Um festzustellen, ob ein Punkt, von dem wir wissen, dass er auf dem Großkreis liegt, innerhalb des Bogens liegt, erstellen Sie zwei Ebenen, sodass: Sie senkrecht zur Ebene des Großkreises stehen, eine den Startpunkt und die andere den Endpunkt enthält und sie sind einander zugewandt. Dann liegt der Punkt auf dem Bogen, wenn er im Inneren beider Ebenen liegt. (Zur Veranschaulichung: Sie haben das Lächeln eines Pac-Mannes erzeugt, als er auf den Bogen kaut. Wenn sich der Testpunkt zwischen seinen Kiefern befindet, liegt er auf dem Bogen, wie wir bereits wissen, dass er sich auf dem großen Kreis befindet.)
  • Um festzustellen, ob sich zwei Bögen kreuzen: Suchen Sie die beiden Schnittpunkte der entsprechenden Großkreise und testen Sie jeden Punkt, um festzustellen, ob er in beiden Bögen liegt.
  • Eine weniger mehrdeutige Definition: Ein Polygon ist eine Sammlung von Punkten, die jeweils durch Kanten verbunden sind, die aus Bögen großer Kreise bestehen. Die Punkte sind so angeordnet, dass, wenn Sie entlang der Erdoberfläche entlang der Kanten des Polygons gehen, das „Innere“ des Polygons zu Ihrer Linken liegt. Lassen wir zunächst komplexe Polygone (Inseln, Löcher und Selbstüberschneidungen) heraus.
  • Anhand des Vorzeichens des Punktprodukts der entsprechenden Vektoren können Sie erkennen, ob sich ein Punkt auf der rechten oder linken Seite einer Ebene befindet. (Dies entspricht, wenn Sie sich um den Großkreis bewegen, ob der Punkt auf der Oberfläche der Kugel zu Ihrer Linken oder zu Ihrer Rechten liegt.)
  • Ein genauer Test, um festzustellen, ob sich ein Punkt innerhalb eines Polygons befindet: Liegt er auf der linken Seite aller Kanten?
  • Jetzt können wir testen, ob sich ein Punkt innerhalb eines Polygons befindet, und Kantenübergänge bestimmen: die Zutaten für Polygon-Polygon-Schnittpunkte! Der Rand dieses Kommentars ist zu klein, um einen vollständigen Algorithmus zu schreiben. Die grundlegenden Schritte bestehen jedoch darin, (a) alle Schnittpunkte zu finden und dann (b) Kanten zu gehen und abwechselnd das Polygon zu wechseln, auf dem Sie laufen, wenn Sie auf Schnittpunkte stoßen.
  • Wenn alles oben Genannte funktioniert, sollten Sie über Indizierungsstrategien nachdenken, um es schneller zu machen, da die Umrisse des Punkt-in-Polygons I O (n) in der Anzahl der Kanten und der Polygonschnitt O (m * n) in der Anzahl von Kanten sind Kanten.

Pheew.

Dieser Ansatz bietet einige große Vorteile: Alle oben genannten Operationen beschränken sich nur auf Multiplikationen und Additionen. (Nach dem Konvertieren von Daten in diese Darstellung: z. B. Lat / Long-Koordinate erfordern einige Trigger, um einen kartesischen XYZ-Vektor zu erhalten.) Es gibt keine Singularitäten an den Polen oder im Koordinatensystem und nicht zu viele Sonderfälle, um die man sich Sorgen machen muss (Beispiel Sonderfall) : parallele überlappende Kanten).

Wenn man sich den Code ansieht, sieht es so aus, als ob das Spheres- Paket, das jemand anderes verlinkt hat, einem Teil dieses Ansatzes folgt, obwohl es auch ein bisschen halbherzig aussieht.

PostGIS verwendet möglicherweise auch einen ähnlichen Ansatz für den geografischen Datentyp , aber ich habe nicht unter die Haube geschaut. Ich weiß, dass sie zumindest für die räumliche Indizierung einen R-Baum über 3D-Kartesisch verwenden.

(Hinweis: Diese Antwort wurde lang genug, dass ich sie wahrscheinlich in einen Blog-Beitrag umwandeln werde ... Feedback / Kommentare sehr willkommen!)

Dan S.
quelle
5

Ich suche tatsächlich nach einer allgemeinen Lösung, die für viele große Polygone funktioniert.

Die meisten Antworten, die ich bisher gelesen habe, konzentrieren sich natürlich darauf, eine geeignete Projektion für Ihre Funktionen zu finden. Dies mag kontraintuitiv sein, aber hier müssen Sie überlegen, wo eine Projektion nicht funktioniert und nicht, wo sie gut funktioniert.

Einige Projektionen funktionieren nur an einem Punkt nicht: Diese gehören zu den sogenannten "globalen" Projektionen. Ein gutes Beispiel ist der azimutale Äquidistant: Der einzige Punkt, den Sie nicht projizieren können, ist der Punkt, der seinem Ursprungsort diametral entgegengesetzt ist. (Eigentlich ist es nicht unbedingt richtig, dass Sie diesen letzten Punkt nicht projizieren können. Ich beschönige ein technisches Problem. Es genügt, dass bei der Projektion am gegenüberliegenden Pol schwerwiegende Probleme auftreten.)

Wenn Sie also einen einzelnen Punkt P außerhalb der Vereinigung aller Polygone finden, mit denen Sie arbeiten möchten, müssen Sie nur innerhalb einer globalen Projektion arbeiten, die nur an diesem einzelnen Punkt fehlschlägt. ( Verwenden Sie z. B. ein schräges äquidistantes Azimuthal, das an dem Punkt diametral gegenüber P zentriert ist .)

(Sie können einen solchen Punkt automatisch finden. Erstellen Sie beispielsweise ein Polygon, das die Erde bedeckt, puffern Sie die Vereinigung aller Ihrer Features leicht, entfernen Sie diesen Puffer aus dem Erdpolygon und wählen Sie einen beliebigen Punkt innerhalb der verbleibenden Punkte aus. Der Grund für die Pufferung ist, dass es am besten ist, sich von den einzelnen Punkten fernzuhalten, an denen eine Projektion fehlschlägt, je weiter desto besser: Dies ist ein Problem der Rechengenauigkeit.)

Es gibt viele globale Projektionen mit wünschenswerten Eigenschaften, einschließlich äquidistanter, konformer ( z. B. stereografischer) und flächengleicher. Daher ist die Auswahl einer globalen Projektion für Ihre Analyse keine schwerwiegende Einschränkung.

Die Hauptbeschränkung besteht darin, dass die beliebteste Software, ArcGIS, nicht viele schräge Versionen globaler Projektionen unterstützt! :-(

whuber
quelle
3

Sie können dies mit Projektionen lösen!

Nur weil die "Enden" der Skalen bei -180 und +180 liegen, heißt das nicht, dass Sie die Karte mit diesen Enden berücksichtigen müssen.

Verwenden Sie EPSG: 32761 ( http://spatialreference.org/ref/epsg/32761/ )

Konvertieren Sie Ihre Breiten- und Längengrade in Koordinaten im USV-Raum. Anschließend können Sie mithilfe regelmäßiger Geometrieberechnungen feststellen, ob sich Ihr Flugzeug über der Antarktis befindet.

Weitere Informationen zum Universal Polar Stereographic-Koordinatensystem finden Sie hier: http://en.wikipedia.org/wiki/Universal_Polar_Stereographic_coordinate_system

Wenn Sie nur versuchen, einen Schnittalgorithmus für gigantische Polygone zu bestimmen, behandeln Sie die Erde als Sphäroid (oder sogar als Kugel) und analysieren Sie Ihre Polygone mithilfe der sphärischen Geometrie .

Ich habe eine Bibliothek namens Spheres gefunden , die GPLed ist, sodass Sie ihren Code überprüfen können, um genau zu sehen, welche Algorithmen sie für die Bogenschnittstelle und "Punkte innerhalb eines sphärischen Polygons" verwenden.

Wenn Sie diese Polygone jedoch auf einer ebenen Fläche (wie Papier oder Ihrem Bildschirm) zeichnen möchten, müssen Sie eine Projektion oder zumindest einen Projektionsstil auswählen .

Walker
quelle
Aber wie würde ich bestimmen, welche Projektion verwendet werden soll? Ich bin damit einverstanden, dass Sie im speziellen Fall der Antarktis UPS verwenden können, aber was ist, wenn Sie nur ein beliebiges Polygon haben, das die Hemisphären kreuzt? Vielleicht ist ein besseres Beispiel als die Antarktis ein Polygon, das darstellt, wo derzeit Tageslicht ist.
Sean
"Die richtige Projektion auswählen" wäre wahrscheinlich ein großartiges Community-Wiki-Thema. Wenn Sie eine Projektion ausgewählt haben, liegt dies im Allgemeinen daran, dass Sie bereits zu einem bestimmten Fall gelangt sind. Wenn Sie nicht wissen, wo auf der Welt sich Ihr Polygon befindet, empfehle ich, den Ursprung der Erde in den Mittelpunkt Ihres Polygons zu transformieren und dann die kartesischen Koordinaten zu verwenden.
Walkerer
Ich suche nach einer allgemeinen Lösung, die mit vielen großen Polygonen funktioniert. Wenn ich für jedes Polygon eine benutzerdefinierte Projektion verwenden muss, scheint die Leistung langsam zu sein.
Sean
1

Während es Probleme bei der Verwendung kartesischer Koordinaten gibt, sehe ich die Reibung in dem von Ihnen bereitgestellten Beispiel nicht. So sieht es für mich aus, mit einem Polygon, das die Antarktis bedeckt und von den Rändern des Koordinatenraums begrenzt wird:

Antarktis
(Quelle: geohack.net )

Der Schnittpunkt zwischen Ebene und Kontinent kann ohne komplexe Transformationen sicher berechnet werden.

scw
quelle
Dies funktioniert in der Tat und Sie müssen nicht in eine andere Projektion konvertieren. Sie müssen nur das Polygon "Antarktis" entlang der drei Kanten schließen, die es kreuzt. Ich würde auch den Schwerpunkt des Flugzeugpolygons verwenden, da es einfacher ist zu bestimmen, ob sich ein Punkt innerhalb eines Polygons befindet, als eine Polygonkreuzung durchzuführen.
Walkerer
Ich mag diesen Ansatz nicht besonders. Es treten starke Verzerrungen auf, und der Schnittpunkt kann sich als wahr herausstellen, obwohl dies nicht der Fall ist.
George Silva
Ich kann mir einige Fälle vorstellen, in denen es wichtig wäre , dies zu projizieren , z. B. die Berechnung von Großkreisentfernungen oder die Berechnung der Fläche in der Nähe der Pole. Vielleicht fehlt mir etwas, aber ich sehe nicht, wie sich die hier angeforderte Kreuzung unter diesen Bedingungen ändern würde.
Scw
Ich denke, dass das Flugzeugbeispiel nicht perfekt war; Ein besseres Beispiel könnte eine große Wolke sein.
Sean
Eigentlich, scw, bist du richtig. Für die Punktanalyse wäre dies überhaupt nicht wichtig, nur für Entfernungs- und Flächenberechnungen. Mein Fehler: P
George Silva
0

Eine Lösung für dieses Problem besteht darin, anstelle eines zweidimensionalen Koordinatensystems ein dreidimensionales Koordinatensystem für die Erde zu verwenden. Wenn den 2D-Polygonen eine große virtuelle Tiefe zugewiesen wird, z. B. auf halbem Weg zum Erdmittelpunkt, müssen die Polygone nicht als Kurven auf der Erdoberfläche modelliert werden. Es sollte ausreichen, dass sich nur die Verticies an der Oberfläche befinden.

Sean
quelle
0

Eine mögliche Lösung besteht darin, das Polygon an allen Kanten (180 Meridian und Pole) in Unterpolygone zu schneiden. Die Software könnte sogar eine Art Referenz zwischen dem ursprünglichen Polygon und den Unterpolygonen beibehalten.

Sean
quelle
0

Nehmen Sie den Schwerpunkt des Flugzeugs und verwenden Sie ihn als Tangentialpunkt für eine azimutale Projektion . Sie müssen die Geometrien auch verdichten, bevor Sie sie für Fälle wie den Bereich der Erde projizieren, in dem derzeit Tageslicht zu sehen ist.

Kirk Kuykendall
quelle
Das klingt teuer. Was ist, wenn ich viele Flugzeuge oder schlimmer noch Wolken habe? Bitte erläutern Sie auch, was unter "verdichten" zu verstehen ist.
Sean
Ja, das wäre teuer. Um dies global zu tun, ist Planar vielleicht doch nicht der richtige Weg. Es scheint, als könnte ein sphärischer Trigonometrieansatz als Grundlage für eine Vektorpolygonüberlagerung von weltweiten Datensätzen verwendet werden (z. B. wie viel Prozent der Erde von Tageslicht und Wolken bedeckt sind). Alle mir bekannten Polygon-Overlays sind jedoch planar. en.wikipedia.org/wiki/Spherical_trigonometry
Kirk Kuykendall