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

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

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

• ось 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 приблизительный радиус Земли (например 6371КМ).

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

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

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

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

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

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

вот что я использую:

• долгота в GPS(WGS84) и декартовых координатах одинаковы.
• широта должна быть преобразована параметрами эллипсоида WGS 84 полу-большая ось 6378137 m, и
• взаимное сплющивание 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.

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

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

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

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

The 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 будет выполнять преобразование без необходимости сначала настроить объект проекции.

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<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;


}