Use PostGIS to get Web Mercator tiles covering a polygon

Question

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.

Answer 1

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.