Built-in 7-param datums

Sep 20, 2014 at 7:24 PM
Edited Sep 20, 2014 at 7:26 PM
Hello,

today I was playing around with some datums (especially the OSG36 datum) and found the following issue:
Datum 1
# OSGB 1936 / British National Grid 
+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs
Datum 2
# OSGB 1936 / British National Grid 
+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy 
+towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894 +units=m +no_defs
These datum definitions should produce the same results, however they do not. When converting from OSGB36 X/Y to WGS84 LAT/LON, datum1 produces incorrect results, datum2 is OK.

Looking into the source code I found that when a 7-param datum is initialized through the proj.4 string, its towgs84 parameters are internally modified by the InitializeToWgs84 method:
// Transform from arc seconds to radians
_toWgs84[3] *= SEC_TO_RAD;
_toWgs84[4] *= SEC_TO_RAD;
_toWgs84[5] *= SEC_TO_RAD;
// transform from parts per millon to scaling factor
_toWgs84[6] = (_toWgs84[6] / 1000000.0) + 1;
However, the internal OSGB36 datum is initialized without these modifications to the parameters. If I modify the code to:
case "osgb36":
        _toWgs84 = new[] { 446.448, -125.157, 542.060, 0.1502, 0.2470, 0.8421, -20.4894 };
        _toWgs84[3] *= SEC_TO_RAD;
        _toWgs84[4] *= SEC_TO_RAD;
        _toWgs84[5] *= SEC_TO_RAD;
        _toWgs84[6] = (_toWgs84[6] / 1000000.0) + 1;
        _spheroid = new Spheroid(Proj4Ellipsoid.Airy_1830);
        _description = "Airy 1830";
        _datumtype = DatumType.Param7;
        break;
then datum1 and datum2 produce the same results.

Could somebody please confirm my findings? I think the same applies for the internal ire65 and nzgd49 datums.