Menu
  • About us
    • Company
    • Team
    • Career
    • References
    • Partners
    • Contact
  • Products
    • Earth Observation metadata enhancer
    • Temperaturdaten
    • LST-Data Client
    • Trainings
  • Services
    • actinia – geoprocessing in the cloud
    • Geospatial and EO data analysis
    • CORONA spy satellite data
    • maps.mundialis
    • Web Map Services
    • GIS Development
    • Open Source GIS
  • Markets
    • Remote Sensing for Agriculture
    • Satellite Images for Forest and Land Cover Restoration
    • Glass fibre route planning – FTTH
    • Copernicus and Sentinel
  • News & Blog
    • News
    • Blog
    • Satellite image of the month
  • English English
  • Deutsch Deutsch
  • About us
    • Company
    • Team
    • Career
    • References
    • Partners
    • Contact
  • Products
    • Earth Observation metadata enhancer
    • Temperaturdaten
    • LST-Data Client
    • Trainings
  • Services
    • actinia – geoprocessing in the cloud
    • Geospatial and EO data analysis
    • CORONA spy satellite data
    • maps.mundialis
    • Web Map Services
    • GIS Development
    • Open Source GIS
  • Markets
    • Remote Sensing for Agriculture
    • Satellite Images for Forest and Land Cover Restoration
    • Glass fibre route planning – FTTH
    • Copernicus and Sentinel
  • News & Blog
    • News
    • Blog
    • Satellite image of the month
  • English English
  • Deutsch Deutsch

Update for our maps.mundialis application solves dateline wrap issues


6. September 2017 | Category blog, General

The update of our web application maps.mundialis.de includes new features which are described subsequently. To directly jump to the section which handles the dateline wrap click here.

 

New Features

  • Location Search
  • Composite Swipe slider
  • Sentinel-2 datepicker
  • Add WMS
  • Restructured menu
  • Complete composites
  • Display features correctly which wrap the dateline

 

Location Search


With the new location search it is now easy to find a location on the map. If you enter a placename it is send to the Nominatim API to search in the OSM database. Suitable results are listed in a table. While hovering with the mouse over a result row the respective feature will be highlighted on the map. A single click on the row will center the map around the feature and a double click will also zoom in.
It is also possible to enter coordinates in lat/long format separated by a comma. This way the Nominatim API is used reverse and locations which are at these coordinates are listed. An example request which the client sends to the API looks like this:

1
2
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

 

Composite Swipe slider

When maps.mundialis is initially loaded into your browser, it shows the slider with which you can choose the color composite from already processed Sentinel-2 scenes. Now there is a button on the left side which lets you switch the slider. The new slider enables you to swipe between two different composites which you can select in combo boxes. This way you can directly compare e.g. a true color composite with a color infrared composite.

 

Sentinel-2 datepicker


Prior to this update we used a slider to choose the start-date and end-date to select Sentinel-2 footprints. As time goes by the dates which can be selected from the slider also grow and it becomes hard to select an exact date. So we exchanged the slider with a datepicker for a comfortable usage.

 

Add WMS


This component is added using an existing component from BasiGX. It allows to load a Web Map Service into the application on the fly. As it doesn’t use a proxy at the moment, only external services can be used which allow Cross-Origin Resource Sharing (CORS).

 

Restructured menu


With more and more old and new features the menu grew and became confusing. With this update we tidied up the menu so that it is now structured into “Background Layers”, “Foreground Layers”, “Map Tools”, “Location Search” and “Fullscreen”.

 

Complete composites

Earlier some Sentinel-2 scenes were processed in some color composites and others in all five different composites which are displayed in the application. This allowed a smooth sliding between composites for most scenes but not for all of them. Then the scene was not visible when sliding to a certain composite simply because it was not processed. Therefore we reprocessed many scenes so they are available in all five color composites. The reprocessing is still ongoing but complete for most scenes allowing a smooth sliding between the color composites.

 

Display features correctly which wrap the dateline

In the earlier version of maps.mundialis.de the dateline (where 180W and 180E “meet”) caused problems in case of small Sentinel-2 tiles lying on top of it. Essentially the tile was wrapped around the globe on the long (wrong) side rather than being displayed properly (see figure below). We had to implement a solution for this well known issue. Read on how we got there or read directly how we solved it.

To display available Sentinel-2 scenes, we request the metadata of a scene along with its footprint and store it in a PostGIS database. Then we publish it with a GeoServer instance as a WMS to display and query it in maps.mundialis. While leaving it like that, some features which cross the dateline wouldn’t actually cross it but turn around to stretch over the whole world to draw all points the polygon contains. This looks like this:


On a first superficial glance the problem might be solved using openlayers’ “wrapX” set to false but this only prevents a feature from being drawn over “multiple worlds”:


So we needed to dig deeper. How about cutting the feature at the dateline and render it as a multipolygon? This way the client doesn’t need to handle the dateline wrapping.

To test this we connected to the database where the features are stored and picked a polygon which crosses the dateline. Then we tried with the PostGIS function ST_Split:

PgSQL
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
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'
;

The query result is a GeometryCollection with two Polygons:

PgSQL
1
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)))

We could quickly check the result with a QGIS Plugin called QuickWKT which allows to paste Well Known Text (WKT). After pasting our WKT which we received from the PostGIS query result from above, we saw that we are not quite there:

The Polygon did split at the dateline, but still doesn’t wrap around. This leads to the question: does the feature itself know how it should be drawn? Some research later we found out that PostGIS itself is clueless but has a function to handle such cases: ST_ShiftLongitude.


So let’s teach the feature to cross the dateline!

PgSQL
1
2
3
4
5
SELECT
    ST_ShiftLongitude(geom)
FROM wkt
WHERE uuid='90981983-044d-44e4-80df-7f5c3252b353'
;

Query result:

PgSQL
1
0103000020E610000001000000080000000000000000806640D0C57403612931C0DAE52807B2886640C41B62417F2A31C06FFC0E3737886640F104A0AC472832C00000000000806640D1E42C242A2732C0071D9A520D676640DE989A15C72332C02F69534DB6676640B05DFA7C412631C00000000000806640D0C57403612931C00000000000806640D0C57403612931C0

And as a matter of fact, the feature now knows as we can see from a preview in QGIS with QuickWKT:

This works at first sight but is no longer working when reprojecting to Web Mercator which we use in maps.mundialis. How come? The description of ST_ShiftLongitude tells that the function:

“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”

So it looses its effect in Web Mercator projection.


Can we combine the two functions ST_Split and ST_ShiftLongitude?

PgSQL
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
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'
;

Query result:

PgSQL
1
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)))

The preview in QGIS with QuickWKT looks promising:


We create a SQL-view in GeoServer to see the result:

PgSQL
1
2
3
4
5
6
7
8
9
10
11
12
13
SELECT
    id,
    ST_Split(
        ST_ShiftLongitude(geom),
        ST_SetSRID(
            ST_MakeLine(
                ST_MakePoint(180,-90),
                ST_MakePoint(180,90)
            ),
            4326
        )
    )
FROM wkt

When requesting the layer we get an exception containing the message “Relate Operation called with a LWGEOMCOLLECTION type. This is unsupported”. We can solve this by using the PostGIS function ST_CollectionExtract which returns a MultiPolygon:

PgSQL
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
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'

This finally seem to work. But wait. What is this?


We applied ST_ShiftLongitude to all features. This way the features, which actually cross the dateline, are split correctly but we literally shift the problem to features which cross the 0° longitude. We need to find a way to distinguish between the polygons which should wrap around the dateline and the others.

When transforming only the ones to which the following WHERE clause applies, we don’t catch all affected polygons because without ST_ShiftLongitude some really don’t intersect the 180° longitude (see first attempts with ST_Split above).

PgSQL
1
2
3
4
5
6
7
8
9
10
11
12
WHERE
    ST_Intersects(
        geom,
        ST_SetSRID(
            ST_MakeLine(
                ST_MakePoint(180,-90),
                ST_MakePoint(180,90)
            ),
            4326
        )
    )
 

When we use ST_ShiftLongitude and try the intersection, we do not only catch the nicely transformed ones which wrap around the dateline but also the ones which wrap around the 0° longitude because now they are shifted and cross the 180° longitude as well (remember that ST_ShiftLongitude adds 360 if the longitude coordinate is below 0.

PgSQL
1
2
3
4
5
6
7
8
9
10
11
12
WHERE
    ST_Intersects(
        ST_ShiftLongitude(geom),
        ST_SetSRID(
            ST_MakeLine(
                ST_MakePoint(180,-90),
                ST_MakePoint(180,90)
            ),
            4326
        )
    )
 

We need to define more than one condition. As one Sentinel-2 scene covers less than half of the world (if you are familiar with S2 data just image it would) we can safely pick the affected features with the following WHERE clause:

PgSQL
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
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 finally we can put together our SQL-query which we directly store in a materialized view.

 

The solution

The solution is to create a materialized view out of our table containing the affected features. These need to be shifted and split to be displayed correctly.

PgSQL
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
 
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
 
;
 

And now we publish the materialized view in GeoServer as a layer to use it in our application. Below you can see a comparison between “before” and “after”:


Leave a Reply Cancel reply

Your email address will not be published. Required fields are marked *

Follow Us

Contact

mundialis GmbH & Co. KG
Kölnstrasse 99
53111 Bonn

Phone: +49 (0)228 / 387 580 80
Fax: +49 (0)228 / 962 899 57

E-Mail: info [at] mundialis [dot] de

Latest News

  • ++ Eine neue Heimat für Actinia ++

    Tuesday January 24th, 2023
  • Positive towards the future – Management changes at mundialis

    Tuesday October 4th, 2022
  • Recap of FOSS4G 2022 in Florence

    Friday September 16th, 2022
  • FOSS4G 2022 in “La Bella” – Greetings from Florence!

    Thursday July 14th, 2022
  • VALE project successfully completed

    Tuesday June 21st, 2022
View all News

Blog

  • Satellite image of the month – January – Yalu River (People’s Republic of China and Democratic People’s Republic of Korea) 2. January 2023
    Yalu River – People’s Republic of China and Democratic People’s ...
  • Satellite image of the month – December – Super Pit gold mine (Australia) 1. December 2022
    Super Pit gold mine – Australia, recorded by the Sentinel-2A ...
  • Satellite image of the month – November – Kufra-Oases (Libya) 1. November 2022
    Kufra-Oases – Libya, recorded by the Sentinel-2A satellite on March ...
View all Blog Posts

Copyright © 2015-2021 mundialis GmbH & Co. KG

Imprint Privacy

Theme created by PWT. Powered by WordPress.org