maplibre / martin

Blazing fast and lightweight PostGIS, MBtiles and PMtiles tile server, tile generation, and mbtiles tooling.
https://martin.maplibre.org
Apache License 2.0
2.23k stars 210 forks source link

Tiles along the antimeridian are always empty for table sources with a positive buffer and a SRID of 4326 #982

Open dhwa-kyle opened 11 months ago

dhwa-kyle commented 11 months ago

I'm working with a polygonal table with srid 4326 in PostGIS, whose extents span the non-contiguous United States. I'd like to be able to visualize the entire dataset at once, which means I need to generate tiles at a very low zoom level. I have Martin set up with the simplest table source configuration imaginable, specifically something like the following:

    table_id:
      schema: public
      table: table_name
      geometry_column: geometry
      srid: 4326
      minzoom: 2
      maxzoom: 13
      geometry_type: POLYGON

Martin is impressively fast at generating a zoom level 2 tile for the eastern US (tile (2, 1, 1)), but generates an empty tile for the western US (tile (2, 0, 1)), despite there being quite a bit of data within the extents of the tile.

In looking into it, I've found that seemingly all tiles along the antimeridian (those for which the tile coordinate x = 0 or x = 2^n - 1) end up empty when a positive buffer is specified on a table source with SRID 4326. The crux of the issue is that the search bounding box for such tiles ends up crossing the antimeridian when there is a positive buffer. But SRID 4326 polygons that cross the antimeridian don't operate well with bounding box intersections due to the circular nature of longitudes. So when Martin tries to get geometries within the bounding box, it actually ends up getting features only outside of the actual tile area.

One can find the erroneous code in table_source.rs. Specifically, the buffered bounding box is calculated as something akin to: ST_TileEnvelope(z, x, y, margin => buffer / extent) . The features that go into tile are then queried via something akin to SELECT ... FROM table WHERE geometry && ST_Transform(buffered_bounding_box, srid)

For the case of tile (2, 0, 1) above, the default buffer of 64, and the default extent of 4096, the post-transform buffered bounding box ends up being POLYGON((178.59375 -1.406108835435152,178.59375 67.06743335108297,-88.59375 67.06743335108297,-88.59375 -1.406108835435152,178.59375 -1.406108835435152)), which does not intersect the actual tile in SRID 4326. Indeed, the result of
SELECT ST_Transform(ST_TileEnvelope(2, 0, 1, margin => 0.0), 4326) && ST_Transform(ST_TileEnvelope(2, 0, 1, margin => 64.0 / 4096.0), 4326); is false, which implies that any geometry within the expected tile will not make its way into the resulting tile.

A workaround is to use a buffer of zero, which enforces that the tile envelope doesn't cross the antimeridian. But that is unideal for some purposes, so I am requesting a more permanent solution that properly handles this case.

One possible solution is to use the ST_ShiftLongitude function when working with a tile of zoom greater than one along the antimeridian. Another might be to transform the feature geometry into SRID 3857 (Web Mercator) instead of converting the bounding box to the table's SRID, but I'm unsure of if that would work in all cases.

nyurik commented 11 months ago

@dhwa-kyle thanks for such a detailed bug report, I wish more issues were like this :) I suspect most users do use pre-transformed tables or materialized views in 3857 - as Martin has been doing this for many years and noone reported any issues. That said, we clearly should figure out a good way to solve it for non-3857 case.

I glanced at the pg_tileserv but i don't see any ST_TileEnvelope usage there (might be in a different package?) - just in case they have solved it in some different way. Ironically, I am also the person who implemented the margin param in the ST_TileEnvelope, so I may need to add some comments to the relevant docs to warn people about it.

So I guess we may need to introduce a special case when the detected (or default) SRID=4326, and perhaps try to convert 0-margin to 4326, followed by some growth factor, e.g. ST_Expand(ST_Transform(ST_Envelope(...), 4326), growth_percent)?

With all that, I am not really an expert like @pramsey or @Komzpa - just in case Paul or Darafei know a better solution on how to add ST_TileEnvelope margin for a non-web-mercator data.

sharkAndshark commented 11 months ago

pg_tileserv made it locally, not on pg.
https://github.com/CrunchyData/pg_tileserv/blob/0a999b10adbd8733489d62815c6c2cd533434472/tile.go#L66-L90

nyurik commented 11 months ago

So, what do we want to do about this? I would prefer to keep all geo computations on the PostGIS side for now (although a few people mentioned that it would be awesome to do the geometry->MVT conversion in Rust instead of the database...

sharkAndshark commented 11 months ago

I'm still thinking of this issue.
We can perform the geometry-to-MVT conversion locally, but only when we are dealing with local files such as GeoPackage, Esri GDB, etc.

It might be beneficial to perform the tile bounding box calculation locally, especially if we intend to extend our support to gridsets that aren't 3857.