[GIS] Seeking accurate transformation library to convert between EPSG 4326 and 3857 for C#

ccoordinate systemepsgsoftware-recommendations

I've used Proj.Net library to convert X, Y in EPSG 3857 to latitude and longitude coordinates in WGS84 and vice-versa, but the conversion of latitude is not accurate compared to http://spatialreference.org/ref/sr-org/6864/.

For example, if the input is:

5676207.469887, 4291510.258087 then outputs are:

  • 50.990239260333, 36.116815741834 using spatialreference.org and
  • 50.990239260328906 36.116603266233213 using Proj.Net

Below is my code:

private const string Wkt3857 = "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.0174532925199433, AUTHORITY[\"EPSG\", \"9102\"]], " +
                                       "AXIS[\"E\", EAST], AXIS[\"N\", NORTH], " +
                                       "AUTHORITY[\"EPSG\",\"4055\"]], " +
                                       "PROJECTION[\"Mercator\"], " +
                                       "PARAMETER[\"False_Easting\", 0], " +
                                       "PARAMETER[\"False_Northing\", 0], " +
                                       "PARAMETER[\"Central_Meridian\", 0], " +
                                       "PARAMETER[\"Latitude_of_origin\", 0], " +
                                       "UNIT[\"metre\", 1, AUTHORITY[\"EPSG\", \"9001\"]], " +
                                       "AXIS[\"East\", EAST], AXIS[\"North\", NORTH], " +
                                       "AUTHORITY[\"EPSG\",\"3785\"]]";


        public static void XyFromLatLng(double lat, double lng, out double x, out double y)
        {
            var gcsWgs84 = GeographicCoordinateSystem.WGS84;
            var cf = new CoordinateSystemFactory();
            var f = new CoordinateTransformationFactory();
            var sys3857 = cf.CreateFromWkt(Wkt3857);

            var transformTo3857 = f.CreateFromCoordinateSystems(gcsWgs84, sys3857);
            double[] fromPoint = { lng, lat };
            double[] output = transformTo3857.MathTransform.Transform(fromPoint);
            x = output[0];
            y = output[1];
        }

        public static void LatLngFromXY(double x, double y, out double lat, out double lng)
        {
            var gcsWgs84 = GeographicCoordinateSystem.WGS84;
            var cf = new CoordinateSystemFactory();
            var f = new CoordinateTransformationFactory();
            var sys3857 = cf.CreateFromWkt(Wkt3857);

            var transformToWgs84 = f.CreateFromCoordinateSystems(sys3857, gcsWgs84);
            double[] fromPoint = { x, y };
            double[] output = transformToWgs84.MathTransform.Transform(fromPoint);
            lat = output[1];
            lng = output[0];
        }

So I need a library in C# to transform these two CRSs to each other as accurate as possible.

I used DotSpatial for projection and the values are accurate compared to PostGIS. But I don't find out the reason of latitude difference required using Proj.Net.

Best Answer

New GDAL shows different OGC WKT than spatialreference.org which is not extremely well up-to-data nowadays. AUTHORITY["EPSG","3785"] Perhaps you should have a try with this:

gdalsrsinfo epsg:3857

PROJ.4 : '+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'

OGC WKT :
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"]]

However, gdaltransform utility from GDAL version 2.2 that is using the projection definitions which are stored into the "epsg" file gives yet another result:

gdaltransform -s_srs epsg:3857 -t_srs epsg:4326
5676207.469887 4291510.258087
50.9902392603289 35.9337569752097 0

Perhaps you must first make some more evaluation about what is the expected and correct result.

Related Question