Update for our maps.mundialis application solves dateline wrap issues
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:
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:
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!
1 2 3 4 5 |
SELECT ST_ShiftLongitude(geom) FROM wkt WHERE uuid='90981983-044d-44e4-80df-7f5c3252b353' ; |
Query result:
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?
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:
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:
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:
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).
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.
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:
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.
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”: