webentwicklung-frage-antwort-db.com.de

Konvertierung von Längengrad\Breitengrad in kartesische Koordinaten

Ich habe einige erdzentrierte Koordinatenpunkte, die als Längen- und Breitengrad angegeben sind ( WGS-84 ).

Wie kann ich sie in kartesische Koordinaten (x, y, z) umwandeln, wobei der Ursprung im Erdmittelpunkt liegt?

86
daphshez

Ich habe kürzlich etwas Ähnliches mit der "Haversine Formula" auf WGS-84-Daten gemacht, eine Ableitung des "Law of Haversines" mit sehr befriedigenden Ergebnissen. 

Ja, WGS-84 geht davon aus, dass die Erde ein Ellipsoid ist, aber ich glaube, Sie erhalten nur einen durchschnittlichen Fehler von 0,5%, wenn Sie einen Ansatz wie die "Haversine-Formel" verwenden, was in Ihrem Fall eine akzeptable Fehlermenge sein kann. Sie werden immer ein gewisses Maß an Fehlern haben, es sei denn, Sie sprechen von ein paar Metern Entfernung, und sogar dann gibt es theoretisch eine Krümmung der Erde ... Wenn Sie einen starreren WGS-84-kompatiblen Ansatz benötigen, überprüfen Sie die "Vincenty-Formel".

Ich verstehe, wo starblue herkommt, aber bei guter Softwareentwicklung geht es oft um Kompromisse, daher hängt alles von der Genauigkeit ab, die Sie für das, was Sie tun, benötigen. Zum Beispiel kann das aus "Manhattan Distance Formula" berechnete Ergebnis gegenüber dem Ergebnis aus der "Distance Formula" für bestimmte Situationen besser sein, da es rechnerisch weniger teuer ist. Denken Sie "welcher Punkt ist am nächsten?" Szenarien, in denen Sie keine genaue Distanzmessung benötigen.

In Bezug auf die "Haversine-Formel" ist sie einfach zu implementieren und ist Nizza, da sie "Sphärische Trigonometrie" anstelle eines auf "Law of Cosines" basierenden Ansatzes verwendet, der auf zweidimensionaler Trigonometrie basiert. Daher erhalten Sie eine Nizza-Balance der Genauigkeit über Komplexität.

Ein Gentlemen mit dem Namen Chris Veness hat eine großartige Website unter http://www.movable-type.co.uk/scripts/latlong.html , die einige der Konzepte, an denen Sie interessiert sind, und demonstrieren verschiedene programmatische Implementierungen; Dies sollte auch Ihre X/Y-Konvertierungsfrage beantworten.

39
bn.

Hier ist die Antwort, die ich gefunden habe:

Um die Definition im kartesischen Koordinatensystem zu vervollständigen:

  • die x-Achse verläuft durch long, lat (0,0), sodass Längengrad 0 auf den Äquator trifft. 
  • die y-Achse geht durch (0,90);
  • und die z-Achse verläuft durch die Pole. 

Die Konvertierung ist:

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

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

z = R *sin(lat)

Wobei R der ungefähre Radius der Erde ist (z. B. 6371 KM).

Wenn Ihre trigonometrischen Funktionen Radiant erwarten (was wahrscheinlich der Fall ist), müssen Sie zuerst Längen- und Breitengrad in Radiant umrechnen. Offensichtlich benötigen Sie eine dezimale Darstellung, nicht Grad\Minuten\Sekunden (siehe hier zum Beispiel über die Konvertierung).

Die Formel für die Rückkonvertierung:

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

asin ist natürlich arc sinus. lese über atan2 in Wikipedia . Vergessen Sie nicht, von Radiant in Grad umzuwandeln. 

Diese Seite gibt den C # -Code dafür (beachten Sie, dass er sich sehr von den Formeln unterscheidet) und auch einige Erklärungen und ein Nizza-Diagramm, warum dies richtig ist.

107
daphshez

Theorie für die Umwandlung von GPS(WGS84) in Kartesische Koordinatenhttps://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_geodetic_to_ECEF_coordinates

Folgendes verwende ich:

  • Längengrad in GPS (WGS84) und kartesische Koordinaten sind die gleichen. 
  • Latitude muss mit den Ellipsoidparametern WGS 84 konvertiert werden, wobei die Semi-Major-Achse 6378137 m beträgt 
  • Kehrwert der Abflachung ist 298.257223563.

Ich habe einen VB Code angehängt, den ich geschrieben habe:

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

Bitte beachten Sie, dass die h die Höhe über dem WGS 84 ellipsoid liegt.

Normalerweise wird GPS uns H von über MSL Höhe geben. Die MSL-Höhe muss mithilfe des geopotential -Modells WGS 84 ellipsoid (Lemoine et al., 1998) in die Höhe h über dem EGM96 konvertiert werden.
Dies geschieht durch Interpolieren eines Gitters der Datei mit der Geoidhöhe mit einer räumlichen Auflösung von 15 Bogenminuten.

Oder wenn Sie über ein Niveau professionalGPS verfügen, das über die Höhe _H (msl, Höhe über dem mittleren Meeresspiegel) und UNDULATION die Beziehung zwischen der geoid und der ellipsoid (m) der gewählten Datumsausgabe von intern verfügt Tabelle. Sie können h = H(msl) + undulation bekommen 

Nach XYZ durch kartesische Koordinaten:

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

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

z = R *sin(lat)
6
Howie

Die proj.4 - Software stellt ein Befehlszeilenprogramm bereit, das die Umwandlung durchführen kann, z.

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

Es bietet auch eine C API . Insbesondere führt die Funktion pj_geodetic_to_geocentric die Konvertierung durch, ohne zuerst ein Projektionsobjekt einrichten zu müssen.

4
Brian Hawkins

Warum etwas implementieren, das bereits implementiert und getestet wurde?

C # hat zum Beispiel NetTopologySuite , den .NET-Port der JTS Topology Suite.

Insbesondere haben Sie einen schwerwiegenden Fehler in Ihrer Berechnung. Die Erde ist keine perfekte Kugel, und die Annäherung des Erdradius schneidet sie möglicherweise nicht für genaue Messungen.

Wenn es in einigen Fällen akzeptabel ist, Homebrew-Funktionen zu verwenden, ist GIS ein gutes Beispiel für ein Gebiet, in dem es sehr bevorzugt wird, eine zuverlässige, testbewährte Bibliothek zu verwenden.

4
Yuval Adam

Wenn Sie Koordinaten aus einem Ellipsoid und nicht aus einer Kugel erhalten möchten, schauen Sie unter http://en.wikipedia.org/wiki/Geodetic_system#From_geodetic_to_ECEF nach. Sie erhalten die Formeln sowie die WGS84-Konstanten Notwendigkeit für die Konvertierung.

Die dortigen Formeln berücksichtigen auch die Höhe relativ zur Referenzellipsoidoberfläche (nützlich, wenn Sie Höhendaten von einem GPS-Gerät erhalten).

3
Stjepan Rajko

In python3.x kann Folgendes verwendet werden:

# 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
1
Mayank Kumar
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);
1
yang jun

Sie können dies auf Java tun.

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;


}
0