1

The transformation functions for the Web Mercator projection EPSG:3857 for latitude and longitude in radians to projected coordinates in meters are defined as

R = 6378137
x = R * longitude
y = R * ln(tan(pi/4 + latitude/2))

I'm unable to find corresponding functions for EPSG:3395, the so-called World Mercator projection.

Running a WMS that supports both, x values appear to be identical in both projections. There is however a deviation in y values that increases with increasing distance from the equator.

So how would the function for calculating y have to be written?

1

EPSG:3395 apparently uses the elliptical version of the Marcator projection, written as

a = 6378137              // WGS84 semi-major axis
b = 6356752.3142         // WGS84 semi-minor axis
e = sqrt(1 - b^2 / a^2)  // ellipsoid eccentricity

c = pow((1 - e*sin(latitude)) / (1 + e*sin(latitude)), e/2)

y = a * ln(tan(pi/4 + latitude/2) * c)
x = a * longitude;

For reference, see Map Projections - A Working Manual, (7-7), p.44, where also the inverse formula is shown.

1

You have discovered the reason for the following information extracted from an NGA Advisory Notice on "Web Mercator".

"The NGA Geomatics Office has assessed the use of Web Mercator and other non-WGS 84 spatial reference systems may cause geo-location / geo-coordinate errors up to 40,000 meters. This erroneous geospatial positioning information poses an unacceptable risk to global safety of navigation activities, and department of defense, intelligence community, and allied partner systems, missions, and operations that require accurate and precise positioning and navigation information. The NGA Geomatics Office reminds the community to use DoD approved World Geodetic System 1984 (WGS 84) applications for all mission critical activities."

EPSG:3395 is WGS 84 compliant. EPSG:3857 is NOT! You are seeing the reason why.

I didn't realize that either until I was most of the way thru a project and I was finally told that I couldn't use EPSG:3857 for the above reasons. I scrambled to find a replacement and came up with EPSG:3395. It was a relatively painless conversion as all I had to do was change the projection declaration in the code. However, a major issue then came up when I ran into the problem of conversions between EPSG:3395 and the other projection I was using, EPSG:4326 which is in LAT/LON rather than meters. Openlayers only supports conversions between EPSG:4326 and EPSG:3857 out of the box. If you want to use any other projection, such as EPSG:3395, you need to use proj4.js and OpenLayers knows to automatically look into proj4.js for the routines it needs. However, you need to prime the pump, so to speak. You need to declare the projection definition string, an example for EPSG:3395 follows.

proj4.defs("EPSG:3395","+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs");

This one can be found at https://epsg.io/3395. I've run into another problem with this though that I am trying solve. The projections don't work the same. I have documented the problem here and I am looking for any information anyone has as to why.

https://stackoverflow.com/questions/51528830/epsg3395-projection-not-providing-map-wrapping

0

The algorithm is the same. There is just a difference in how to interpret the spheroid. EPSG:3857 (Pseudo Mercator) use globe with the same radius for the semi major and minor axis (see parameter a and b in proj4 EXTENSION tag), EPSG:3395 use a ellipsoid of rotation see inverse flatting value in SPHEROID tag.

For the complete algorithm see the OpenStreetMap Wiki under Elliptical (true) Mercator Projection

EPSG:3857 := PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG", 7030]], TOWGS84[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], AUTHORITY["EPSG", 6326]], PRIMEM["greenwich", 0.0, AUTHORITY["EPSG", 8901]], UNIT["degree", 0.0174532925199433, AUTHORITY["EPSG", 9122]], AUTHORITY["EPSG", 4326]],PROJECTION["MERCATOR_1SP"],PARAMETER["false_easting", 0.0], PARAMETER["false_northing", 0.0], PARAMETER["central_meridian", 0.0], PARAMETER["scale_factor", 1.0], UNIT["metre", 1.0, AUTHORITY["EPSG", 9001]],EXTENSION["proj4", "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"], AXIS["X", EAST], AXIS["Y", NORTH],AUTHORITY["EPSG", 3857]]

EPSG:3395 := PROJCS["WGS 84 / World Mercator",GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG", 7030]], TOWGS84[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], AUTHORITY["EPSG", 6326]], PRIMEM["greenwich", 0.0, AUTHORITY["EPSG", 8901]], UNIT["degree", 0.0174532925199433, AUTHORITY["EPSG", 9122]], AUTHORITY["EPSG", 4326]],PROJECTION["MERCATOR_1SP"],PARAMETER["false_easting", 0.0], PARAMETER["false_northing", 0.0], PARAMETER["central_meridian", 0.0], PARAMETER["scale_factor", 1.0], UNIT["metre", 1.0, AUTHORITY["EPSG", 9001]], AXIS["Easting", EAST], AXIS["Northing", NORTH],AUTHORITY["EPSG", 3395]]

  • The ellipsoidal and spherical algorithms are not the same. Clemens's question has the equations for working with a sphere (note that only the radius / semimajor axis is used). The ellipsoidal Mercator equations are more complicated. – mkennedy Oct 20 '17 at 20:14
  • @Clemens I think I'm seeing an error in the OSM doc, but I'm not 100%. Try using Snyder instead. You might also try changing the projection name for 3395 to just "Mercator" to see if you get different values. – mkennedy Oct 20 '17 at 20:15
  • @mkennedy The OSM methods seem to be wrong where they do ts = tan(0.5 * (M_PI * 0.5 - phi)) / con, which should actually be ts = tan(0.5 * (M_PI * 0.5 - phi)) * con – Clemens Oct 23 '17 at 10:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.