Inaccurate Transformations - OSGB 1936 / British National Grid to WGS84 with Proj.Net

Question

I am new to all this GIS work, so if I have done something silly or just plain wrong please just let me know.

I m trying to covert from coordinates from OSGB 1936 / British National Grid to WGS84 but I am seeing different results to other tools.

My example coordinates X,Y is; 362895, 155602

My conversions appear to be nearly correct however they are always a little out. Here is a link showing the conversion on epsg.io. If you switch the result to decimal you will see the expected results is;

Long=-2.5335813°, Lat = 51.2983258°

https://epsg.io/transform#s_srs=27700&t_srs=4326&x=362895.0000000&y=155602.0000000

I have also done the same conversion on https://mygeodata.cloud/cs2cs/ with the same results as epsg.io.

However in Proj.Net my result is;

-2.53225739401083, 51.2985133921956

Here is a code example to show my issues;

class Program
    {

        public const string wktOsgb1936 = "PROJCS[\"OSGB 1936 / British National Grid\"," +
        "GEOGCS[\"OSGB 1936\"," +
        "DATUM[\"OSGB_1936\"," +
        "SPHEROID[\"Airy 1830\",6377563.396,299.3249646,AUTHORITY[\"EPSG\",\"7001\"]]," +
        "AUTHORITY[\"EPSG\",\"6277\"]]," +
        "PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]]," +
        "UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]]," +
        "AUTHORITY[\"EPSG\",\"4277\"]]," +
        "UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]," +
        "PROJECTION[\"Transverse_Mercator\"]," +
        "PARAMETER[\"latitude_of_origin\",49]," +
        "PARAMETER[\"central_meridian\",-2]," +
        "PARAMETER[\"scale_factor\",0.9996012717]," +
        "PARAMETER[\"false_easting\",400000]," +
        "PARAMETER[\"false_northing\",-100000]," +
        "AUTHORITY[\"EPSG\",\"27700\"]," +
        "AXIS[\"Easting\",EAST]," +
        "AXIS[\"Northing\",NORTH]]";

        static void Main(string[] args)
        {
            // The orginal 27700 OSG co-ords
            var x = 362895;
            var y = 155602;

            var res = TransformSpatial(wktOsgb1936, x, y);

            // the long/lat which I am expecting as confirmed on other systems
            var shouldx = -2.5335813;
            var shouldy = 51.2983258;

            // A small calc to see the difference
            var difx = res[0] - shouldx;
            var dify = res[1] - shouldy;

            Console.WriteLine("From    : {0}, {1}", x, y);
            Console.WriteLine("conv    : {0}, {1}", res[0], res[1]);
            Console.WriteLine("Should  : {0}, {1}", shouldx, shouldy);
            Console.WriteLine("Diff    : {0}, {1}", difx, dify);

            Console.ReadLine();
        }

        public static double[] TransformSpatial(string wktFrom, double x, double y)
        {
            var cf = new CoordinateSystemFactory();
            var f = new CoordinateTransformationFactory();
            var sysFrom = cf.CreateFromWkt(wktFrom);
            var transformTo = f.CreateFromCoordinateSystems(sysFrom, ProjNet.CoordinateSystems.GeographicCoordinateSystem.WGS84);
            var ret = transformTo.MathTransform.Transform( x, y );
            return new double[] {
                ret.x, ret.y
            };
        }

    }

Please help and let me know what I am doing wrong, this is doing my head in!!

Is there a problem in my "From" mappings, my "To" mapping, is there another step I am missing or is there a bug?

Thanks in advance

Answer 1

PROJ6/7 (via cs2cs) returned the result as below.

If we convert all those results to sexagesimal, we can see that they all concur except for Proj.Net.

To convert a pair of OSGB 1936 planimetric coordinates to "WGS84" geodetic coordinates, the steps are as follow :

Step 1:
Reverse project the OSGB 1936 easting and northing coordinates.
Output is a pair of geodetic coordinates on the Airy ellipsoid.

Step 2:
Convert the geodetic coordinates to 3D cartesian coordinates.
Output is a set of (geocentric ECEF) XYZ coordinates (in the OSGB 36 reference frame ???).

Step 3:
Apply 7-parameter Helmert transformation.
Output is a set of XYZ coordinates in the "WGS84" reference frame.
The parameter values are as below.

Step 4:
Convert the transformed 3D cartesian coordinates to geodetic coordinates on the WGS84 ellipsoid.

PROJ automatically does all the steps above if the its coordinate system database (proj.db) is present, or if you specify it explicitly e.g., via the "+towgs84" parameter.

It appears that (at this time of writing) Proj.Net does not (yet) do auto coordinate transformation although all the codes are already there, which means Coders with some knowledge in coordinate system will be able to make them work. For now, it stops at Step 2. Thus Proj.Net's output is actually a pair of geodetic coordinates on the Airy ellipsoid (in the OSGB 36 reference frame???)

Proj.Net is a fantastic project for the .NET community, and hopefully we will have more contributors to bring it up to speed.

A caveat:
While PROJ, Mygeodata Cloud, and epsg.io may concur on the output, and be methodologically and mathematically correct, but this doesn't necessarily mean they are authoritatively correct unless someone can validate that the Helmert seven-parameter (shown above) are correct. (I haven't done any UK Ordnance projects, so I can't tell.)