gravitystorm / openstreetmap-carto

A general-purpose OpenStreetMap mapnik style, in CartoCSS
Other
1.53k stars 819 forks source link

Virtual borders on antimeridian #3504

Open kocio-pl opened 5 years ago

kocio-pl commented 5 years ago

There are virtual borders on antimeridian, which of course should not be there:

https://www.openstreetmap.org/way/239992336#map=6/67.306/-177.726

screenshot_2018-11-09 openstreetmap 1

https://www.openstreetmap.org/way/42394837#map=9/51.9933/179.8860

screenshot_2018-11-09 openstreetmap

kocio-pl commented 5 years ago

This problem is visible also on HOT layer, but is not present on French fork:

http://tile.openstreetmap.fr/?zoom=9&lat=51.87829&lon=-180.54689&layers=B000000FFFFFF https://github.com/cquest/osmfr-cartocss

screenshot_2018-11-09 tile openstreetmap fr

jeisenbe commented 4 years ago

This bug is still visible at https://www.openstreetmap.org/#map=9/51.9933/179.8860

@imagico - I recall you mentioned a way to solve a similar problem for archipelagos. How can we fix this problem for admin boundary relations?

pitdicker commented 3 years ago

Isn't it the case that all data is split up on the antimeridian? I see roads and rivers that don't even visually attach when crossing the boundary.

mboeringa commented 3 years ago

Isn't it the case that all data is split up on the antimeridian?

Yes, in most cases. However, sometimes multipolygon relations are defined that "connect" the two parts in an "administrative" way. However, these MP relations cause issues with many GIS systems and PostGIS functions designed for planar coordinate systems, as there is no inherent knowledge of the antimeridian in a planar coordinate system like Web Mercator.

E.g. if you try to calculate a bounding box for such a multipolygon spanning the antimeridian, it will cover the globe: https://github.com/openmaptiles/openmaptiles/issues/1122#issuecomment-913147529

pitdicker commented 3 years ago

Then I wonder if there is any solution possible here except pre-processing admin boundaries (and maybe other geometry).

kocio-pl commented 3 years ago

Does anyone know how the French fork deals with it?

pitdicker commented 3 years ago

That would be this part of their style? https://github.com/cquest/osmfr-cartocss/blob/master/osmfr.yml#L1139:L1160

TheRealJackUK commented 2 years ago

It also affects routing in OSMAnd and any other navigation using OSM, there is a primary road which crosses over this line, navigation will not work over the line, the road is an important road for truckers.

sommerluk commented 2 years ago

I cannot reproduce this issue locally (using Kosmtik) neither for our last release tag 5.6.1 (117df85a0897e4d75ee6a3291dbaaae4703467de) nor for the current master (56a7a8fecdd62265e4d4fa799c7b6d516790cd3b). On my machine, it renders correctly, without these artefacts.

sommerluk commented 1 year ago

So if I understand correctly, this is the status quo:

OSM data model and mapping praxis

The OSM database defines limits for the coordinate space of nodes.

So, lat 0°, lon −180° and lat 0°, lon 180° represent the same point on earth, but are different points in the coordinate space. OSM ways cannot go beyond the coordinate space. Therefore, OSM ways cannot cross the antimeridian at ±180°. In recent discussions about an evolution of the OSM data model, like in Jochen Topf’s paper, the problem with antimeridian are mentioned, but considered as unimportant, so future improvements are not likely. The current common mapping praxis is:

Software support

The problem extends beyond the OSM database: It seems that most GIS software does not care much about proper antimeridian support. This includes Mapnik. (There were plans for antimeridian support in Mapnik within a GSOC project that apparently has not been realized.) Also the “Nearby objects” search at openstreetmap.org does not work across antimeridian.

Rendering at openstreetmap.org

So what are the results? A good test place are the Fiji Islands. At openstreetmap.org we see the borders at antimeridian because these borders actually are in the database. Also note that the country name “Vity” is not centered within the whole country, but centered within the bigger (left) part of the multipolygon: grafik

The twice-mapped Dateline monument mentioned earlier renders fine at z19: grafik

It renders bad at z17, I suppose because when rendering, the overlap data around the actually rendered tile, which normally is included in the rendering process to avoid artefacts on tile borders, is not available across the antimeridian: grafik

This leads also to artefacts on roads rendering at places like https://www.openstreetmap.org/#map=19/-16.79356/180.00000 where the roads do not join up correctly because of the lack of the round border. This is more visible, the bigger the angle at this joint. grafik Here for demonstration propose a local rendering of the same place, with a red road and line-cap: butt; line-join: bevel;, making the problem more obvious: demo1

Local rendering with Kosmtik

Local rendering with Kosmtik works just as at openstreetmap.org if your OSM data is complete. Note that popular database extracts from Geofabrik and download.openstreetmap.fr do not handle this will. See here for appropriate test data.

imagico commented 1 year ago

The core of the issue lies in the fact that almost no software working with geodata is actually able to work with geographic coordinates but instead treats them as coordinates in equirectangular projection. This has lead to mapping in OSM being done in equirectangular projection as well.

There are in principle mainly two solutions to this issue:

  1. Render tiles near the 180 degree meridian in a different projection with a shifted central meridian. The generated tiles then could be shifted back to the right place. Since in the shifted coordinate system is continuous across the 180 degree meridian this would avoid the issues.
  2. Extend the data beyond the 180 degree meridian with a copy of the data from the other side shifted by +/-360 degrees. This is more work on the data processing side but can use existing software as is - as long as it is able to deal with coordinates beyond +/-180 degrees (or beyond +/- 20037508.34 X-Coordinate in Mercator) without troubles.

Both approaches will require elaborate and error prone merging and de-duplication at the 180 degree meridian to convert the mapping in equirectangular projection into a projection with a different central meridian. This would be easier if mapping was actually in geographic coordinates and not in equirectangular projection.

To illustrate how this works practically here a short demo for method 1, first with a split geometry (like in OSM):

echo '{"type": "FeatureCollection","crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },"features": [{ "type": "Feature", "geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ -175.0, 2.0 ], [ -175.0, -1.0 ], [ -180.0, -2.0 ], [ -180.0, 1.0 ], [ -175.0, 2.0 ] ] ], [ [ [ 180.0, 1.0 ], [ 180.0, -2.0 ], [ 175.0, -1.0 ], [ 175.0, 2.0 ], [ 180.0, 1.0 ] ] ] ] } }]}' > split.geojson
ogr2ogr -a_srs "EPSG:4326" -t_srs "+proj=latlong +datum=WGS84" -s_srs "+proj=latlong +datum=WGS84 +pm=180" -f GEOJSON split2.geojson split.geojson
ogr2ogr -f GEOJSON -dialect SQLITE -sql "SELECT ST_Union(GEOMETRY) FROM split" merged.geojson split2.geojson

And then the simpler version with an unsplit geometry plainly mapped in geographic coordinates (the data will look like a large globe spanning polygon in equirectangular projection - like when you open it in QGIS, but is actually a smaller polygon with 10 degree extent in longitude):

echo '{"type": "FeatureCollection","features": [{ "type": "Feature", "geometry": { "type": "Polygon", "coordinates": [ [ [ -175.0, 2.0 ], [ 175.0, 1.0 ], [ 175.0, -2.0 ], [ -175.0, -1.0 ], [ -175.0, 2.0 ] ] ] } }]}' > unsplit.json
ogr2ogr -a_srs "EPSG:4326" -t_srs "+proj=latlong +datum=WGS84" -s_srs "+proj=latlong +datum=WGS84 +pm=180" -f GEOJSON unsplit2.geojson unsplit.geojson

The problem is of course that implementing either of these solution without builtin support in our toolchain is hard.

What you observe in kosmtik is probably that some clipping of geometries is done there that removes anything directly on the 180 degree meridian. That is a fairly simple method to remove the most obvious artefacts but as you already indicate it causes problems with any real geometries directly on the 180 degree meridian or with styling that depend on continuity of geometries across it.

sommerluk commented 1 year ago

The different local rendering that I was seeing was indeed an issue of the database extract I was using (Geofabrik). Things at 180° or -180° longitude are simply missing.

The database extracts of Fiji from Geofabrik and also those from download.openstreetmap.fr both lack some data on antimeridian.

Also for customized extracts with https://extract.bbbike.org/ and https://protomaps.com/downloads/osm it seems difficult to get data around the antimeridian

Therefore, I've downloaded with JOSM a small region around the antimeridian at Fiji. Yes, this works (though JOSM issue an error message and does not mark in the file the downloaded region at the other side of the antimeridian as downloaded region). Here, it is bundled with the extract from Geofabrik: Both files together have complete data for Fiji, including antimeridian, and can therefore be used for local testing:

fiji.zip

sommerluk commented 1 year ago

@imagico What you propose is a high-quality solution. I like the idea.

I suppose it isn’t obvious that we can implement this high-quality solution in the short term. But perhaps, in short term, we can do something more simple.

It seems to be consistent mapping praxis to construct administrative boundary relations so that the part that is on the anti-meridian is represented by different ways than the “actual” part: ways that have a longitude of either -180° or 180° for all their nodes and do not extend into the “normal” space. All these ways are tagged with closure_segment=yes. This applies to Fiji, Tuvalu, USA, Russia. This is also mentioned in the wiki.. However, for this protected-area boundary crossing the antimeridian this is apparently not the case.

At mentioned earlier in the discussion, neither the French style nor the Hot style show these boundary artefacts. How do they do that?

The French style relies on the ways with boundary=administrative tag instead of the boundary relation that we use in our style. At low zoom, instead of WHERE osm_id < 0 it uses WHERE osm_id > 0. And combined with the above-mentioned tagging practice and the default WHERE way && !bbox! this eliminates the border. And on high zoom, I think it works similar (though the query is more complex and involves relations, it seems to still rely basically on the OSM ways themselves). Note that the border would become visible again if a single way has one part on the antimeridian, but also another part in the “normal” space. It renders without artefacts (note that the antimeridian is near to the cursor): Screenshot_20220927_145614

The Hot style does differently. It calculates an intersection: Instead of ( SELECT way, it uses ( SELECT ST_Intersection(way, !bbox!) AS way,. That means also that the boundary dash-array would “restart” at each meta-tile boundary. And it leads to small artefacts at the antimeridian visible from z18 and higher: Screenshot_20220927_145720

What could we do? Maybe for mapper feedback and for clarity, the Hot style approach is better. Its two rendering artefacts could both be avoided by using this code instead:

(SELECT
    ST_Intersection(
        way,
        ST_MakeEnvelope(-20037508.342789242, -20037508.342789244, 20037508.342789242, 20037508.342789244, 3857)
    ) AS way,
    -- …

The floating point numbers used in this snippet come from this comment and I have reduced the longitude step by step until it eliminates the boundary rendering on the antimeridian, while leaving the latitude as-is. This renders almost without visible artefacts, even on z21: Screenshot_20220927_150625

mboeringa commented 1 year ago

@sommerluk ,

The work on the new osm2pgsql 'flex' style for openstreetmap-carto by Paul Norman may be relevant here as well:

https://github.com/gravitystorm/openstreetmap-carto/pull/4431

It will feature a new 'planet_osm_admin' Line geometry table, which allows rendering administrative boundaries based on Line geometry instead of Polygon geometry. There may be more opportunities to fix this using that new table, possibly in combination with further adjustments to what is included in 'planet_osm_admin' by adjusting the 'flex' style to deal with these issues.

@pnorman may be able to give some more specific comments.

mboeringa commented 1 year ago

One other prove that changes to the Lua 'flex' style in development may help here, is that in my personal adjusted version of that osm2pgsql Lua flex style of Paul, I already skip in Lua any ways that are tagged 'natural=coastline', so they do not end up in the 'planet_osm_admin' table.

That leads to the following result (blue lines), where the antimeridian is visible for the parts over sea (that don't have 'natural=coastline' tagged), but not for the land areas that have a combination of 'closure_segment=yes' and 'natural=coastline' tagged.

In fact, this now gives me the new idea that I should also check for the 'closure_segment=yes' tag, and skip adding those as well for the new 'planet_osm_admin' table in my variant of Paul's work.

This all requires 'stage 2' osm2pgsql flex processing by the way.

afbeelding

mboeringa commented 1 year ago

@sommerluk,

After adjusting my personal variant of the carto style to exclude 'closure_segment=yes' tagged features, I get this result based on the 'planet_osm_admin' Line geometry table that solves the boundary rendering on the anti-meridian. The dark blue line is the relevant line in the 'planet_osm_admin' table for 'admin_level=2', note that the anti-meridian part of the same 'admin_level=2' boundary polygon, is now correctly omitted.

Note that the data of the 'planet_osm_admin' table though, cannot be easily used for boundary 'name' labeling, parallel to the line of the boundary, as the information of which country is on each side of the border is lost.

This could potentially be solved by using the Polygons for labeling only, and setting symbology invisible for them, and use the Lines for the display of the boundary. However, that may cause "orphaned" boundary name labels for the 'closure_segment=yes' tagged elements, where there is no boundary line from the 'planet_osm_admin' table, but there is a boundary line label from the 'planet_osm_polygon' table.

Maybe that with future enhancements to osm2pgsql or the Lua flex style of carto, it is possible to add 'right' and 'left' name labels to each 'planet_osm_admin' Line record, so as to allow name labeling along the boundary without relying on Polygons. This would solve the potential mismatch between labels and line symbology.

afbeelding