4

Let's say I run the following query on a PostgreSQL database with the PostGIS extension:

SELECT ST_AsText((ST_Transform('SRID=4326;POINT(-34.91148366443491 -8.077603972529051)',3857)));

The result is:

POINT(-3886328.584362583 -902188.317568294)

However, when I try to do it the other way around:

SELECT ST_AsText(ST_Transform('SRID=3857;POINT(-3886328.584362583 -902188.317568294)', 4326 ));

The result is:

POINT(-34.91148366443491 -8.077603972529053)

Even though it was the same conversion the last result is different from the original one. That brings me some problems if I try to compare the last result with the first. How should I be handling this? I thought about diminishing the precision of the final and original results in order to get comparisons right, but I don't think that's the right way of handling it. How can I deal with precision problems generated by the ST_Transform function?

raylight
  • 1,241
  • 4
  • 19
  • 2
    see this answer https://gis.stackexchange.com/questions/8650/measuring-accuracy-of-latitude-and-longitude/8674#8674 for an explanation of why this doesn't matter - or this more succinct xkcd https://xkcd.com/2170/ – Ian Turton Apr 22 '22 at 19:31
  • @IanTurton I see... but still, if I try to compare the first and last result inside a WHERE statement I'll get no results. It's not clear to me how I'd handle it. if I consider that 5 decimal places are enough for me (1 meter precision), what would be the most straighforward way of handling it? Would remaking my original geometry column with 5 decimal places be the most appropriate way of handling it? Or should I just use a function that diminishes its precision when I want to make a comparison with =? – raylight Apr 22 '22 at 20:11
  • 2
    there is the https://postgis.net/docs/ST_ReducePrecision.html 'Returns a valid geometry with all points rounded to the provided grid tolerance.' – Mapperz Apr 22 '22 at 20:16
  • 1
    Also the conversion from double to text and back is lossy. Compare with SELECT ST_AsText( ST_Transform( (ST_Transform( ST_GeomFromText('POINT(-3886328.584362583 -902188.317568294)', 3857 ),4326)),3857)); I recommend to use ST_SnapToGrid https://postgis.net/docs/ST_SnapToGrid.html for dealing with tolerances. – user30184 Apr 22 '22 at 20:22
  • @Mapperz I've tried using this function before but I get the error Precision reduction requires GEOS-3.9 or higher... I'm still trying to figure out a way of installing it on Ubuntu 20.04. – raylight Apr 22 '22 at 20:29

1 Answers1

6

To summarize and expand on the comments:

  • coordinate transformation is inherently imprecise
  • testing for equality between original and reprojected data must be done by limiting precision (down to say 5 decimal places)
  • For Points and LineStrings, ST_SnapToGrid can be used and is fast
  • For polygonal geometry, ST_ReducePrecision should be used, since it preserves geometry validity

The example case:

WITH data(geom) AS (
  VALUES ('SRID=4326;POINT(-34.91148366443491 -8.077603972529051)'::geometry)
),
reproj AS (SELECT geom AS original,
                  ST_Transform( ST_Transform( geom,3857), 4326) reprojected
          FROM data
)
SELECT    ST_SnapToGrid( original,    5)
        = ST_SnapToGrid( reprojected, 5) AS is_equal,
        original, reprojected
FROM  reproj;

is_equal |                   original                   |                 reprojected

----------+----------------------------------------------+---------------------------------------------
 t        | POINT(-34.91148366443491 -8.077603972529051) | POINT(-34.91148366443491 -8.07760397252905)
dr_jts
  • 5,038
  • 18
  • 15