5

I have a polygon in a PostGIS database and I need a list of the Web Mercator tile X/Y pairs that cover the polygon at a particular zoom level. Let's say the polygon is the boundary of the France and the zoom level is 15. For an example of what I mean by X/Y pairs, the slippy map on this site presents the concept visually.

What I've tried so far is using ST_SquareGrid() to create a level-15 grid/graticule for the entire world, and then intersecting it with the boundary of France. The grid took forever to generate and I didn't realize it would be ~1 billion polygons. So far, any operation against that grid has been very expensive and slow, even after indexing. I'm wondering if there's a more straightforward way of doing this that I'm missing.

3
  • 3
    Not a PostGIS solution, so only commenting, but I have found github.com/mapbox/mercantile handy for doing this (Python/CLI). Commented Nov 20, 2022 at 21:01
  • I highly recommend using mercantile.readthedocs.io/en/stable/cli.html#tiles instead of a geometric database operation. Commented Nov 20, 2022 at 21:01
  • This looks like a fantastic library, thank you @alphabetasoup and @bugmenot123! I haven't used Python in a while but I can dust it off and try this out.
    – bertday
    Commented Nov 20, 2022 at 22:01

1 Answer 1

2

While this can be highly optimized algorithmically in a custom low-level implementation, if your Polygons of Interest exist in PostgreSQL/PostGIS already I dare say don't bother - just do it right:

  • ST_SubDivide your POIs; since indexation is useless here, you can run this on-the-fly
  • create a coverage grid per POI and filter for ST_Intersects; since the vertex count per geometry is highly reduced above, intersection checks are much cheaper
WITH
  -- parameter injection, for convenience
  zoom(lvl, csize) AS (
    VALUES ( <ZOOM>, (2*PI()*6378137)/POW(2, <ZOOM>) )
  ),

  -- subdivide your polygons to minimize per-geometry vertex count
  poi AS (
    SELECT
      id, sdv AS geom
    FROM
      <POIs> AS ply,
      LATERAL ST_SubDivide(
        ST_Transform(ply.geom, 3857),
        <MAX_VERTICES>
      ) AS sdv
  )

-- get all covering tile indices for each POI
SELECT
  t.id AS poi_id,
  (grid.i, grid.j, z.lvl) AS index
FROM
  zoom as z,
  poi AS t,
  LATERAL ST_SquareGrid(z.csize, t.geom) AS grid

-- filter for those that actually intersect any of the subdivisions
WHERE
  ST_Intersects(t.geom, grid.geom)

-- return uniques only; much faster than a DISTINCT for multi-column
GROUP BY
  t.id, (grid.i, grid.j, z.lvl)
;

This takes ~6s on a 220k vertex boundary of France for <ZOOM> = 15 & <MAX_VERTICES> = 64.


If you plan to do this over and over, ST_SubDivide (with lower <MAX_VERTICES> value) your POIs into a separate table and run the main query on those.

1
  • 1
    This is incredible, thank you! One thing I'm noticing is that the origin of these tile indices is Web Mercator (0, 0), whereas vector tiles tend to index starting in the upper-left corner and increase rightward and downward, like a bitmap. (Google Maps and Mapbox do this, for instance.) In hindsight I should have been clearer in my question! I'm wondering how I might adapt this for an upper-left origin...
    – bertday
    Commented Nov 22, 2022 at 18:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.