I am converting OSM coordinates into X and Y values with accurate distances in meters between all points using c++. The calculations/code I am using is resulting in an incorrect ratio/scale in x and y. Where have I gone wrong in my calculation or code? I know there is an error in ratio due to roundabouts that are in reality round are showing as elliptical. I know the scale is wrong because ring roads are showing as being smaller than they should be relative to UK DSM Lidar data.
What I have tried so far is step 1. Convert (geodetic) latitude/longitude to ECEF geocentric cartesian coordinates from here https://www.movable-type.co.uk/scripts/latlong-os-gridref.html.
I have also checked the formula against this site https://gssc.esa.int/navipedia/index.php/Ellipsoidal_and_Cartesian_Coordinates_Conversion.
I have also tried applying the Helmert transform after the initial ECEF conversion to transform being the data more in line with OSGB36.
I have tried applying a different degrees to radians calculation only for my lattitude to correct the aspect ratio issue. This works to an extent in a rough way, but lacks accuracy and doesn't solve the underlying problem with my calculations or the scale.
I originally had height set as 0, but my calculation results differed from those on http://www.apsalin.com/convert-geodetic-to-cartesian.aspx, through trial an error, I adjusted the height value in my code until it produced results that matched the online converter to something like 4 or 5 decimal places.
WGS84 to Cartesian ECEF just appears (from my results) to be an orthographic projection of the WGS84 coordinates, rather than a projection that compensates for projecting onto a flat surface to give accurate distance between points.
(I have completely changed the way I have worded and presented the original question to reflect the further work I have been doing so solve my issue.)
//Constants for WGS84 to ECEF double height = (27716.480814599*2); double SemiMajorAxisA = 6378137.0; double RadiusOfCurvature; double FlatteningFactoroftheEarth = (1/298.257223563); double EccentricitySquared = 0.00669437999014; const double DEG_TO_RAD = 0.01745329251994330; const double SEC_TO_RAD = 0.000004848; //Helmerts constants for OSGB86// double cx = -446.448; double cy = 125.157; double cz = -542.06; double s = 0.0000204894; double rx = -0.1502*SEC_TO_RAD; double ry = -0.247*SEC_TO_RAD; double rz = -0.8421*SEC_TO_RAD; //WGS84 to ECEF for maximum lat and lon MaxLatdouble = MaxLatdouble * DEG_TO_RAD; MaxLondouble = MaxLondouble * DEG_TO_RAD; RadiusOfCurvature = SemiMajorAxisA / (sqrt(1 + (EccentricitySquared * (sin(MaxLatdouble) * sin(MaxLatdouble) )))); double MaxLatdoubleA = (RadiusOfCurvature + height) *(cos(MaxLatdouble) * cos(MaxLondouble)); double MaxLondoubleA = (RadiusOfCurvature + height)* (cos(MaxLatdouble) * sin(MaxLondouble)); //Apply Helmert transform to max lat and lon MaxLatdoubleA = cx + (MaxLatdoubleA * (1 + s)) + (-rz * MaxLondoubleA) + (ry * height); MaxLondoubleA = cy + (rz * MaxLatdoubleA) + (MaxLondoubleA * (1 + s)) + (-rx * height); //Now perform the same calculations used above for each OSM Lat Lon. // .. same code as above in a loop for each coordinate // final lat and lon for each is (maxlat - indivudual lat) and (maxlon - individual lon)