locked
Converting OSGB36 to WGS84 with ProjNet; Doing My Head In RRS feed

  • Question

  • I realise this isn't really related to SQL Server spatial objects, but I do want to put the finished CLR assembly into SS, so I'm hoping someone will take pity on me.

    I started here:

    http://www.ordnancesurvey.co.uk/oswebsite/gps/osnetfreeservices/furtherinfo/spreadsheet.html

    which has a spreadsheet that converts between OS easting/northing and WGS84 long/lat. The coordinates I used were

        easting = 651409.903, northing = 313177.271

    These are equivalent to (according to the spreadsheet)

        latitude = 1.71605201472, longitude 52.65797861194

    Try as I might, I cannot get the conversion right using ProjNet. I have tried all manor of WKT and built-in coordinate systems, projections and so on, but the best I can get is a point 134m away from where it should be. I did manage to get the right answer using something called GeoCoordConversion.dll, but that lacks support for VE/Bing X/Y (which I also need). And besides, I ought to be able to get it right with ProjNet.

    One attempt used this:

                ICoordinateSystemFactory cFac = new CoordinateSystemFactory();

                IEllipsoid ellipsoid = cFac.CreateFlattenedSphere("Airy 1830", 6377563.396, 299.32496, LinearUnit.Metre);

                IHorizontalDatum datum = cFac.CreateHorizontalDatum("Airy 1830", DatumType.HD_Geocentric, ellipsoid, null);

                IGeographicCoordinateSystem gcs;
                ICoordinateSystem coordsys;

                // Create the transform and projection for use in converting between WGS84 polar and OSGB36 grid coordinates
                gcs = cFac.CreateGeographicCoordinateSystem("Airy 1830", AngularUnit.Degrees, datum, PrimeMeridian.Greenwich,
                    new AxisInfo("Lon", AxisOrientationEnum.East),
                    new AxisInfo("Lat", AxisOrientationEnum.North));

                List<ProjectionParameter> parameters = new List<ProjectionParameter>(5);

                parameters.Add(new ProjectionParameter("latitude_of_origin", 49));
                parameters.Add(new ProjectionParameter("central_meridian", -2));
                parameters.Add(new ProjectionParameter("scale_factor", 0.9996012717));
                parameters.Add(new ProjectionParameter("false_easting", 400000));
                parameters.Add(new ProjectionParameter("false_northing", -100000));

                IProjection projection = cFac.CreateProjection("Transverse Mercator", "Transverse_Mercator", parameters);

                coordsys = cFac.CreateProjectedCoordinateSystem("OSGB 1936 / British National Grid", gcs, projection, LinearUnit.Metre, new AxisInfo("East", AxisOrientationEnum.East), new AxisInfo("North", AxisOrientationEnum.North));

                ICoordinateTransformation = wgs84PolarToOsgb36GridTransform = new CoordinateTransformationFactory().CreateFromCoordinateSystems(gcs, coordsys);

    followed by

                double[] d = _wgs84PolarToOsgb36GridTransform.MathTransform.Transform(new double[] { longitude, latitude });

    but this was 134m out. I think I can see what is wrong with this attempt:it doesn't use WGS84 anywhere, but I have tried all manner of WKT for WGS84 as well and still can't get it.

    Could someone please tell me how to do this simply before I go insane?

    Charles

     

    Friday, September 16, 2011 5:25 PM

Answers

  • We're obviously misunderstanding each other somewhere.... so let's tackle this from a different approach.

    Why not try defining the source and target coordinate systems from which you're creating your coordinatetransformation factory using the CoordinateSystemWktReader.Parse() method instead? That way, you just have to pass in the Well-Known Text of each coordinate system rather than manually setting the projection parameters as you are doing currently.

    There's an example of how to define the WGS84 geographic coordinate system from WKT in the Proj.NET documentation, at the top of the following page: http://projnet.codeplex.com/wikipage?title=CreateCSByWKT

    Then just repeat this process for the EPSG:27700 projected coordinate system, using WKT as follows:

    PROJCS["OSGB 1936 / British National Grid",  GEOGCS["OSGB 1936",    DATUM["OSGB 1936",      SPHEROID["Airy 1830", 6377563.396, 299.3249646, AUTHORITY["EPSG","7001"]],      TOWGS84[446.448, -125.157, 542.06, 0.15, 0.247, 0.842, -4.2261596151967575],      AUTHORITY["EPSG","6277"]],    PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],    UNIT["degree", 0.017453292519943295],    AXIS["Geodetic longitude", EAST],    AXIS["Geodetic latitude", NORTH],    AUTHORITY["EPSG","4277"]],  PROJECTION["Transverse Mercator"],  PARAMETER["central_meridian", -2.0],  PARAMETER["latitude_of_origin", 49.0],  PARAMETER["scale_factor", 0.9996012717],  PARAMETER["false_easting", 400000.0],  PARAMETER["false_northing", -100000.0],  UNIT["m", 1.0],  AXIS["Easting", EAST],  AXIS["Northing", NORTH],  AUTHORITY["EPSG","27700"]]

    Then, create the transformation between these two coordinate systems as normal.

     


    twitter: @alastaira blog: http://alastaira.wordpress.com/
    • Marked as answer by ckl42 Sunday, September 18, 2011 11:45 PM
    Sunday, September 18, 2011 3:42 PM
    Answerer
  • Every spatial reference system is based on an underlying model of the earth - the datum. The datum models used in different spatial reference systems may vary in size and shape.
    To convert coordinates defined using one spatial reference system into another, you need to account for any differences in the datums they use. The seven TOWGS84 parameters you list give the rotations (in 3 axes), translations (in 3 axes), and scaling factor that can be applied to transform from the OSGB36 datum to the WGS84 datum ("To WGS84" => "TOWGS84") using a transformation called the Helmert transformation, which is what Proj.NET (and many other programs use)

    If you didn't provide these TOWGS84 parameters, that would explain why your original result was 134m out because, without them, Proj.NET would not be able to perform the datum transformation and you'd get results defined on the same OSGB36 datum on which they were originally stated.

    Note that datum transformations (whether via the Helmert transformation or using other methods) are only ever approximate, and you will generally lose precision with every transformation. Any differences following the 6th d.p. are probably just inherent rounding errors introduced at some stage of the calculation.


    twitter: @alastaira blog: http://alastaira.wordpress.com/
    • Marked as answer by ckl42 Sunday, September 18, 2011 11:45 PM
    Sunday, September 18, 2011 6:51 PM
    Answerer

All replies

  • I think the point that you believe to be "134m out" is actually correct. Your problem lies in a misunderstanding, here:

    "....[the Ordnance Survey] has a spreadsheet that converts between OS easting/northing and WGS84 long/lat..."

    That spreadsheet does convert between OS easting/northing and long/lat, but it's not WGS84 long/lat - it's long/lat measured relative to the OSGB36 datum (i.e. the datum on which the National Grid is based.

    The EPSG for that system is 4277, and if you use Proj.NET with that reference you should get very similar results to in the spreadsheet.

    DECLARE @g geometry = geometry::Point(651409.903, 313177.271, 27700);
    SELECT Spatial.dbo.GeometryToGeography(@g, 4277).ToString();                       <-- POINT (1.7179216092420577 52.657570309931266)

    However, if you want WGS84 long/lat instead then use EPSG:4326, which gives the following:

    DECLARE @g geometry = geometry::Point(651409.903, 313177.271, 27700);
    SELECT Spatial.dbo.GeometryToGeography(@g, 4326).ToString();                   <-- POINT (1.7160520484753845 52.6579756312633)

    I imagine the reason for you thinking your results are 134m out is because, if you treat both these points as being defined relative to the WGS84 datum, they are 134m apart:

    declare @k geography = 'POINT (1.7179216092420577 52.657570309931266)';
    declare @j geography = 'POINT (1.7160520484753845 52.6579756312633)';
    select @k.STDistance(@j);         <-- 134.30742475263
    



    twitter: @alastaira blog: http://alastaira.wordpress.com/
    Friday, September 16, 2011 6:46 PM
    Answerer
  • Except that on the "UD Funcs Transformation Example" tab in the spreadsheet, it appears very clear that it is WGS84. It even refers to a handheld GPS receiver, which would be WGS84.

    Also, where is the GeometryToGeography function, I don't see it in ProjNet?

    Saturday, September 17, 2011 8:40 AM
  • That tab on the spreadsheet confirms exactly what I said - for the point in question:

    52.657570304134, 1.717921583102 is lat/long relative to the OSGB36 datum. 52.657978605604, 1.716051956076 is lat/long relative to the WGS84 datum.

    That spreadsheet is pretty horribly laid out, but look at the section on the "UD Funcs" tab beginning at row 55 that says OSGB36 to WGS84. There are 4 steps involved:

    • The first step is to convert from OSGB36 Eastings/Northings to OSGB36 latitude/longitude. This gives the values 52.657570304134, 1.717921583102 (cells R61C6 and R62C6)
    • You then convert to geocentric Cartesian X/Y/Z, which gives 3874938.849, 116218.624, 5047168.208
    • Then apply the Helmert transformation to get X/Y/Z relative to WGS84 - which gives 3875311.472, 116103.230, 5047602.299
    • Finally, convert back from X/Y/Z to geographic lat/long in WGS84, which gives 52.657978605604, 1.716051956076 (cells R98C6 and R99C6)

    twitter: @alastaira blog: http://alastaira.wordpress.com/
    Saturday, September 17, 2011 5:37 PM
    Answerer
  • I agree with everything you said there, but that was the point of my original question. The code I posted produces a result of 52.657570304134, 1.717921583102 when I am trying to produce  52.657978605604, 1.716051956076. I have tried all sorts of WKT and built-in conversions from ProjNet but haven't been able to get the right answer yet. I'm obviously missing a step somewhere, but can't work it out.

    Charles

    Sunday, September 18, 2011 11:20 AM
  • We're obviously misunderstanding each other somewhere.... so let's tackle this from a different approach.

    Why not try defining the source and target coordinate systems from which you're creating your coordinatetransformation factory using the CoordinateSystemWktReader.Parse() method instead? That way, you just have to pass in the Well-Known Text of each coordinate system rather than manually setting the projection parameters as you are doing currently.

    There's an example of how to define the WGS84 geographic coordinate system from WKT in the Proj.NET documentation, at the top of the following page: http://projnet.codeplex.com/wikipage?title=CreateCSByWKT

    Then just repeat this process for the EPSG:27700 projected coordinate system, using WKT as follows:

    PROJCS["OSGB 1936 / British National Grid",  GEOGCS["OSGB 1936",    DATUM["OSGB 1936",      SPHEROID["Airy 1830", 6377563.396, 299.3249646, AUTHORITY["EPSG","7001"]],      TOWGS84[446.448, -125.157, 542.06, 0.15, 0.247, 0.842, -4.2261596151967575],      AUTHORITY["EPSG","6277"]],    PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],    UNIT["degree", 0.017453292519943295],    AXIS["Geodetic longitude", EAST],    AXIS["Geodetic latitude", NORTH],    AUTHORITY["EPSG","4277"]],  PROJECTION["Transverse Mercator"],  PARAMETER["central_meridian", -2.0],  PARAMETER["latitude_of_origin", 49.0],  PARAMETER["scale_factor", 0.9996012717],  PARAMETER["false_easting", 400000.0],  PARAMETER["false_northing", -100000.0],  UNIT["m", 1.0],  AXIS["Easting", EAST],  AXIS["Northing", NORTH],  AUTHORITY["EPSG","27700"]]

    Then, create the transformation between these two coordinate systems as normal.

     


    twitter: @alastaira blog: http://alastaira.wordpress.com/
    • Marked as answer by ckl42 Sunday, September 18, 2011 11:45 PM
    Sunday, September 18, 2011 3:42 PM
    Answerer
  • Eureka! I have it ... with a lot of help :-). That works.

    I actually used

        IGeographicCoordinateSystem gcs = GeographicCoordinateSystem.WGS84;

    for the WGS84 coordinate system, and

        ICoordinateSystem coordsys = cFac.CreateFromWkt(epsg27700);

    using your WKT, above, and it producse pretty accurate results; by which I mean compared to the spreadsheet. I had tried something like that before, but one particular element I was missing was

        TOWGS84[446.448, -125.157, 542.06, 0.15, 0.247, 0.842, -4.2261596151967575]

    Could you explain what that does, exactly?

    Also - and your're going to think I'm being picky now - why aren't the answers I get now exactly the same as those from the spreadsheet? I had presumed that all the constants involved would be universal for these various conversions, but presumably they vary slightly? The answer I get converting OSGB36 to WGS84 is 52.6579756312633, 1.71605204847538. Working to 6 decimal places for long/lats the latitude is exact, but the longitude is 3 points out in the 6th place. Which 'constant' would be responsible for this variation, or could I tweak to bring the answers closer still?

    Charles

     

    Sunday, September 18, 2011 5:11 PM
  • Every spatial reference system is based on an underlying model of the earth - the datum. The datum models used in different spatial reference systems may vary in size and shape.
    To convert coordinates defined using one spatial reference system into another, you need to account for any differences in the datums they use. The seven TOWGS84 parameters you list give the rotations (in 3 axes), translations (in 3 axes), and scaling factor that can be applied to transform from the OSGB36 datum to the WGS84 datum ("To WGS84" => "TOWGS84") using a transformation called the Helmert transformation, which is what Proj.NET (and many other programs use)

    If you didn't provide these TOWGS84 parameters, that would explain why your original result was 134m out because, without them, Proj.NET would not be able to perform the datum transformation and you'd get results defined on the same OSGB36 datum on which they were originally stated.

    Note that datum transformations (whether via the Helmert transformation or using other methods) are only ever approximate, and you will generally lose precision with every transformation. Any differences following the 6th d.p. are probably just inherent rounding errors introduced at some stage of the calculation.


    twitter: @alastaira blog: http://alastaira.wordpress.com/
    • Marked as answer by ckl42 Sunday, September 18, 2011 11:45 PM
    Sunday, September 18, 2011 6:51 PM
    Answerer
  • Thanks Alastair, a very comprehensive answer.

    Charles

    Sunday, September 18, 2011 11:44 PM
  • Hi guys,

    Last night I have written an article so I hope it makes things clearer to understand:

    http://osedok.wordpress.com/2012/01/17/conversion-of-british-national-grid-wkid27700-to-wgs84wkid4326-and-then-to-webmercator-wkid102100/

    Regards,

    Andrzej

    Wednesday, January 18, 2012 1:42 AM
  • I tried several WKT combination already and they all gave the 130m error. Following these instructions exactly gives the same answer as ArcGIS - so this answer helped me too.
    Thursday, October 5, 2017 12:38 PM
  • Thank you for posting this. Surprisingly this was harder to find than I expected, and your response has enabled me to validate (and match) using the ProjNet NuGet against an old 32-bit conversion in-house routine we've built.
    Wednesday, June 6, 2018 12:25 PM