How to convert MVT geometry back to lat,lng (4326) using postgis functions?

Question

Given some arbitrary geometry in the standard ESPG-4326 (latitude & longitude), that was converted to an MVT tile geometry with ST_AsMvtGeom() postgis function (in web mercator 3857), I need to convert it back to the 4326 using postgis functionality.

select st_astext(
   st_asmvtgeom(
      st_transform(ST_GeomFromText('POINT(-73.985130 40.748817)', 4326), 3857),
      st_tileenvelope(8, 75, 96),
      extent := 4096, buffer := 0, clip_geom := true));

The above returns POINT(1591 890). I need to write a similar SQL statement to convert it back to POINT(-73.985130 40.748817) using the inputs of 'POINT(1591 890)', 8, 75, 96 (geometry + tile coordinates). The result would obviously be slightly different from the original.

Answer 1

This seems to do the job. The st_tileenvelope() is a bit of an overkill, because it is only needed to compute the width/height of a tile in meters (i.e. 40075000/(2^z)). Also, the last st_transform() step is not really needed - the resulting geometry is already in a usable format and has SRID attached to it.

CREATE OR REPLACE FUNCTION decode_mvt_geom(geom geometry, z int, x int, y int, extent int) RETURNS geometry AS $$
SELECT st_transform(
  st_scale(
    st_translate(geom, extent * (x - 2 ^ (z - 1)), extent * (y - 2 ^ (z - 1)))
      , (st_xmax(st_tileenvelope(z, x, y)) - st_xmin(st_tileenvelope(z, x, y))) / extent
      , -(st_ymax(st_tileenvelope(z, x, y)) - st_ymin(st_tileenvelope(z, x, y))) / extent
    )
  , 4326)
$$ COST 1 LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE;

A simpler function without the final projection transformation could be this:

CREATE OR REPLACE FUNCTION decode_mvt_geom(geom geometry, z int, x int, y int, extent int) RETURNS geometry AS $$
SELECT st_scale(
  st_translate(geom, extent * (x - 2 ^ (z - 1)), -extent * (y - 2 ^ (z - 1))),
  40075016.6855785 / (2 ^ z) / extent,
  40075016.6855785 / (2 ^ z) / extent
)
$$ COST 1 LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE;