1

Summary

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)
  • The general concept you're reaching for is called "projection". You need to transform your coordinates from WGS84 to a target projection. Can you give more information about the X/Y/Z target coordinate system? Once identified, typically you'd use tools that do the correct mathematical transformations for you. – Richard Law Apr 29 at 22:37
  • Hi Richard, the target coordinate space is simply 3D space X, Y, Z in cm relative to a central point zero. I just need all of the coordinate points in the Osm to be placed at accurate X and Y distance in meters from each other. Because the data is in osm xml I'll probably need to do the calculations in my own c++ osm parser. – James Apr 29 at 23:08
  • PROJ4 can be used with C++, and I'd personally lean on that. I think you'd need to define a custom projection (e.g. gis.stackexchange.com/questions/97589/…) where you give it the central location, extent, vertical datum, etc. I've never had to actually do that so I'll leave it to someone who has to compose a real answer, but hopefully that's helpful. – Richard Law Apr 30 at 0:28
  • Thanks Richard, it's good to have options. I'm suspecting that my original code would yeild a less banana shaped projection if I used the correct number of meters per degree fo longitude and latitude above, as I live in England. I could use maxlon-(maxlon-minlon), and maxlat-(maxlat-minlat) and use the middle of the map lat and lon find degrees per meter at that point. It's the issue of trying to get that degrees per meter for an ellipsoid. – James Apr 30 at 7:23
  • If the target coordinate space is as @James describes then there is no need for a projection. Simply convert the lat/lon values to XYZ Earth Centred Earth Fixed (ECEF). full maths here Assume input height of 0 if you only have lat/lon values available. – JimT Apr 30 at 7:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.