Conversion de longitude latitude en coordonnées cartésiennes

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

Juste pour faire de la définition complète, dans le système de coordonnées Cartésiennes:

• l'axe des x passe par long, lat (0,0), donc la longitude 0 rencontre l'Équateur;
• l'axe des y passe par (0,90);
• et l'axe z passe à travers 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 rayon approximatif de la terre (par exemple 6371KM).

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

la formule pour la conversion en arrière:

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

asin est de cours arc sinus. lire à propos de atan2 dans wikipédia . N'oubliez pas de convertir en degrés le retour de radians.

Cette page donne le code c# (à noter qu'il est très différent de la formule), et aussi une explication et une belle diagramme de pourquoi c'est correct,

théorie pour convertir GPS（WGS84) en coordonnées cartésiennes https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_geodetic_to_ECEF_coordinates

ce qui suit est ce que j'utilise:

• la Longitude dans le GPS WGS84 le WGS84 le WGS84 le WGS84 le WGS84 et les coordonnées cartésiennes sont les mêmes.
• Latitude doivent être convertis par WGS 84 paramètres ellipsoïdes demi-grand axe est 6378137 m, et
• réciproque d'aplatissement est 298.257223563.

j'ai joint un code VB 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 le h est une altitude au-dessus du WGS 84 ellipsoid .

habituellement GPS nous donnera H de plus de MSL hauteur. La hauteur MSL doit être converti en hauteur h au-dessus du WGS 84 ellipsoid en utilisant le géopotentiel modèle EGM96 ( Lemoine et al, 1998 ).

Ceci est fait en interpolant une grille du fichier de hauteur géoïde avec une résolution spatiale de 15 arc-minutes.

ou si vous avez un certain niveau professionnel GPS a de L'Altitude H ( msl,heigh au-dessus de la mer moyenne niveau ) et UNDULATION , la relation entre le geoid et le ellipsoid (m) de la sortie choisie sortie de référence de la table interne. vous pouvez obtenir h = H(msl) + undulation

À XYZ par coordonnées cartésiennes:

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

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

z = R *sin(lat)

Pourquoi mettre en œuvre quelque chose qui a déjà été mis en œuvre et testé?

C#, pour un, a la NetTopologySuite qui est le .NET port de la STC de la Topologie de Suite.

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

si dans dans certains cas, il est acceptable d'utiliser des fonctions maison, GIS est un bon exemple d'un domaine dans lequel il est préférable d'utiliser une bibliothèque fiable, testée.

si vous vous souciez d'obtenir des coordonnées basées sur un ellipsoïde plutôt qu'une sphère, jetez un oeil à 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 prennent également en compte l'altitude relative à la surface ellipsoïde de référence (utile si vous obtenez des données d'altitude à partir d'un GPS).

Le proj.4 logiciel fournit un programme de 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 un C API . En particulier, la fonction pj_geodetic_to_geocentric effectuera la conversion sans avoir à configurer d'abord un objet de projection.

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;


}