Почему мой код C ++ создает неправильные пропорции / масштаб при преобразовании OSM WGS84 в ECEF и OSGB36? - PullRequest
Новая
       55

Почему мой код C ++ создает неправильные пропорции / масштаб при преобразовании OSM WGS84 в ECEF и OSGB36?

0 голосов
/ 11 мая 2019

Резюме

Я конвертирую координаты OSM в значения X и Y с точными расстояниями в метрах между всеми точками, используя c ++. Расчет / код, который я использую, приводит к неправильному отношению / шкале в x и y. Где я ошибся в своем расчете или коде? Я знаю, что есть ошибка в соотношении из-за кольцевых развязок, которые в действительности круглые показывают как эллиптические. Я знаю, что шкала неверна, потому что кольцевые дороги имеют меньшие размеры, чем они должны быть относительно британских данных DSM Lidar.

До сих пор я пытался выполнить шаг 1. Преобразуйте (геодезическую) широту / долготу в геоцентрические декартовы координаты ECEF отсюда https://www.movable -type.co.uk / scripts / latlong-os-gridref.html .

Я также проверил формулу по этому сайту https://gssc.esa.int/navipedia/index.php/Ellipsoidal_and_Cartesian_Coordinates_Conversion.

Я также пытался применить преобразование Гельмерта после первоначального преобразования ECEF, чтобы преобразовать данные, которые в большей степени соответствуют OSGB36.

Я пытался применить различные градусы к радианам только для моей широты, чтобы исправить проблему с соотношением сторон. Это работает в некоторой степени грубым образом, но не хватает точности и не решает основную проблему с моими вычислениями или масштабом.

У меня изначально была высота, равная 0, но мои результаты вычислений отличались от тех, которые были получены на http://www.apsalin.com/convert-geodetic-to-cartesian.aspx,, из-за пробной ошибки, я корректировал значение высоты в своем коде до тех пор, пока он не привел к результатам, которые соответствовали онлайн-конвертеру чему-то вроде 4 или 5 десятичных знаков.

WGS84 для декартовой ECEF просто (по моим результатам) представляет собой ортографическую проекцию координат WGS84, а не проекцию, которая компенсирует проецирование на плоскую поверхность. 

//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)
...