Conversion from 4326 to 3857: Y coordinate wrong

[skrawn]

I'm trying convert GPS coordinates to the Pseudo-Mercator system (because it works worldwide) and the Y coordinate is always off; the X coordinate matches all of the calculators I've tried.

CoordinateSystemFactory c = new CoordinateSystemFactory();
ICoordinateSystem target = c.CreateFromWkt(@"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.0174532925199433,
     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],
  PARAMETER[""latitude_of_origin"",0],
  UNIT[""metre"",1,
    AUTHORITY[""EPSG"",""9001""]],
  AXIS[""X"", EAST],
  AXIS[""Y"", NORTH],
  AUTHORITY[""EPSG"", ""3857""]]");
            
ICoordinateSystem source = c.CreateFromWkt(@"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.0174532925199433,
      AUTHORITY[""EPSG"", ""9122""]],
    AUTHORITY[""EPSG"", ""4326""]]");

CoordinateTransformationFactory trf = new CoordinateTransformationFactory();
ICoordinateTransformation tr = trf.CreateFromCoordinateSystems(source, target);

double[] pointOrigin = new double[] { -86.070885, 39.521519 };
double[] pointDesintation = new double[] { -88.025821, 41.979628 };

// Gives { X: -9581367.0903264079, Y: 4769456.3827050393 }
double[] originConvert = tr.MathTransform.Transform(pointOrigin);
// Gives { X: -9798989.5703798477, Y: 5129340.707539645 }
double[] destConvert = tr.MathTransform.Transform(pointDesintation);

If I use the calculator on https://epsg.io, the X value is correct but the Y value is off:
https://epsg.io/transform#s_srs=4326&t_srs=3857&x=-86.0708850&y=39.5215190
https://epsg.io/transform#s_srs=4326&t_srs=3857&x=-88.0258210&y=41.9796280

This changes the distance between these two points from the correct value of 318.99km to 420.57km, a pretty substantial difference. Any ideas about why the calculation is off by so much?

[terribletim]

Instead of instantiating 3857 from WKT, try using ProjNet.CoordinateSystems.ProjectedCoordinateSystem.WebMercator;

WKT is an ellipsoidal definition, while WebMercator uses spherical earth. See more here

[skrawn]

Perfect, thank you! I am now getting the correct Y values.

[skrawn]

Actually, despite getting the correct Y values, I still get an incorrect distance when I use Point.Distance. It's coming out to 421.76km now. Did I miss something else?

[skrawn]

Actually, despite getting the correct Y values, I still get an incorrect distance when I use Point.Distance. It's coming out to 421.76km now. Did I miss something else?

Nevermind, this all moot. I can just use the GeoCoordinate and GetDistanceTo functions.

[ca0v]

I understand you've made a distinction that WKT is not intended for 3857 because it assume a spherical earth. And yet epsg.io/3857.wkt does exist as does ProjectedCoordinateSystem.WebMercator.WKT. Why is it that ProjNet cannot read this epsg.io version?

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.0174532925199433,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],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"],AUTHORITY["EPSG","3857"]]
Or maybe the question is, why doesn't epsg.io use this version found in ProjNet:

PROJCS["WGS 84 / Pseudo-Mercator", GEOGCS["WGS 84", DATUM["World Geodetic System 1984", SPHEROID["WGS 84", 6378137, 298.257223563, AUTHORITY["EPSG", "7030"]], AUTHORITY["EPSG", "6326"]], PRIMEM["Greenwich", 0, AUTHORITY["EPSG", "8901"]], UNIT["degree", 0.0174532925199433, AUTHORITY["EPSG", "9102"]], AUTHORITY["EPSG", "4326"]], UNIT["metre", 1, AUTHORITY["EPSG", "9001"]], PROJECTION["Popular Visualisation Pseudo-Mercator", AUTHORITY["EPSG", "3856"]], PARAMETER["latitude_of_origin", 0], PARAMETER["central_meridian", 0], PARAMETER["false_easting", 0], PARAMETER["false_northing", 0], AUTHORITY["EPSG", "3857"], AXIS["East", EAST], AXIS["North", NORTH]]

From https://alastaira.wordpress.com/2011/01/23/the-google-maps-bing-maps-spherical-mercator-projection/ it seems proj.4 solved it by introducing "nadgrids". Maybe you can read the EXTENSION data from the epsg version since it has the "nadgrids"?

[terribletim]

Maybe you can read the EXTENSION data from the epsg version since it has the "nadgrids"?

I think the nadgrids hack forced the semi-minor axis to be read from EXTENSION data. Also seems that only 3857 and deprecated 3785 use extension data, so little benefit implementing EXTENSION data for single coordinate reference system.

@airbreather / @FObermaier is there some way these re-occurring issues with 3857 can be resolved? It seems this is the intention of the CoordinateSystemServices.DefaultInitializer, but doesn't help when reading a WKT for 3857, and also gets over-written when loading multiple WKT definitions (eg SharpMap). Can I suggest a small change in CoordinateSystemWktReaderReadProjectedCoordinateSystem:

Replace:

IProjection projection = new Projection(projectionName, paramList, projectionName, authority, authorityCode, String.Empty, String.Empty, string.Empty);

With:

IProjectedCoordinateSystem projectedCS = null;
if (authorityCode == 3857 && projection.ClassName == "Mercator_1SP")
    projectedCS = ProjectedCoordinateSystem.WebMercator;
else
    projectedCS = new ProjectedCoordinateSystem(geographicCS.HorizontalDatum, geographicCS, unit as LinearUnit, projection, axes, name, authority, authorityCode, String.Empty, String.Empty, String.Empty);

[terribletim]

Duplicate of #37