|
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 |
|
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 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 '19, 10:02)
Terry Wilkinson
Replies hidden
Have you noticed the Remarks section of the ALTER SPATIAL REFERENCE SYSTEM statement?
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 '19, 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 '19, 11:43)
Terry Wilkinson
Looks like that ALTER remark should be included with CREATE OR REPLACE also. I'll make a note.
(16 Apr '19, 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 '19, 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
(02 May '19, 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 '19, 10:55)
JBSchueler
|
Maybe the thing you discovered is related to IEEE 754? It seems that the precision was lost during whatever calculations...
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?
You need all the digits, or...
You say that could mean "serious loss" not limited to precision? :)