Kontakt
terrestris

Aktualisierung unserer maps.mundialis-Anwendung löst Problem der Datumsgrenzen-Überschreitung

Das Update unserer Webapplikation maps.mundialis.de beinhaltet neue Features, die nachfolgend beschrieben werden. Um direkt zum Abschnitt zu springen, der die Datumsgrenzen-Überschreitung behandelt, hier klicken.

 

Neue Funktionen

 

Ortssuche


Mit der neuen Standortsuche ist es jetzt einfach, einen Standort auf der Karte zu finden. Wenn man einen Ortsnamen eingibt, wird dieser an die Nominatim API gesendet, um in der OSM Datenbank zu suchen. Die Ergebnisse werden in einer Tabelle aufgelistet. Während man mit der Maus über eine Ergebniszeile fährt, wird das jeweilige Feature auf der Karte hervorgehoben. Ein einfacher Klick auf die Zeile zentriert die Karte auf das Feature und ein Doppelklick vergrößert die Karte.
Es ist auch möglich, Koordinaten im Lat/Long-Format durch Komma getrennt einzugeben. Auf diese Weise wird die Nominatim API umgekehrt verwendet und Orte, die sich an diesen Koordinaten befinden, werden aufgelistet. Ein Beispiel für eine Anfrage, die der Client an die API sendet, sieht so aus:

https://nominatim.openstreetmap.org/reverse?
q=&format=json&limit=100&viewboxlbrt=&bounded=1&polygon_text=1&lat=50.7409291&lon=7.0968584&page=1&start=0&maxFeatures=25

 

Kompositionen-Schieberegler

Wenn maps.mundialis initial lädt, zeigt es den Schieberegler, mit dem man verschiedene Farbkompositionen von bereits prozessierten Sentinel-2 Szenen auswählen kann. Jetzt ist links daneben ein Button, mit dem man den Schieberegler wechseln kann. Der neue Schieberegler erlaubt es, von zwei Kompositionen, die man über eine Kombobox auswählen kann, die Schnittkante zu verschieben. So kann man direkt z.B. eine Echtfarben-Komposition mit einer Infrarot-Falschfarben-Komposition vergleichen.

Sentinel-2 Datumsauswahl


Vor dieser Aktualisierung hatten wir einen Schieberegler benutzt, um Anfangs- und Enddatum auszuwählen, mit denen Sentinel-2 Szenen-Umringe selektiert wurden. Im Laufe der letzten Monate kamen immer mehr auswählbare Zeitstempel zusammen – damit wurde es wird immer schwieriger, auf dem Schieberegler das gewünschte Datum auszuwählen. Deshalb haben wir diesen Schieberegler mit Datumsfeldern ausgetauscht, um eine komfortable Bedienung zu ermöglichen.

Hinzufügen von WMS


Diese Komponente basiert auf einer existierenden Komponente aus der BasiGX Software von terrestris. Sie ermöglicht es, einen Web Map Service bei Laufzeit in die Applikation zu laden. Weil hierfür momentan kein Proxy benutzt wird können nur externe Dienste eingebunden werden, die Cross-Origin Resource Sharing (CORS) erlauben.

Neu strukturiertes Menü


Mit immer mehr alten und neuen Funktionen wuchs das Menü und wurde unübersichtlich. Mit dieser Aktualisierung haben wir aufgeräumt und die Menüpunkte in „Background Layers“ (Hintergrunddienste), „Foreground Layers“ (Layer anzeigen), „Map Tools“ (Kartenwerkzeuge), „Location Search“ (Ortssuche) und „Fullscreen“ (Vollbild) aufgeteilt.

Vollständig berechnete Farb-Komposite

Früher wurden einige Sentinel-2-Szenen in einigen Farbkompositen und andere in allen fünf verschiedenen Kompositen verarbeitet, die in der Anwendung angezeigt werden. Dies ermöglichte ein sanftes Gleiten zwischen den Kompositen für die meisten Szenen, aber nicht für alle. Dann war die Szene nicht sichtbar, wenn man zu einem bestimmten Komposit rutschte, nur weil es nicht prozessiert wurde. Deshalb haben wir viele Szenen nachbearbeitet, so dass sie in allen fünf Farbkompositen zur Verfügung stehen. Die Nachbearbeitung ist noch im Gange, aber für die meisten Szenen ist sie bereits abgeschlossen, so dass ein reibungsloses Gleiten zwischen den Farbkompositen möglich ist.

Korrekte Anzeige von Features, die die Datumsgrenze überschreiten

In der früheren Version von maps.mundialis.de bereitete die Datumsgrenze (wo sich 180W und 180E „treffen“) Probleme bei kleinen Sentinel-2 Polygonen, die oben drauf lagen. Im Wesentlichen wurde das Polygon auf der langen (falschen) Seite um den Globus gewickelt, anstatt richtig dargestellt zu werden (siehe Abbildung unten). Wir mussten eine Lösung für dieses bekannte Problem implementieren. Weiter unten kann man weiterlesen, wie wir die Lösung gefunden haben, oder auch direkt zur Lösung springen.

Um die verfügbaren Sentinel-2-Szenen anzeigen zu können, fragen wir die Metadaten einer Szene mitsamt ihrer Footprint-Anzeige ab und speichern sie in einer PostGIS Datenbank. Dann veröffentlichen wir sie mit einer GeoServer Instanz als WMS zur Anzeige und Abfrage in maps.mundialis. Belässt man es dabei, würden einige Features, die die Datumsgrenze kreuzen, sie nicht wirklich überqueren, sondern umdrehen, um die über die ganze Welt zu verlaufen, um alle Punkte des Polygon zu zeichnen. Das sieht so aus:

Auf den ersten oberflächlichen Blick könnte das Problem mit Setzen von openlayers‘ „wrapX“ auf false gelöst werden, aber das verhindert nur, dass ein Feature über „mehrere Welten“ gezeichnet wird:

Also mussten wir tiefer graben. Wie wäre es, wenn das Feature an der Datumsgrenze geteilt und als Multipolygon gerendert wird? Auf diese Weise muss der Client nicht mit der Überschreitung der Datumsgrenze umgehen.

Um dies zu testen, haben wir uns mit der Datenbank verbunden, in der die Features gespeichert sind und ein Polygon ausgewählt, das die Datumsgrenze kreuzt. Dann versuchten wir es mit der PostGIS-Funktion ST_Split:

SELECT
    ST_AsText(
        ST_Split(
            geom,
            ST_SetSRID(
                ST_MakeLine(
                    ST_MakePoint(180,-90),
                    ST_MakePoint(180,90)
                ),
                4326
            )
        )
    )
FROM wkt
WHERE uuid='90981983-044d-44e4-80df-7f5c3252b353'
;

Das Abfrageergebnis ist eine GeometryCollection mit zwei Polygonen:

GEOMETRYCOLLECTION(POLYGON((180 -17.1616365585413,179.240749367346 -17.161645777118,179.241003668701 -17.1494367705579,180 -17.1616365585413)),POLYGON((180 -18.1529867753605,-179.743259878855 -18.1573436632797,-179.728268070337 -17.1660042633555,179.240749367346 -17.161645777118,179.220376301768 -18.1397565367159,180 -18.1529867753605)))

Das Ergebnis konnten wir schnell mit einem QGIS Plugin namens QuickWKT überprüfen, das es erlaubt, Well Known Text (WKT) einzufügen. Nach dem Einfügen unseres WKT, das wir vom obigen Ergebnis der PostGIS-Abfrage erhalten haben, stellten wir fest, dass wir noch nicht ganz fertig sind:

Das Polygon hat sich an der Datumsgrenze geteilt, aber es wird immer noch nicht umgebrochen. Daraus ergibt sich die Frage: Weiß das Feature selbst, wie es gezeichnet werden soll? Einige Untersuchungen später fanden wir heraus, dass PostGIS selbst ahnungslos ist, aber eine Funktion zur Handhabung solcher Fälle hat: ST_ShiftLongitude.

Bringen wir also dem Feature bei, die Datumsgrenze zu überschreiten!

SELECT
    ST_ShiftLongitude(geom)
FROM wkt
WHERE uuid='90981983-044d-44e4-80df-7f5c3252b353'
;

Abfrageergebnis:

0103000020E610000001000000080000000000000000806640D0C57403612931C0DAE52807B2886640C41B62417F2A31C06FFC0E3737886640F104A0AC472832C00000000000806640D1E42C242A2732C0071D9A520D676640DE989A15C72332C02F69534DB6676640B05DFA7C412631C00000000000806640D0C57403612931C00000000000806640D0C57403612931C0

Und tatsächlich „weiß“ es das Feature jetzt, wie wir aus einer Vorschau in QGIS mit QuickWKT sehen können:

Dies funktioniert auf den ersten Blick. Es funktioniert aber nicht mehr, wenn wir nach Web Mercator reprojizieren, was wir in maps.mundialis verwenden. Wieso nicht? Die Beschreibung von ST_ShiftLongitude sagt aus, dass die Funktion:

„Reads every point/vertex in every component of every feature in a geometry, and if the longitude coordinate is <0, adds 360 to it. The result would be a 0-360 version of the data to be plotted in a 180 centric map“

Übersetzt bedeutet das, dass die Funktion jeden Punkt/Wertex in jeder Komponente jedes Features in eine Geometrie liest, und wenn die Längengradkoordinate So verliert sie ihre Wirkung in der Web-Mercator-Projektion.

Können wir die beiden Funktionen ST_Split und ST_ShiftLongitude kombinieren?

SELECT
    ST_AsText(
        ST_Split(
            ST_ShiftLongitude(geom),
            ST_SetSRID(
                ST_MakeLine(
                    ST_MakePoint(180,-90),
                    ST_MakePoint(180,90)
                ),
                4326
            )
        )
    )
FROM wkt
WHERE uuid='90981983-044d-44e4-80df-7f5c3252b353'
;

Abfrageergebnis:

GEOMETRYCOLLECTION(POLYGON((180 -17.1616365585413,180.271731929663 -17.1660042633555,180.256740121145 -18.1573436632797,180 -18.1529867753605,180 -17.1616365585413)),POLYGON((180 -18.1529867753605,179.220376301768 -18.1397565367159,179.241003668701 -17.1494367705579,180 -17.1616365585413,180 -18.1529867753605)))

Die Vorschau in QGIS mit QuickWKT sieht vielversprechend aus:

Wir erstellen eine SQL-Ansicht im GeoServer, um das Ergebnis zu sehen:

SELECT
    id,
    ST_Split(
        ST_ShiftLongitude(geom),
        ST_SetSRID(
            ST_MakeLine(
                ST_MakePoint(180,-90),
                ST_MakePoint(180,90)
            ),
            4326
        )
    )
FROM wkt

Bei der Abfrage des Layers erhalten wir eine Exception mit der Meldung „Relate Operation called with a LWGEOMCOLLECTION type. This is unsupported“. Dies können wir mit der PostGIS-Funktion ST_CollectionExtract lösen, die ein MultiPolygon zurückgibt:

SELECT
    ST_AsText(
        ST_CollectionExtract(
            ST_AsText(
                ST_Split(
                    ST_ShiftLongitude(geom),
                    ST_SetSRID(
                        ST_MakeLine(
                            ST_MakePoint(180,-90),
                            ST_MakePoint(180,90)
                        ),
                        4326
                    )
                )
            ),
            3
        )
    )
FROM wkt
WHERE uuid='90981983-044d-44e4-80df-7f5c3252b353'

Dieses scheint schließlich zu funktionieren. Aber Moment mal. Was ist denn das?

Wir haben ST_ShiftLongitude auf alle Features angewendet. Auf diese Weise werden die Features, die tatsächlich die Datumsgrenze überschreiten, korrekt geteilt, aber wir verschieben das Problem buchstäblich auf Features, die den 0° Längengrad überschreiten. Wir müssen einen Weg finden, um die Polygone zu unterscheiden, die die Datumsgrenze überschreiten.

Bei der Transformation von nur denjenigen Polygonen, auf die die folgende WHERE-Klausel zutrifft, finden wir nicht alle betroffenen Polygone, da ohne ST_ShiftLongitude einige den 180°-Längengrad nicht wirklich schneiden (siehe erste Versuche mit ST_Split oben).

WHERE
    ST_Intersects(
        geom,
        ST_SetSRID(
            ST_MakeLine(
                ST_MakePoint(180,-90),
                ST_MakePoint(180,90)
            ),
            4326
        )
    )

Wenn wir ST_ShiftLongitude verwenden und ST_Intersects anwenden, finden wir nicht nur die korrekt transformierten Features, die die Datumsgrenzen überschreiten, sondern auch diejenigen, die den 0° Längengrad überschreiten, weil sie nun verschoben sind und auch den 180° Längengrad überqueren (man beachte, dass ST_ShiftLongitude 360 addiert, wenn der Längengrad kleiner als 0° ist).

WHERE
    ST_Intersects(
        ST_ShiftLongitude(geom),
        ST_SetSRID(
            ST_MakeLine(
                ST_MakePoint(180,-90),
                ST_MakePoint(180,90)
            ),
            4326
        )
    )

Wir müssen mehr als eine Bedingung definieren. Da eine Sentinel-2-Szene weniger als die Hälfte der Welt abdeckt (wenn man S2-Daten kennt, kann man sich einfach mal vorstellen, dass sie es würde…), können wir die betroffenen Features mit der folgenden WHERE-Klausel sicher auswählen:

WHERE
(
    ST_Intersects(
        geom,
        ST_SetSRID(
            ST_MakeLine(
                ST_MakePoint(-90,-90),
                ST_MakePoint(-90,90)
            ),
            4326
        )
    )
    AND
    ST_Intersects(
        geom,
        ST_SetSRID(
            ST_MakeLine(
                ST_MakePoint(90,-90),
                ST_MakePoint(90,90)
            ),
            4326
        )
    )
)

So können wir endlich unsere SQL-Abfrage zusammenstellen, die wir direkt in einer Materialized View speichern.

Die Lösung

Die Lösung besteht darin, aus unserer Tabelle eine Materialized View mit den betroffenen Features zu erzeugen. Diese müssen verschoben und geteilt werden, um korrekt dargestellt zu werden.


CREATE MATERIALIZED VIEW wkt_shiftdateline AS

    -- Splits geometries which cross the dateline (these are stretched around
    -- the whole world, so intersect -90 and 90 longitude).
    -- To split, shift to "next world" (ST_ShiftLongitude adds 360 if value is
    -- below 0) and split at dateline.
    -- Transform to MultiPolygon with ST_CollectionExtract and set SRID
    -- to avoid geometry errors.

    SELECT
        id,
        name,
        uuid,
        size,
        cloud,
        cast(substring(name from '...............$') as timestamp) as date,
        ST_SetSRID(
            ST_CollectionExtract(
                ST_AsText(
                    ST_Split(
                        ST_ShiftLongitude(geom),
                        ST_SetSRID(
                            ST_MakeLine(
                                ST_MakePoint(180,-90),
                                ST_MakePoint(180,90)
                            ),
                            4326
                        )
                    )
                ),
                3
            ),
            4326
        ) geom
    FROM wkt
    WHERE (
        ST_Intersects(
            geom,
            ST_SetSRID(
                ST_MakeLine(
                    ST_MakePoint(-90,-90),
                    ST_MakePoint(-90,90)
                ),
                4326
            )
        )
        AND
        ST_Intersects(
            geom,
            ST_SetSRID(
                ST_MakeLine(
                    ST_MakePoint(90,-90),
                    ST_MakePoint(90,90)
                ),
                4326
            )
        )
    )
    AND
        geom IS NOT NULL

    UNION

    -- Leaves geometries as they are.
    -- Selects all geometries which do not intersect both -90 AND 90 longitude.

    SELECT
        id,
        name,
        uuid,
        size,
        cloud,
        cast(substring(name from '...............$') as timestamp) as date,
        geom
    FROM wkt
    WHERE (
        uuid NOT IN (
            SELECT uuid
            FROM wkt
            WHERE
                ST_Intersects(
                    geom,
                    ST_SetSRID(
                        ST_MakeLine(
                            ST_MakePoint(-90,-90),
                            ST_MakePoint(-90,90)
                        ),
                        4326
                    )
                )
                AND
                ST_Intersects(
                    geom,
                    ST_SetSRID(
                        ST_MakeLine(
                            ST_MakePoint(90,-90),
                            ST_MakePoint(90,90)
                        ),
                        4326
                    )
                )
            )
    )
    AND
        geom IS NOT NULL

    --the order by is applied to the complete resultset
    ORDER BY date ASC

;

Und jetzt veröffentlichen wir die Materialized View im GeoServer als Layer, um sie in unserer Anwendung zu verwenden. Unten sieht man einen Vergleich zwischen „vorher“ und „nachher“:

Kontakt
terrestris
Github Linkedin Xing Twitter