Lost polygon after #239

[QuentinC]

Hello,

During my tests on this software, I have seen a lost polygon. It displays fine before #239 but disappears after... This is the lost polygon: https://www.openstreetmap.org/relation/2980911

Well it might be more complicated, here are some screenshots taken before #239 : Zoom 12: Zoom 13: And Zoom 14:

As you can see some parts of the polygon are lost depending of the zoom level, and the whole polygon disappear when using latest commit... I have dont my tests using the France>Bretagne export from there: https://download.geofabrik.de/europe/france/bretagne-latest.osm.pbf with a custom process.lua that includes military landuse: process-openmaptiles.lua.txt (L496)

Thanks again for the support !

[toshelp]

is_valid checks geometry in the sense of OCG. I think this is the result of that check. https://www.boost.org/doc/libs/1_76_0/libs/geometry/doc/html/geometry/reference/algorithms/is_valid/is_valid_2_with_failure_value.html

I found a workaround to uncheck failure_self_intersections. But, it may cause other problems.

diff --git a/include/osm_lua_processing.h b/include/osm_lua_processing.h
index 5911f66..8c1cfa9 100644
--- a/include/osm_lua_processing.h
+++ b/include/osm_lua_processing.h
@@ -125,7 +125,7 @@ public:
         } else if (isWay && !geom::is_valid(geom,failure)) {
             if (verbose) std::cout << "Way " << originalOsmID << " has " << boost_validity_error(failure) << std::endl;
         }
-        if (failure==boost::geometry::failure_self_intersections || failure == boost::geometry::failure_few_points) 
+        if (failure == boost::geometry::failure_few_points) 
                        return false;
                if (failure==boost::geometry::failure_spikes)
                        geom::remove_spikes(geom);

[kleunen]

Indeed what toshelp is saying, i changed the behaviour of the is_valid check to not process self_intersecting polygons. On my map i had an issue that a tile was rendered incorrectly because of a self-intersecting polygon. I am not sure if the OSM data really has this many self-intersecting polygons, or boost is too picky. Because boost also seems to reject polygons which appear to be valid.

The area is a multipolygon, so maybe it would be good to try to "salvage" the geometry. Only keeping the polygons which are valid. I think the pier area has self-intersection, so possibly by removing this, the rest can be rendered.

[kleunen]

Try also run with --verbose, to see if the polygon is mentioned

[kleunen]

Also on openstreetmap.org itself, you see rendering issues with this relation:

[systemed]

Yes, that's definitely a bad multipolygon - at https://www.openstreetmap.org/relation/2980911#map=19/48.36403/-4.51177 and at https://www.openstreetmap.org/relation/2980911#map=19/48.36266/-4.51813 you can see the outers overlapping each other.

[kleunen]

But would a union of these overlapping polygons fix this?

i tried this but it does not seem to fix the polygon.

[QuentinC]

Yes, that's definitely a bad multipolygon - at https://www.openstreetmap.org/relation/2980911#map=19/48.36403/-4.51177 and at https://www.openstreetmap.org/relation/2980911#map=19/48.36266/-4.51813 you can see the outers overlapping each other.

Should I report it to OSM ?

[toshelp]

I tried a union test. The overlapping areas are merged. Is there a rounding problem?

test code

#include <iostream>
#include <vector>
#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <boost/geometry/geometries/box.hpp>
#include <boost/assign/list_of.hpp>

namespace bg = boost::geometry;

typedef bg::model::d2::point_xy<double> point;
typedef bg::model::polygon<point> polygon;
typedef bg::model::box<point> box;

int main()
{
    const box bx(point(5, 0), point(10, 5));

    polygon poly;
    bg::exterior_ring(poly) = boost::assign::list_of<point>
        (0, 0)
        (0, 5)
        (5, 5)
        (5, 0)
        (0, 0)
        ;

    std::vector<polygon> out;
    bg::union_(bx, poly, out);
    
    std::cout << "poly result: " << bg::dsv(out[0]) << std::endl;
}

result

poly result: (((10, 5), (10, 0), (5, 0), (0, 0), (0, 5), (10, 5)))

[dougkpowers]

I'm experiencing a problem since #239 as well as many polygons I've been importing from OSM have suddenly dropped.

For the sake of rendering a map, why do we care if two outer polygons in a multi-polygon overlap? Instead, can we simply verify that each polygon that makes up the multi-polygon is valid?

[kleunen]

Change this: https://github.com/kleunen/tilemaker/blob/master/include/osm_lua_processing.h#L128

to return true, and see if this helps. Many polygons are self-intersecting, but it seems also many do not cause any rendering issue. So it is difficult to detect when a polygon will cause rendering issues.

[kleunen]

I have this way causing issues: https://www.openstreetmap.org/way/261930721#map=19/52.25461/6.20323

If you enable it, the tile looks like this:

[dougkpowers]

@kleunen I had already tried your suggestion and it did work. I found that the multipolygon was "self-intersecting" because two polygons touched each other (but they didn't overlap). That seems very legitimate.

This seems to be behaving like an experimental feature. Perhaps there could be a command line flag to enable/disable it.

[kleunen]

Well, actually. Many polygons contain some self-intersection, also due to the clipping and coordinate conversion that is happening. You can see that if you run with --verbose. But many times, this is not actually an issue. Mapbox can handle some level of self-intersection on the polygons.

I was thinking of maybe to restore the original behaviour. And add a possibility to the config file to block certain geometries. So if certain geometries cause issues on the map, you can blacklist them from inside the config file.

[systemed]

I have this way causing issues: https://www.openstreetmap.org/way/261930721#map=19/52.25461/6.20323

I was able to fix this particular case with:

Polygon fixed;
geom::convex_hull(mp, fixed);
geom::assign(mp, fixed);

which seems to act as a quick-and-dirty 'dissolve'. Ideally we should run this on each individual ring of a multipolygon, I guess.

[kleunen]

I tried also the buffer approach here: https://wandbox.org/permlink/iUbX79o2m3yeNQWe

From this:

it generates this:

[kleunen]

The convex hull gets generated like this:

[kleunen]

This is when i first do a buffer with distance 0.000003, and then a buffer with distance -0.000003

[kleunen]

And with distance 0.00001 and -0.00001:

it is a valid polygon after this https://wandbox.org/permlink/BJo8aUCp4e3RKGlm

[kleunen]

Maybe we should try this buffer approach on: https://www.openstreetmap.org/relation/2980911

#include <boost/format.hpp>
if(originalOsmID == 261930721) {
    for(auto const &p: mp[0].outer())
        std::cout << (boost::format("(%.10f,%.10f),") % p.x() % p.y()).str() << std::endl;
}
if(!CorrectGeometry(mp)) return;