Serious loss of precision transforming from EPSG:4326 to EPSG:3857 and back

[former_SQLA_member1694873]

I have the following example code:

BEGIN
    DECLARE @lng DOUBLE;
    DECLARE @lat DOUBLE;
    SET @lng = -79.0;
    SET @lat = 43.0;
    SELECT
        NEW ST_Point(@lng,@lat,4326),
        NEW ST_Point(@lng,@lat,4326).ST_Transform(3857).ST_Transform(4326),
        NEW ST_Point(@lng,@lat,4326).ST_Transform(2163).ST_Transform(4326),
        Value
    FROM sa_eng_properties() WHERE PropName = 'MessageText' ORDER BY PropName;
END

The output columns contain:

Point (-79 43)
Point (-78.999999999779291 42.999701923506009)
Point (-78.999999999989356 42.99999999999968)
SQL Anywhere Network Server Version 17.0.9.4935

In the second output line you can see that after converting from EPSG:4326 (round earth) to EPSG:3857 (projected) and back, the longitude is accurate to around the 10th decimal digit, but the latitude is out in the 4th digit - quite a lot when a degree of latitude is around 111,000 metres.

For comparison, I included a similar conversion from EPSG:4326 to EPSG:2163 (also projected) and back, and both the longitude and the latitude are good to around the 10th digit.

Is there something I'm missing? Is there something wrong with the built-in PROJECTION definition for 3857?

FYI, I'm running SA 17.0.9.4935 and see the same behaviour on both Windows and Linux platforms.

Terry

[Vlad]

Maybe the thing you discovered is related to IEEE 754? It seems that the precision was lost during whatever calculations...

[former_SQLA_member1694873]

That is what I thought at first, but why then are the other transformations so accurate? I didn't show it in my example, but the the EPSG:3857 value for point in question is Point (-8794239.7726 5282806.7292) - note that the x and y values each have 7-digits before the decimal and 4-digits after. Unless I'm not understanding IEEE 754 correctly?

[Breck_Carter]

You need all the digits, or...

[VolkerBarth]

You say that could mean "serious loss" not limited to precision? 🙂

[jack_schueler · Accepted solution]

My spatial colleagues here at SAP Waterloo suggest updating the SRS definition as follows (since the supplied version in st_geometry_config.tgz is out-of-date):

create or replace spatial reference system [WGS 84 / Pseudo-Mercator]
identified by 3857
organization 'EPSG' identified by 3857
definition 'PROJCS["Popular Visualisation CRS / Mercator",GEOGCS["Popular Visualisation CRS",DATUM["Popular_Visualisation_Datum",SPHEROID["Popular Visualisation Sphere",6378137,0,AUTHORITY["EPSG","7059"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6055"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4055"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],AUTHORITY["EPSG","3785"],AXIS["X",EAST],AXIS["Y",NORTH]]'
transform definition '+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +nadgrids=@null +wktext +no_defs'
type planar
coordinate x between -20037508.342789 and 20037508.3427892
coordinate y between -19929191.766855 and 19929191.7668548

[former_SQLA_member1694873 · Accepted solution]

Thanks Jack. That fixes the issue. However, I had to create a new database to test it because when I tried it on my existing database I got the message "Could not execute statement. SRID 3857 is referenced by column 'MY_COLUMN' of table 'MY_TABLE' SQLCODE=-1474". It seems I can't redefine the spatial reference system on an existing database. Is there a work around for that? Am I stuck with doing a unload/reload?

[VolkerBarth · Accepted solution]

Have you noticed the Remarks section of the ALTER SPATIAL REFERENCE SYSTEM statement?

Remarks

You cannot alter a spatial reference system if there is existing data that references it. For example, if you have a column declared as ST_Point(SRID=8743), you cannot alter the spatial reference system with SRID 8743. This is because many spatial reference system attributes, such as storage format, impact the storage format of the data. If you have data that references the SRID, create a new spatial reference system and transform the data to the new SRID.

---

I don't know whether it would be fitting to transform those spatial data based on SRID 3857 to another SRID, then alter SRID 3857 and transform it back to that?

[former_SQLA_member1694873 · Accepted solution]

Thanks for the reference Volker - I had not looked at ALTER, only CREATE. I'll try doing the transformations (with appropriate backups of course :-))

[jack_schueler · Accepted solution]

Looks like that ALTER remark should be included with CREATE OR REPLACE also. I'll make a note.

[former_SQLA_member1694873 · Accepted solution]

In case anyone wants to know, I was able to change the spatial reference definition for EPSG:3857 in my database by first changing the SRID for the column to another projection with ALTER TABLE MYTABLE ALTER MY_COLUMN ST_Geometry(SRID=2163), then executing the statement from JBSchueler above, and then using another ALTER TABLE stsatement to set the SRID back to 3857. I need to do a bit more testing but it seems to have worked just fine.

[former_SQLA_member1694873 · Accepted solution]

It turns out that I still encountered precision problems using the spatial reference definition discussed above (from SAP), so after a whole lot of trial and error, I ended up using the following:

create or replace spatial reference system [WGS 84 / Pseudo-Mercator]
identified by 3857
organization 'EPSG' identified by 3857
definition 'PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","3857"]]'
transform definition '+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs '
type planar
coordinate x between -20037508.342789 and 20037508.3427892
coordinate y between -19929191.766855 and 19929191.7668548;

The transform definition clauses are different but I'm not sure what the differences mean. I'm no expert at this so if anyone could offer more explanation, I would really appreciate it.

[jack_schueler · Accepted solution]

I am no expert either but I did learn that the relationship between +a=6378137 +b=6378137 and coordinate x is pi (pi*6378137=20037508.3427892430765884088807‬). But I have no idea what it is between coordinate y. Given that there are no other interesting numbers in the transform definition, I can't imagine why your change would make a difference. We need a "projections" expert.

Also, it is interesting to note that the between x and y uses 6 digits of accuracy for x and 7 digits of accuracy for y. Why?