Use PostGIS to get OSM tile coordinates covering a geometry linestring

Question

I have a linestring in SRID 4326 in my table. I would like to know which tiles the linestring crosses regarding a zoom level on OpenStreetMap. The solution found here looks to come close. Use PostGIS to get Web Mercator tiles covering a polygon

But the problem with the given answer, is that the coordinates do not match the z/x/y coordinates that a tileserver is using.

I'm still learning about different projections, it might be that that's where its going wrong?

Let's say I would like to have the coordinates for zoom level 14, This would be the SQL I'm running.

WITH
  -- parameter injection, for convenience
  zoom(lvl, csize) AS (
    VALUES ( 14, (2*PI()*6378137)/POW(2, 14) )
  ),

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

-- 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)
;

But where i'm hoping to get tiles from the netherlands around(z=14 x=8480 y=5319), i'm getting (i=293 j=2875 z=14).

Is using a different SRID and size for st_squaregrid the solution. Working with a size of 0.022 (length of tile at zoomlevel 14 according to https://wiki.openstreetmap.org/wiki/Zoom_levels) and SRID 4326 gives the same results...

Answer 1

In the first post you reference, the same issue seems to happen as well (see the asker's comment to the answer). Here is a solution for you.

Background:

1. The center coordinate of the main Web Map Tiling systems is the top left corner of level 0.

2. ST_SquareGrid() returns tile numbers with an origin centered on the CRS Origin, in your case, with CRS 4326 or 3857, this is located at Lat = 0 & Lon = 0.

3. ST_SquareGrid coordinates will be (0,0) where the Web Map Tile index will be (2^(level-1) , 2^(level-1)). 2^(level-1) being half the vertical/horizontal amount of tiles at the zoom level) which brings you to the center of the overall tiling system.

4. This means that you need compensations of 2^(level-1) somewhere , depending on which Lat/Lon quadrant your geoetries are See interactive tile numbering

In your case, ff all your polygons are East from the prime meridian and in the Northern hemisphere, you should be good with these formulas :

• X = 8192+i =8192+293= 8485
• Y = 8192-j = 8192-2875= 5317

This seem to match with your current results

Visual representation of both tiling systems: