NetTopologySuite / ProjNet4GeoAPI

.NET Spatial Reference and Projection Engine
GNU Lesser General Public License v2.1
270 stars 81 forks source link

Inaccurate Transformations - OSGB 1936 / British National Grid to WGS84 #67

Open lakesol opened 4 years ago

lakesol commented 4 years ago

Hi all, 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)
        {
            var x = 362895;
            var y = 155602;

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

            var shouldx = -2.5335813;
            var shouldy = 51.2983258;

            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

Metritutus commented 4 years ago

I experienced something similar very recently when experimenting with this, but was able to resolve it. The key part appeared to be the WKT string for OSGB 1936.

Your WKT string doesn't quite match the one here (from OGS WKT): https://epsg.io/27700

To be specific, it looks like the part you appear to be missing is TOWGS84[446.448,-125.157,542.06,0.15,0.247,0.842,-20.489].

This site explains the TOWGS84 parameter (and the rest), although I confess I've not yet found the need to read into these: https://www.bgs.ac.uk/data/webservices/convertForm.cfm.

I hope this helps.

pchilds commented 3 years ago

The TOWGS84 string is a hack to provide conversions between datums. It is left off most Wkt lists because it forces modern GIS libraries to use an inaccurate method of transformation in going via WGS84 and ECEF coordinates. ProjNET; however, does need this hack and these strings will have to be added in until it updates to be based on a more recent version of Proj (https://github.com/NetTopologySuite/ProjNet4GeoAPI/issues/89)