ST_AsMVT vector tiles distort at low zoom levels

Question

I'm serving up a data from PostGIS by converting it into a tile and then bytes and then onto a leaflet map. It all works, with one exception, the latitudes stretch ever greater towards the poles when zoomed out, but this effect disappears when zooming in. Both pink and blue lines are from the same data set.

When zoomed in the dataset come together/align suggesting there isn't a problem with the initial dataset.

When zooming out to levels of 5 and below the tiles produced seem to distort and stretch ever more towards the poles. The tiles near the equator are effected a lot less.

This is the code python that is used to produce the tiles from the server:

#

def tile_ul(x, y, z):
    n = 2.0 ** z
    lon_deg = x / n * 360.0 - 180.0
    lat_rad = arctan(sinh(pi * (1 - 2 * y / n)))
    lat_deg = degrees(lat_rad)
    return  lon_deg,lat_deg
############################################################
def get_tile(z,x,y):
    xmin, ymin = tile_ul(x, y, z)
    xmax, ymax = tile_ul(x + 1, y + 1, z)
    query = """SELECT ST_AsMVT(tile) FROM (
               SELECT id, ST_AsMVTGeom(geom, ST_MakeEnvelope( %s, %s, %s, %s ,4326), 4096, 256, false ) geom 
               FROM seaarea ) as tile"""
    cursor = db.connection.cursor()
    cursor.execute(query,(xmin,ymin,xmax,ymax))
    tile = bytes(cursor.fetchone()[0])
    cursor.close()
    return tile

the client code which is used before adding the tile layer to a leaflet map:

var url = "http://localhost:5000/tiles/{z}/{x}/{y}"
var mapillaryLayer = L.vectorGrid.protobuf(url, {vectorTileLayerStyles: {
            fillColor: 'red'}}).addTo(mymap);

I believe the issue may lay in the PostGIS query or possibly the bounds created for the tiles when at low zoom levels. A similar issue is report here

Any suggestions on how to make the tiles not distort at low zoom levels?

Stay Safe All.

Answer 1

The tile_ul function computes lon/lat values (EPSG:4326), assuming that its inputs are tile coordinates in the Pseudo-Mercator projection (EPSG:3857). The code then uses the lon/lat values to clip the geometries that have been read from the database, and returns them directly.

So the data returned by this server uses a wrong projection.

You have to convert the geometries from EPSG:4326 to EPSG:3857 before outputting the tile. And it would be easier to let PostGIS compute the tile bounds; something like this:

WITH webmercator(envelope) AS (
  SELECT ST_TileEnvelope(%s, %s, %s)   -- z, x, y
),
wgs84(envelope) AS (
  SELECT ST_Transform((SELECT envelope FROM webmercator), 4326)
),
geometries(id, geom) AS (
  SELECT id, ST_Transform(geom, 3857)
  FROM seaarea
  WHERE geom && (SELECT envelope FROM wgs84)   -- do not read objects outside the tile
)
SELECT ST_AsMVT(tile) FROM (
  SELECT id, ST_AsMVTGeom(geom, (SELECT envelope FROM webmercator))
  FROM geometries
) AS tile

---

NOTE: improved answers with clipping here and here.