Invalid polygons after clipping at tile boundaries

[Nakaner]

While using a fresh new set of Shortbread vector tiles to render a client's style with our Mapnik-GDAL vector-to-raster toolchain, I found huge water polygons across the size of tile on zoom level 6 to 8.

Example (map style rendering only inland water bodies, glaciers and sea):

The affected OSM ways have in common that they cross tile boundaries and, are long and narrow polygons (rivers/canals). The artefacts occur when their narrow "tube" becomes a row of multiple small polygons due to simplification/snapping onto the grid of the vector tiles' internal coordinate system. Examples:

• https://www.openstreetmap.org/way/75687836
• https://www.openstreetmap.org/way/854303758

You can use the following example to reproduce the issue:

{
    "layers": {
        "water_polygons": {
            "minzoom": 4,
            "maxzoom": 14
        }
    },
    "settings": {
        "minzoom": 0,
        "maxzoom": 14,
        "basezoom": 14,
        "include_ids": false,
        "combine_below": 14,
        "name": "Water polygon test",
        "version": "1.0",
        "description": "valid geometries desired",
        "compress": "gzip",
        "filemetadata": {
            "tilejson": "2.0.0",
            "scheme": "xyz",
            "type": "baselayer",
            "format": "pbf",
            "tiles": [
                "https://example.com/liechtenstein/{z}/{x}/{y}.pbf"
            ]
        },
        "metadata": {
          "author": "OpenStreetMap contributors, Geofabrik GmbH",
          "license": "Open Database License 1.0"
        }
    }
}

function zmin_for_area(min_square_pixels, way_area)
    -- Return minimum zoom level where the area of the way/multipolygon is larger than
    -- the provided threshold.
    local circumfence = 40052725.78
    local zmin = (math.log((min_square_pixels * circumfence^2) / (2^16 * way_area))) / (2 * math.log(2))
    return math.floor(zmin)
end

function process_water_polygons(way_area)
    local natural = Find("natural")
    if natural ~= "water" then
        return
    end
    mz = math.max(4, zmin_for_area(0.1, way_area))
    Layer("water_polygons", true)
    MinZoom(mz)
    Attribute("osm_id", Id())
end

function way_function()
    local area = Area()
    if area > 0 then
        process_water_polygons(area)
    end
end

function init_function()
end
function exit_function()
end

node_keys = { }
function node_function()
end

Filter the planet dump:

osmium tags-filter -o planet-water-polygons.osm.pbf -v current-planet.osm.pbf a/natural=water a/waterway a/landuse=reservoir a/landuse=basin a/natural=glacier

Clip the thematic extract:

osmium extract -b 5.32,52.29,5.92,52.68 --strategy smart -o flevoland-water.osm.pbf planet-water-polygons.osm.pbf

Run Tilemaker: tilemaker --input flevoland-water.osm.pbf --output polygon-test.mbtiles --bbox 5.32,52.29,5.92,52.68 --config config-polygon-test.json --process process-polygon-test.lua (needs about 700 MB RAM and took a few seconds)

If you open the resulting vector tile set with QGIS on zoom level 8, you will see a broken polygon like this:

MapLibre GL JS and OpenLayers render the polygon as expected. tippecanoe-decode crashes with Polygon begins with an inner ring.

ogr2ogr -f GeoJSON poly.json polygon-test.mbtiles -oo ZOOM_LEVEL=8 -where "osm_id='1187867289'"'" polygon-test.mbtiles reports following GeoJSON:

{
"type": "FeatureCollection",
"name": "water_polygons",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::3857" } },
"features": [
{ "type": "Feature", "properties": { "osm_id": "1187867289" }, "geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ 624949.143259599339217, 6775187.094627310521901 ], [ 624987.361773742013611, 6775187.094627310521901 ], [ 625025.58028788457159, 6775148.876113167032599 ], [ 624949.143259599339217, 6775187.094627310521901 ] ] ], [ [ [ 626172.135712162242271, 6731350.458905761130154 ], [ 624753.464467189158313, 6731350.458905761130154 ], [ 624753.464467189158313, 6887893.492833802476525 ], [ 626172.135712162242271, 6887893.492833802476525 ], [ 626172.135712162242271, 6774843.12800002656877 ], [ 626172.135712162242271, 6731350.458905761130154 ] ], [ [ 625981.043141449335963, 6774957.783542454242706 ], [ 625789.950570736313239, 6775072.439084881916642 ], [ 625522.420971738174558, 6775110.657599024474621 ], [ 625216.672858597477898, 6775148.876113167032599 ], [ 625140.235830312361941, 6775148.876113167032599 ], [ 625293.109886882710271, 6775110.657599024474621 ], [ 625369.546915167826228, 6775110.657599024474621 ], [ 625484.202457595616579, 6775072.439084881916642 ], [ 625484.202457595616579, 6775110.657599024474621 ], [ 625751.73205659375526, 6775072.439084881916642 ], [ 625904.60611316410359, 6774996.002056596800685 ], [ 625981.043141449335963, 6774919.565028311684728 ], [ 626095.698683877009898, 6774881.346514169126749 ], [ 626172.135712162242271, 6774843.12800002656877 ], [ 625981.043141449335963, 6774957.783542454242706 ] ] ] ] } },
{ "type": "Feature", "properties": { "osm_id": "1187867289" }, "geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ 626286.791254593175836, 6774804.909485884010792 ], [ 626325.009768735733815, 6774804.909485884010792 ], [ 626363.228282878291793, 6774843.12800002656877 ], [ 626401.446797020966187, 6774843.12800002656877 ], [ 626363.228282878291793, 6774804.909485884010792 ], [ 626286.791254593175836, 6774766.690971741452813 ], [ 626248.572740450617857, 6774804.909485884010792 ], [ 626172.135712165385485, 6774843.12800002656877 ], [ 626286.791254593175836, 6774804.909485884010792 ] ] ] ] } }
]
}

What does not seem to have any influence:

• building Tilemaker with LuaJIT or plain Lua
• Boost 1.71 vs 1.74
• simplification enabled/disabled in Tilemaker configuration
• Tilemaker 2.4.0 vs. 3.0.0 vs. master (as of yesterday)

Using different OSM datasets containing the same features will make different polygons to break. You can reproduce this by skipping the osmium extract step and running Tilemaker with planet-water-polygons.osm.pbf instead of the Flevoland .osm.pbf file. While converting GDAL reports a couple of self-intersections and ring self-intersections.

[geoneutrino]

seems to be the related/same effect as in https://github.com/systemed/tilemaker/issues/580 Long story shorts, it is difficult ... Pull request 611 solved a lot of the edge cases but this already went into the 3.0 master

although i dont have this effect at this zoom level

[koufopoulosf]

Hello all!

@geoneutrino It's weird that you had this issue #580 with Athens on November/December 2023, as half year earlier it was working for me #495.

@Nakaner If I remember correctly, the discrepancies were because I was using the release of v2.4.0, and not the master branch of github repo directly, which was more updated. @geoneutrino maybe you were using the v2.4.0 release as well?

I haven't checked the v3.0.0 as it includes breaking changes and I need to refactor my stuff now 😛.

In the meantime, let me include this case, as it breaks my heart seeing the map broken.

In Malta (lat: 35.903860, long: 14.503326) at around zoom level 12-13, there's the following issue:

@systemed Regarding the shortbread schema, I remember there were some instructions here: make CONFIG=-DFLOAT_Z_ORDER. [Edit: Found it here]. Is it still the case? It's not required anymore?

Thanks in advance!

[systemed]

Essentially the only sane way to render river/canal polygons at low zoom levels is to skeletonise them, and that's not something tilemaker can yet do[^1]. I usually drop in Natural Earth shapefiles at these zoom levels which are perfect for this use case.

That aside, a significant part of the coastline fix in #611 was to add area filtering to the coastline layer config. Previously, small inners were intersecting with the outers when rescaled to output coordinates, and MBGL really doesn't like self-intersections. Config file area filtering works on the level of individual inners/outers, vs the Lua code you have here which looks at the total area of the multipolygon.

So my first suggestion here would be to experiment with area filtering, e.g.

"filter_below": 12, "filter_area": 0.5

@koufopoulosf -DFLOAT_Z_ORDER is no longer required - tilemaker 3.0 can cope with much wider Z-order values out of the box.

[^1]: I'd love to support it, but boost::geometry doesn't have a skeletonise operation. I think there's a boost::polygon implementation but porting that is beyond my geometry-fu. But more specialised simplify operations like this are definitely an interesting area for future development.

[koufopoulosf]

@systemed Thank you very much for your detailed response.

Is there any chance that there is a documentation/guide on how to export the whole globe correctly for production deployments? Also do you suggest export to mbtiles or in pmtiles format for the whole globe?

Thanks in advance!

[systemed]

The intention is that, if you use the latest code with the default OMT-compatible profile/config, it will generate the planet correctly. If you find significant departures from this then please do report them here.

If you're using your own profile/config then you will of course have to experiment until you get a pleasing and correct result. 🙂

mbtiles vs pmtiles is really a question of how you want to deploy the tiles, and that's outside tilemaker's scope. Personally I use mbtiles because I have a custom server built around that, but pmtiles is an amazing stack and allows you to go "serverless".

[koufopoulosf]

@systemed That's indeed interesting, but going serverless would make sense for personal projects to me, rather than bigger projects as I am not sure if that would be cost effective at a bigger scale😄. Thanks once again for your support.

[geoneutrino]

@koufopoulosf i started last year summer with master branch which was basically v2.4 + some commits not affecting clipping stuff as far as i see in comments. I also use the natural earth shapefile as well as an own config/lua, but still with v3 have some effects Mainly visible in denmark at z5 #5.87/55.147/11.37 And now your Malta :) The small island you mentioned works for me but the southern part of Malta is completely flooded at #9.89/35.8069/14.5662

[koufopoulosf]

@geoneutrino oh no, that's not good news, I was hoping v3 to have it fixed actually😅 (whole island, in general).

So at the end of the day, we agree that we should generate tiles country by country rather than whole globe or whole europe? So that we fix each country individually by every new version of tilemaker?😄

[koufopoulosf]

Hi everyone! 👋

@systemed, could we please get an update on these issues?

I'm currently working on a project where I'm considering setting a maximum zoom level for the map to avoid encountering certain issues. It seems that some of these issues were resolved in previous versions, which makes their persistence a bit perplexing.

If these issues persist and can't be resolved, I'd like to know at which zoom level it would be advisable to stop before encountering them. This way, I can ensure a smoother user experience without running into these specific problems.

Thank you in advance for your assistance! ✌️

[systemed]

Could you confirm whether this happens with the default OMT-compatible schema?

[koufopoulosf]

@systemed I'm using the shortbread schema

[koufopoulosf]

Update: After conducting further testing, it appears that the issue of "missing" land over certain zoom levels may be attributed to how the Lua processing script handles water polygons (shortbread schema related lua). Currently, these polygons overlap significant portions of the land. So probably, it seems that the issue is not directly related to Tilemaker itself.

@geoneutrino, have you made any progress in resolving this issue?

Update 2: @systemed Upon retesting Greece, I've confirmed that the same issue (#580) persists (coordinates: /#8.7/37.7811/24.4236). Previously, I had resolved this by using the -DFLOAT_Z_ORDER option, which you mentioned is no longer necessary. I've tested both the release source and the main branch on an Ubuntu 22.04. Could you please verify whether this issue persists?

Thanks in advance for your assistance!

[systemed]

tl;dr - the default profile is optimised to minimise coastline issues. Please use the same config for the ocean layer and report if you still have significant issues.

---

Applying major simplification to very complex shapes is always likely to produce invalid geometries. Taking coastlines designed for z14+ and reducing them to small scales is an example of that.

tilemaker jumps through a lot of hoops to minimise invalid geometries but it isn't infallible, and you can improve your chances of success by configuring the layer properly. In particular the default has this:

"filter_below": 12, "filter_area": 0.5

which dynamically filters out small geometries, e.g. little islands which would shrink down to a few pixels at lower zoom levels. By filtering out these islands you are greatly reducing the chance of self-intersections.

For my own purposes I use a different ocean dataset for low zooms.

If there is a bug in the default config/profile then please do report it here but I can't undertake to debug every single config/profile created for tilemaker, I'm afraid. Issues specific to shortbread should probably go to the shortbread repo.

[koufopoulosf]

@systemed Thanks for the update!

In terms of the issue with the islands and the config, it's totally understandable. In fact, I reported it here.

The second issue I raised on update 2, was this one. You are suggesting that this issue does not exist in tilemaker currently?

[systemed]

I'm not suggesting anything other than

If there is a bug in the default config/profile then please do report it here

I'm afraid I'm finding this issue very hard to track because of all the different locations being thrown about, and confusion about what's an issue with shortbread and what's an issue with tilemaker and its default config.

It's very possible there's an issue around Athens - I remember encountering problems there before, and indeed #602 was partly tested there. If you can still reproduce a problem using tilemaker's default config, then please do open a Github issue citing that location. Please don't try and generalise it or guess the solution, just report exactly what and where. Thank you!