Answered by:
Converting OSGB36 to WGS84 with ProjNet; Doing My Head In
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 builtin 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 WellKnown 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 PMAnswerer 
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 PMAnswerer
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 PMAnswerer 
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 PMAnswerer 
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 builtin 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 WellKnown 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 PMAnswerer 
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 PMAnswerer 
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:
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 32bit conversion inhouse routine we've built.Wednesday, June 6, 2018 12:25 PM