Conversion de longitude \ latitude en coordonnées cartésiennes

J'ai quelques points de coordonnées centrés sur la terre donnés en latitude et longitude ( WGS-84 ).

Comment puis-je les convertir en coordonnées cartésiennes (x, y, z) avec l'origine au centre de la terre?

J'ai récemment fait quelque chose de similaire en utilisant la "formule Haversine" sur les données WGS-84, qui est un dérivé de la "loi des havresins" avec des résultats très satisfaisants.

Oui, WGS-84 suppose que la Terre est un ellipsoïde, mais je crois que vous n'obtenez qu'une erreur moyenne d'environ 0,5% en utilisant une approche comme la «formule Haversine», ce qui peut être une quantité d'erreur acceptable dans votre cas. Vous aurez toujours une certaine quantité d'erreur à moins que vous ne parliez d'une distance de quelques mètres et même alors il y a théoriquement une courbure de la Terre ... Si vous avez besoin d'une approche plus rigoureusement compatible WGS-84, consultez la "formule Vincenty".

Je comprends d'où vient starblue , mais une bonne ingénierie logicielle est souvent une question de compromis, donc tout dépend de la précision dont vous avez besoin pour ce que vous faites. Par exemple, le résultat calculé à partir de la «formule de distance de Manhattan» par rapport au résultat de la «formule de distance» peut être meilleur dans certaines situations car il est moins coûteux en calcul. Pensez "quel point est le plus proche?" scénarios où vous n'avez pas besoin d'une mesure de distance précise.

En ce qui concerne la "formule Haversine", elle est facile à mettre en œuvre et est agréable car elle utilise la "trigonométrie sphérique" au lieu d'une approche basée sur la "loi des cosinus" qui est basée sur la trigonométrie bidimensionnelle, vous obtenez donc un bel équilibre de précision sur la complexité.

Un gentlemen du nom de Chris Veness a un excellent site Web à http://www.movable-type.co.uk/scripts/latlong.html qui explique certains des concepts qui vous intéressent et montre diverses implémentations programmatiques; cela devrait également répondre à votre question de conversion x / y.

Voici la réponse que j'ai trouvée:

Juste pour rendre la définition complète, dans le système de coordonnées cartésien:

• l'axe des x passe par long, lat (0,0), donc la longitude 0 rencontre l'équateur;
• l'axe y passe par (0,90);
• et l'axe z passe par les pôles.

La conversion est:

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

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

z = R *sin(lat)

Où R est le rayon approximatif de la Terre (par exemple, 6371 km).

Si vos fonctions trigonométriques attendent des radians (ce qu'elles font probablement), vous devrez d'abord convertir votre longitude et votre latitude en radians. Vous avez évidemment besoin d'une représentation décimale, pas de degrés \ minutes \ secondes (voir par exemple ici à propos de la conversion).

La formule de la rétro-conversion:

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

asin est bien sûr un arc sinus. lisez atan2 sur wikipedia . N'oubliez pas de convertir les radians en degrés.

Cette page donne du code c # pour cela (notez que c'est très différent des formules), ainsi que des explications et un joli diagramme expliquant pourquoi c'est correct,

Théorie de la conversion GPS（WGS84)en coordonnées cartésiennes https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_geodetic_to_ECEF_coordinates

Voici ce que j'utilise:

• La longitude en GPS （WGS84） et les coordonnées cartésiennes sont identiques.
• La latitude doit être convertie par les paramètres de l'ellipsoïde WGS 84.Le demi-grand axe est 6378137 m, et
• La réciproque de l'aplatissement est 298,257223563.

J'ai joint un code VB que j'ai écrit:

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

Veuillez noter que l' haltitude est au-dessus du WGS 84 ellipsoid.

Habituellement GPS, nous donnera une hauteur Hsupérieure MSL. La MSLhauteur doit être convertie en hauteur hau-dessus de la WGS 84 ellipsoiden utilisant le modèle géopotentielEGM96 ( Lemoine et al, 1998 ).
Cela se fait en interpolant une grille du fichier de hauteur du géoïde avec une résolution spatiale de 15 minutes d'arc.

Ou si vous avez un niveau professionnel GPS a Altitude H( msl, hauteur au-dessus du niveau moyen de la mer ) et UNDULATION, la relation entre le geoidet le ellipsoid (m)de la sortie de référence choisie de la table interne. vous pouvez obtenirh = H(msl) + undulation

Vers XYZ par coordonnées cartésiennes:

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

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

z = R *sin(lat)

Le logiciel proj.4 fournit un programme en ligne de commande qui peut faire la conversion, par exemple

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

Il fournit également une API C . En particulier, la fonction pj_geodetic_to_geocentriceffectuera la conversion sans avoir à configurer au préalable un objet de projection.

Dans python3.x, cela peut être fait en utilisant:

# 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

Si vous souhaitez obtenir des coordonnées basées sur un ellipsoïde plutôt qu'une sphère, jetez un œil à http://en.wikipedia.org/wiki/Geodetic_system#From_geodetic_to_ECEF - il donne les formules ainsi que les constantes WGS84 dont vous avez besoin pour la conversion .

Les formules y prennent également en compte l'altitude par rapport à la surface de l'ellipsoïde de référence (utile si vous obtenez des données d'altitude à partir d'un appareil GPS).

Pourquoi implémenter quelque chose qui a déjà été implémenté et testé?

C #, pour sa part, a la NetTopologySuite qui est le port .NET de JTS Topology Suite.

Plus précisément, vous avez un grave défaut dans votre calcul. La terre n'est pas une sphère parfaite et l' approximation du rayon de la terre pourrait ne pas la couper pour des mesures précises.

Si dans certains cas il est acceptable d'utiliser des fonctions homebrew, le SIG est un bon exemple d'un domaine dans lequel il est de loin préférable d'utiliser une bibliothèque fiable et éprouvée.

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

Vous pouvez le faire de cette façon sur 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;


}