Преобразование из долготы \ широты в декартовы координаты

У меня есть координатно-ориентированные координаты с землей, заданные как широта и долгота ( WGS-84 ).

Как я могу преобразовать их в декартовы координаты (x, y, z) с началом в центре Земли?

    Недавно я сделал что-то подобное этому, используя «Формулу Хаверсина» по данным WGS-84, которая является производным от «Закона Хаверсинеса» с очень удовлетворительными результатами.

    Да, WGS-84 предполагает, что Земля является эллипсоидом, но я полагаю, что вы получаете только среднюю ошибку 0,5%, используя подход, подобный «Формуле Хаверсина», которая может быть приемлемой суммой ошибок в вашем случае. У вас всегда будет некоторая ошибка, если вы не говорите о расстоянии в несколько футов, и даже тогда теоретически кривизна Земли … Если вам потребуется более жесткая проверка соответствия с WGS-84 на «Формуле Vincenty».

    Я понимаю, откуда исходит звездная звезда , но хорошая разработка программного обеспечения часто связана с компромиссом, поэтому все зависит от точности, которую вы требуете от того, что вы делаете. Например, результат, рассчитанный по формуле «Формула расстояния Манхэттена» по сравнению с результатом «Формулы расстояния», может быть лучше для определенных ситуаций, поскольку он является расчетно менее дорогостоящим. Подумайте «какой момент ближе?» где вам не требуется точное измерение расстояния.

    Что касается «Формулы Хаверсина», ее легко реализовать и хорошо, потому что она использует «сферическую тригонометрию» вместо подхода «Закон косинуса», основанного на двумерной тригонометрии, поэтому вы получаете хороший баланс точности по сложности.

    Джентльмены по имени Крис Винес имеют отличный веб-сайт по адресу http://www.movable-type.co.uk/scripts/latlong.html, который объясняет некоторые интересующие вас концепции и демонстрирует различные программные реализации; это также должно ответить на ваш вопрос об изменении x / y.

    Вот ответ, который я нашел:

    Чтобы сделать определение полным, в декартовой системе координат:

    • ось x проходит через long, lat (0,0), поэтому долгота 0 встречается с экватором;
    • ось y проходит (0,90);
    • и ось z проходит через полюсы.

    Преобразование:

    x = R * cos(lat) * cos(lon) y = R * cos(lat) * sin(lon) z = R *sin(lat) 

    Где R – приблизительный радиус земли (например, 6371KM).

    Если ваши тригонометрические функции ожидают радианы (что они, вероятно, делают), вам нужно будет сначала преобразовать свою долготу и широту в радианы. Вам, очевидно, нужно десятичное представление, а не gradleусы \ минут \ секунд (см., Например, здесь об конверсии).

    Формула для обратного преобразования:

      lat = asin(z / R) lon = atan2(y, x) 

    asin – это, конечно, синус. Читайте о atan2 в Википедии . Не забудьте перевести назад из радианов в gradleусы.

    Эта страница дает код c # для этого (обратите внимание, что она сильно отличается от формул), а также некоторое объяснение и хорошая диаграмма, почему это правильно,

    Теория преобразования GPS(WGS84) в декартовы координаты https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_geodetic_to_ECEF_coordinates

    Ниже приводится то, что я использую:

    • Долгота в GPS (WGS84) и декартовы координаты одинаковы.
    • Широта должна быть преобразована по параметрам эллипсоида WGS 84, полумалая ось 6378137 м и
    • Взаимное выравнивание равно 298.257223563.

    Я прикрепил код VB, который я написал:

     Imports System.Math 'Input GPSLatitude is WGS84 Latitude,h is altitude above the WGS 84 ellipsoid Public Function GetSphericalLatitude(ByVal GPSLatitude As Double, ByVal h As Double) As Double Dim A As Double = 6378137 'semi-major axis Dim f As Double = 1 / 298.257223563 '1/f Reciprocal of flattening Dim e2 As Double = f * (2 - f) Dim Rc As Double = A / (Sqrt(1 - e2 * (Sin(GPSLatitude * PI / 180) ^ 2))) Dim p As Double = (Rc + h) * Cos(GPSLatitude * PI / 180) Dim z As Double = (Rc * (1 - e2) + h) * Sin(GPSLatitude * PI / 180) Dim r As Double = Sqrt(p ^ 2 + z ^ 2) Dim SphericalLatitude As Double = Asin(z / r) * 180 / PI Return SphericalLatitude End Function 

    Обратите внимание, что h – высота над WGS 84 ellipsoid .

    Обычно GPS даст нам H выше высоты MSL . Высота MSL должна быть преобразована в высоту h над WGS 84 ellipsoid с использованием геопотенциальной модели EGM96 ( Lemoine et al, 1998 ).
    Это делается путем интерполяции сетки файла высоты геоида с пространственным разрешением 15 минут дуги.

    Или, если у вас есть определенный профессиональный GPS , высота H ( msl, высота выше среднего уровня моря ) и UNDULATION , связь между geoid и ellipsoid (m) выбранной исходной точки, выводимой из внутренней таблицы. вы можете получить h = H(msl) + undulation

    К XYZ декартовыми координатами:

     x = R * cos(lat) * cos(lon) y = R * cos(lat) * sin(lon) z = R *sin(lat) 

    Зачем внедрять что-то, что уже было реализовано и проверено на тестирование?

    C #, например, имеет NetTopologySuite, который является портом .NET для топологии серии JTS.

    В частности, у вас есть серьезный недостаток в ваших расчетах. Земля не является идеальной сферой, и приближение радиуса Земли не может вырезать ее для точных измерений.

    Если в некоторых случаях приемлемо использовать функции доморощенного, GIS является хорошим примером поля, в котором гораздо предпочтительнее использовать надежную, проверенную временем библиотеку.

    Если вы хотите получить координаты на основе эллипсоида, а не сферы, посмотрите на http://en.wikipedia.org/wiki/Geodetic_system#From_geodetic_to_ECEF – он дает формулы, а также константы WGS84, необходимые для преобразования ,

    Формулы, которые также принимают во внимание высоту по отношению к поверхности референц-эллипсоида (полезно, если вы получаете данные высоты от устройства GPS).

    Программное обеспечение proj.4 предоставляет программу командной строки, которая может выполнять преобразование, например

     LAT=40 LON=-110 echo $LON $LAT | cs2cs +proj=latlong +datum=WGS84 +to +proj=geocent +datum=WGS84 

    Он также предоставляет API C. В частности, функция pj_geodetic_to_geocentric выполнит преобразование без предварительной настройки проекционного объекта.

     Coordinate[] coordinates = new Coordinate[3]; coordinates[0] = new Coordinate(102, 26); coordinates[1] = new Coordinate(103, 25.12); coordinates[2] = new Coordinate(104, 16.11); CoordinateSequence coordinateSequence = new CoordinateArraySequence(coordinates); Geometry geo = new LineString(coordinateSequence, geometryFactory); CoordinateReferenceSystem wgs84 = DefaultGeographicCRS.WGS84; CoordinateReferenceSystem cartesinaCrs = DefaultGeocentricCRS.CARTESIAN; MathTransform mathTransform = CRS.findMathTransform(wgs84, cartesinaCrs, true); Geometry geo1 = JTS.transform(geo, mathTransform); 

    Вы можете сделать это на Java.

     public List convertGpsToECEF(double lat, double longi, float alt) { double a=6378.1; double b=6356.8; double N; double e= 1-(Math.pow(b, 2)/Math.pow(a, 2)); N= a/(Math.sqrt(1.0-(e*Math.pow(Math.sin(Math.toRadians(lat)), 2)))); double cosLatRad=Math.cos(Math.toRadians(lat)); double cosLongiRad=Math.cos(Math.toRadians(longi)); double sinLatRad=Math.sin(Math.toRadians(lat)); double sinLongiRad=Math.sin(Math.toRadians(longi)); double x =(N+0.001*alt)*cosLatRad*cosLongiRad; double y =(N+0.001*alt)*cosLatRad*sinLongiRad; double z =((Math.pow(b, 2)/Math.pow(a, 2))*N+0.001*alt)*sinLatRad; List ecef= new ArrayList<>(); ecef.add(x); ecef.add(y); ecef.add(z); return ecef; } 
    Давайте будем гением компьютера.