Statistics of MVT tile GIS errors when encoding/decoding with ST_AsMVTGeom

query.sql

SELECT x,
       y,
       ST_Y(p_mid)                                                                          mid_lat,
       ST_Y(p_min)                                                                          min_lat,
       ST_Y(p_max)                                                                          max_lat,
       ST_Y(d_mid)                                                                          mid_lat_decoded,
       ST_Y(d_min)                                                                          min_lat_decoded,
       ST_Y(d_max)                                                                          max_lat_decoded,
       abs(ST_Y(p_mid) - ST_Y(d_mid))                                                       mid_lat_error,
       abs(ST_Y(p_min) - ST_Y(d_min))                                                       min_lat_error,
       abs(ST_Y(p_max) - ST_Y(d_max))                                                       max_lat_error,
       ST_DistanceSpheroid(p_mid, ST_SetSRID(ST_MakePoint(ST_X(p_mid), ST_Y(d_mid)), 4326)) mid_lat_sph_error,
       ST_DistanceSpheroid(p_min, ST_SetSRID(ST_MakePoint(ST_X(p_min), ST_Y(d_min)), 4326)) min_lat_sph_error,
       ST_DistanceSpheroid(p_max, ST_SetSRID(ST_MakePoint(ST_X(p_max), ST_Y(d_max)), 4326)) max_lat_sph_error,
       ST_Y(p_mid_3857)                                                                     mid_y,
       ST_Y(p_min_3857)                                                                     min_y,
       ST_Y(p_max_3857)                                                                     max_y,
       ST_Y(d_mid_3857)                                                                     mid_y_decoded,
       ST_Y(d_min_3857)                                                                     min_y_decoded,
       ST_Y(d_max_3857)                                                                     max_y_decoded,
       abs(ST_Y(p_mid_3857) - ST_Y(d_mid_3857))                                             mid_y_error,
       abs(ST_Y(p_min_3857) - ST_Y(d_min_3857))                                             min_y_error,
       abs(ST_Y(p_max_3857) - ST_Y(d_max_3857))                                             max_y_error,
       abs(ST_Y(p_mid_3857) - ST_Y(d_mid_3857_2))                                           mid_y_error_2,
       abs(ST_Y(p_min_3857) - ST_Y(d_min_3857_2))                                           min_y_error_2,
       abs(ST_Y(p_max_3857) - ST_Y(d_max_3857_2))                                           max_y_error_2
FROM (SELECT *,
             -- and transform back to 4326
             ST_Transform(d_mid_3857, 4326) d_mid,
             ST_Transform(d_min_3857, 4326) d_min,
             ST_Transform(d_max_3857, 4326) d_max
      FROM (SELECT *,
                   -- decode MVT geometries with translate, scale
                   ST_Scale(ST_Translate(v_mid, ext * (x - 2 ^ (z - 1)), ext * (y - 2 ^ (z - 1))),
                            (ST_XMax(env3857) - ST_XMin(env3857)) / ext,
                            -(ST_YMax(env3857) - ST_YMin(env3857)) / ext
                       ) d_mid_3857,
                   ST_Scale(ST_Translate(v_min, ext * (x - 2 ^ (z - 1)), ext * (y - 2 ^ (z - 1))),
                            (ST_XMax(env3857) - ST_XMin(env3857)) / ext,
                            -(ST_YMax(env3857) - ST_YMin(env3857)) / ext
                       ) d_min_3857,
                   ST_Scale(ST_Translate(v_max, ext * (x - 2 ^ (z - 1)), ext * (y - 2 ^ (z - 1))),
                            (ST_XMax(env3857) - ST_XMin(env3857)) / ext,
                            -(ST_YMax(env3857) - ST_YMin(env3857)) / ext
                       ) d_max_3857,
                   ST_Scale(ST_Translate(v_mid, ext * (x - 2 ^ (z - 1)), ext * (y - 2 ^ (z - 1))),
                            40075016.68557849 / (2 ^ z) / ext,
                            -40075016.68557849 / (2 ^ z) / ext
                       ) d_mid_3857_2,
                   ST_Scale(ST_Translate(v_min, ext * (x - 2 ^ (z - 1)), ext * (y - 2 ^ (z - 1))),
                            40075016.68557849 / (2 ^ z) / ext,
                            -40075016.68557849 / (2 ^ z) / ext
                       ) d_min_3857_2,
                   ST_Scale(ST_Translate(v_max, ext * (x - 2 ^ (z - 1)), ext * (y - 2 ^ (z - 1))),
                            40075016.68557849 / (2 ^ z) / ext,
                            -40075016.68557849 / (2 ^ z) / ext
                       ) d_max_3857_2
            FROM (SELECT *,
                         -- encode center, min and max points as MVT geometries
                         ST_AsMVTGeom(p_mid_3857, env3857, ext) v_mid,
                         ST_AsMVTGeom(p_min_3857, env3857, ext) v_min,
                         ST_AsMVTGeom(p_max_3857, env3857, ext) v_max
                  FROM (SELECT *,
                               -- create points at the center of the tile, at a min,min corner, and a max,max corner
                               ST_Transform(p_mid, 3857) p_mid_3857,
                               ST_Transform(p_min, 3857) p_min_3857,
                               ST_Transform(p_max, 3857) p_max_3857
                        FROM (SELECT *,
                                     -- create points at the center of the tile, at a min,min corner, and a max,max corner
                                     ST_Centroid(env4326)                                               p_mid,
                                     ST_SetSRID(ST_MakePoint(ST_XMin(env4326), ST_YMin(env4326)), 4326) p_min,
                                     ST_SetSRID(ST_MakePoint(ST_XMax(env4326), ST_YMax(env4326)), 4326) p_max
                              FROM (SELECT *,
                                           -- create a tile bbox in 3857 and 4326
                                           ST_TileEnvelope(z, x, y)                     env3857,
                                           ST_Transform(ST_TileEnvelope(z, x, y), 4326) env4326
                                    FROM (
                                             -- list of x,y tiles to test
                                             SELECT 14 as z, abs(x - 1) as x, abs(y - 1) as y, 409600 as ext
                                             FROM generate_series(0, (2 ^ 14)::int, (2 ^ 10)::int) as x,
                                                  generate_series(0, (2 ^ 14)::int, (2 ^ 10)::int) as y
                                             -- close all
                                         ) t) t) t) t) t) t) t
order by y, x;