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

103

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?

daphshez
la source
1
Avez-vous réussi à convertir la longitude et la latitude WGS-84 en coordonnées cartésiennes?. J'ai aussi l'élévation. J'ai essayé la réponse acceptée ici, mais cela ne me donne pas la bonne réponse. J'ai comparé mes résultats avec ce site: apsalin.com/convert-geodetic-to-cartesian.aspx .
Yasmin

Réponses:

47

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.

bn.
la source
1
0,5% d'erreur - 0,5% de quoi? Dans le contexte de cette question, il pourrait s'agir du rayon de la terre, donc 0,5% pourrait être de 30 km :)
MarkJ
2
J'ai vérifié votre lien. La citation de 0,5% correspond à une erreur dans la distance du grand cercle entre deux points et n'est donc pas strictement pertinente pour cette question. Je pense que lors de la conversion lat-long en coordonnées cartésiennes avec l'origine au centre de la terre, les erreurs de supposer une terre sphérique pourraient être significatives. On ne sait pas ce que le questionneur veut faire des coordonnées cartésiennes. Soit il est simplement plus pratique d'y travailler pour une raison étrange, soit c'est peut-être une exigence pour l'exportation de données? Dans ce dernier cas, la précision serait importante.
MarkJ
130

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,

daphshez
la source
18
-1 C'est faux. Vous supposez que la terre est une sphère, tandis que WGS-84 suppose un ellipsoïde.
starblue
42
@starblue: Je ne suis pas sûr que vous soyez en mesure d'étiqueter la réponse donnée "bonne" ou "mauvaise". L'approximation sphérique (pour obtenir les coordonnées x, y, z de style ECEF) en utilisant les lat / lngs disponibles (qui sont référencés à WGS-84) est soit "adéquate" pour les besoins de l'affiche originale, soit "non adéquate". Pour les estimations de distance et de relèvement, je parie que cette simple conversion est très bien. S'il lance des satellites, peut-être pas. Après tout, le WGS-84 lui-même est "faux" ... en ce qu'il n'est pas un modèle parfait de la surface de la terre; tous les modèles ellipsoïdaux sont des approximations. Dommage que l'OP ne nous ait pas dit ce qu'il essayait de faire.
Dan H du
11
@Dan H La question demande WGS-84, et si vous répondez à autre chose, vous devriez au moins discuter des différences / erreurs, ce que cette réponse ne fait pas.
starblue
@ daphna-shezaf ne peut pas faire une conversion arrière ... J'ai aussi fait le retour des radians aux degrés, mais le résultat n'est pas le même ...
merci, passez une heure à comprendre pourquoi cela ne fonctionne pas, il s'avère que j'ai échangé certains cos (lat) et sin (lat)
aeroson
6

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)
Howie
la source
Quelle est la valeur de R?
eych
4
Je suppose que c'est le rayon de la sphère, qui est de 6371 km pour la terre.
Matthias
5

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.

Brian Hawkins
la source
5

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
Mayank Kumar
la source
3

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

Stjepan Rajko
la source
Vote positif même si vous n'avez pas publié le contenu du lien ici.
Mad Physicist
2

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.

Yuval Adam
la source
1
+1. Utiliser une bibliothèque fiable est plus précis qu'une fonction homebrew et aussi plus facile .
MarkJ
5
Comment NetTopologySuite convertit-il de long / tard en cartésion?
vinayan
1
NTS n'inclut pas de capacités de conversion de coordonnées, peut-être avez-vous besoin de Proj.NET projnet.codeplex.com
D_Guidi
6
Ridicule, la réponse ne fournit même pas de capacité de conversion.
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);
yang jun
la source
Seriez-vous en mesure d'élaborer? J'ai créé une application simple qui hiérarchise pour transformer une seule coordonnée en utilisant votre approche. Cela échoue toujours car les dimensions de la source (2) et les dimensions de la cible (3) diffèrent, ce qui entraîne une exceptionjava.lang.IllegalArgumentException: dimension must be <= 3
oschrenk
Hmmm ... J'ai regardé JTS pendant un moment. Les lignes jusqu'à et y compris la nouvelle LineString () ressemblent à JTS. Mais je ne vois pas les éléments CRS et Transform dans JTS. Alors: sont-ils là et ils me manquent? Y avait-il, et supprimé dans 1.12? Ou: est-ce une bibliothèque différente?
Dan H
0

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;


}
Ashutosh Chapagain
la source
qu'est-ce que le paramètre alt?
baliman
altitude, que faites-vous même ici si vous ne savez pas comment fonctionne le GPS;)
MushyPeas