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

104

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

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

Дафшез
источник
1
Удалось ли вам преобразовать долготу и широту WGS-84 в декартовы координаты? У меня тоже есть возвышение. Я пробовал здесь принятый ответ, но он не дает мне правильного ответа. Я сравнил свои результаты с этим сайтом: apsalin.com/convert-geodetic-to-cartesian.aspx .
Ясмин

Ответы:

47

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

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

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

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

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

млрд.
источник
1
Ошибка 0,5% - 0,5% чего? В контексте этого вопроса это может быть радиус Земли, поэтому 0,5% может составлять 30 км :)
MarkJ
2
Проверил вашу ссылку. Цитата 0,5% предназначена для ошибки в расстоянии большого круга между двумя точками, поэтому не имеет прямого отношения к этому вопросу. Я думаю, что при преобразовании долготы и долготы в декартовы координаты с началом координат в центре Земли, ошибки от предположения о сферической Земле могут быть значительными. Непонятно, что спрашивающий хочет делать с декартовыми координатами. Либо просто в них удобнее работать по какой-то причудливой причине, либо это какое-то требование для экспорта данных? В последнем случае важна точность.
MarkJ
132

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

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

  • ось x проходит через долгую широту (0,0), так что долгота 0 пересекает экватор;
  • ось Y проходит через (0,90);
  • а ось z проходит через полюса.

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

x = R * cos(lat) * cos(lon)

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

z = R *sin(lat)

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

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

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

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

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

На этой странице приведен код C # для этого (обратите внимание, что он сильно отличается от формул), а также некоторые объяснения и красивую диаграмму того, почему это правильно,

Дафшез
источник
19
-1 Это неправильно. Вы предполагаете, что Земля - ​​это сфера, а WGS-84 - эллипсоид.
starblue
43
@starblue: Я не уверен, что вы в состоянии назвать данный ответ "правильным" или "неправильным". Сферическое приближение (для получения координат x, y, z в стиле ECEF) с использованием доступных значений широты и долготы (которые относятся к WGS-84) либо «адекватно» для нужд оригинального плаката, либо «неадекватно». Готов поспорить, что это простое преобразование подходит для оценки расстояния и пеленга. Если он запускает спутники, может и нет. В конце концов, WGS-84 сам по себе "неправильный" ... в том, что это не идеальная модель земной поверхности; все эллипсоидальные модели являются приближениями. Жаль, что OP не сказал нам, что он пытался сделать.
Dan H
11
@Dan H Вопрос касается WGS-84, и если вы ответите на что-то еще, вы должны хотя бы обсудить различия / ошибки, которых в этом ответе нет.
starblue 06
@ daphna-shezaf не может выполнить обратное преобразование ... Я также сделал обратное преобразование из радианов в градусы, но результат не тот ...
спасибо, потратьте час, выясняя, почему это не работает, оказывается, я поменял местами cos (лат) и sin (лат)
aeroson
6

Теория преобразования 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 и соавт, 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)
Хауи
источник
Какое значение R?
eych 04
4
Я предполагаю, что это радиус сферы, который составляет 6371 км для Земли.
Matthias
6

В python3.x это можно сделать, используя:

# Converting lat/long to cartesian
import numpy as np

def get_cartesian(lat=None,lon=None):
    lat, lon = np.deg2rad(lat), np.deg2rad(lon)
    R = 6371 # radius of the earth
    x = R * np.cos(lat) * np.cos(lon)
    y = R * np.cos(lat) * np.sin(lon)
    z = R *np.sin(lat)
    return x,y,z
Маянк Кумар
источник
5

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

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

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

Брайан Хокинс
источник
3

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

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

Степан Райко
источник
Голосование за, даже если вы не разместили здесь содержание ссылки.
Безумный физик
2

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

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

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

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

Юваль Адам
источник
1
+1. Использование надежной библиотеки точнее, чем функция домашнего приготовления, а также проще .
MarkJ
5
Как NetTopologySuite преобразуется из длинного / позднего в картографическое?
vinayan
1
NTS не включает в себя возможности преобразования координат, возможно, вам понадобится Proj.NET projnet.codeplex.com
D_Guidi
6
Нелепо, но ответ даже не предусматривает возможности преобразования.
Motes
1
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);
Ян Цзюнь
источник
Не могли бы вы уточнить? Я создал простое приложение, которое ярусами для преобразования одной координаты, используя ваш подход. Однако он всегда терпит неудачу, поскольку размеры источника (2) и размеры цели (3) различаются, что приводит к исключениюjava.lang.IllegalArgumentException: dimension must be <= 3
oschrenk
Хммм ... Я немного посмотрел на JTS. Строки до и включая новую LineString () выглядят как JTS. Но я не вижу CRS и Transform в JTS. Итак: они там, и я их скучаю? Были ли там и удалены в 1.12? Или: это другая библиотека?
Dan H
0

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

public List<Double> 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<Double> ecef= new ArrayList<>();
    ecef.add(x);
    ecef.add(y);
    ecef.add(z);

    return ecef;


}
Ашутош Чапагейн
источник
что такое параметр alt?
Baliman
высота, что ты вообще здесь делаешь, если не знаешь, как работает GPS;)
MushyPeas