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

asked 12 Apr, 18:14

Terry%20Wilkinson's gravatar image

Terry Wilkinson
616212643
accept rate: 25%

edited 12 Apr, 18:15

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

(13 Apr, 16:52) Vlad

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?

(13 Apr, 19:29) Terry Wilkinson
Replies hidden
1

You need all the digits, or...

(15 Apr, 09:48) Breck Carter

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

(15 Apr, 10:09) Volker Barth

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
permanent link

answered 15 Apr, 16:58

JBSchueler's gravatar image

JBSchueler
2.4k21242
accept rate: 16%

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?

(16 Apr, 10:02) Terry Wilkinson
Replies hidden

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?

(16 Apr, 11:13) Volker Barth

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

(16 Apr, 11:43) Terry Wilkinson

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

(16 Apr, 14:37) JBSchueler
3

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.

(16 Apr, 20:39) Terry Wilkinson

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.

(02 May, 10:43) Terry Wilkinson
Replies hidden

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?

(02 May, 10:55) JBSchueler
showing 4 of 7 show all flat view
Your answer
toggle preview

Follow this question

By Email:

Once you sign in you will be able to subscribe for any updates here

By RSS:

Answers

Answers and Comments

Markdown Basics

  • *italic* or _italic_
  • **bold** or __bold__
  • link:[text](http://url.com/ "title")
  • image?![alt text](/path/img.jpg "title")
  • numbered list: 1. Foo 2. Bar
  • to add a line break simply add two spaces to where you would like the new line to be.
  • basic HTML tags are also supported

Question tags:

×136
×24

question asked: 12 Apr, 18:14

question was seen: 233 times

last updated: 02 May, 10:59