Conversión de longitud \ latitud a coordenadas cartesianas

103

Tengo algunos puntos de coordenadas centrados en la Tierra dados como latitud y longitud ( WGS-84 ).

¿Cómo puedo convertirlos a coordenadas cartesianas (x, y, z) con el origen en el centro de la tierra?

daphshez
fuente
1
¿Ha logrado convertir la longitud y latitud WGS-84 en coordenadas cartesianas? También tengo elevación. He probado la respuesta aceptada aquí, pero no me da la respuesta correcta. Comparé mis resultados con este sitio web: apsalin.com/convert-geodetic-to-cartesian.aspx .
Yasmin

Respuestas:

46

Recientemente he hecho algo similar a esto usando la "Fórmula de Haversine" en datos de WGS-84, que es un derivado de la "Ley de Haversines" con resultados muy satisfactorios.

Sí, WGS-84 asume que la Tierra es un elipsoide, pero creo que solo obtienes un error promedio de 0.5% usando un enfoque como la "Fórmula Haversine", que puede ser una cantidad aceptable de error en tu caso. Siempre tendrá una cierta cantidad de error a menos que esté hablando de una distancia de unos pocos pies e incluso entonces hay una curvatura teórica de la Tierra ... Si necesita un enfoque más rígidamente compatible con WGS-84, consulte la "Fórmula Vincenty".

Entiendo de dónde viene starblue , pero una buena ingeniería de software a menudo se trata de compensaciones, por lo que todo depende de la precisión que requiera para lo que está haciendo. Por ejemplo, el resultado calculado a partir de la "Fórmula de distancia de Manhattan" frente al resultado de la "Fórmula de distancia" puede ser mejor para determinadas situaciones, ya que es computacionalmente menos costoso. Piense en "¿qué punto está más cerca?" escenarios en los que no necesita una medición de distancia precisa.

En cuanto a la "Fórmula de Haversine", es fácil de implementar y es agradable porque utiliza "Trigonometría esférica" ​​en lugar de un enfoque basado en la "Ley de los cosenos" que se basa en la trigonometría bidimensional, por lo que se obtiene un buen equilibrio de precisión. sobre la complejidad.

Un caballero llamado Chris Veness tiene un gran sitio web en http://www.movable-type.co.uk/scripts/latlong.html que explica algunos de los conceptos que le interesan y demuestra varias implementaciones programáticas; esto también debería responder a su pregunta de conversión x / y.

bn.
fuente
1
0.5% de error - 0.5% de qué? En el contexto de esta pregunta, podría ser el radio de la tierra, por lo que 0.5% podría ser 30 km :)
MarkJ
2
Revisó su enlace. La cotización del 0,5% es por error en la distancia del círculo máximo entre dos puntos, por lo que no es estrictamente relevante para esta pregunta. Yo pensaría que al convertir lat-long a coordenadas cartesianas con el origen en el centro de la tierra, los errores de asumir una tierra esférica podrían ser significativos. No está claro qué quiere hacer el interrogador con las coordenadas cartesianas. ¿O es más conveniente trabajar en ellos por alguna extraña razón, o tal vez es algún requisito para la exportación de datos? Si es lo último, la precisión sería importante.
MarkJ
130

Aquí está la respuesta que encontré:

Solo para completar la definición, en el sistema de coordenadas cartesianas:

  • el eje x pasa por long, lat (0,0), por lo que la longitud 0 se encuentra con el ecuador;
  • el eje y pasa por (0,90);
  • y el eje z atraviesa los polos.

La conversión es:

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

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

z = R *sin(lat)

Donde R es el radio aproximado de la Tierra (por ejemplo, 6371 km).

Si sus funciones trigonométricas esperan radianes (lo que probablemente es así), primero deberá convertir su longitud y latitud a radianes. Obviamente, necesita una representación decimal, no grados \ minutos \ segundos (consulte, por ejemplo, aquí sobre la conversión).

La fórmula para la conversión inversa:

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

asin es, por supuesto, arc sine. leer acerca de atan2 en wikipedia . No olvide convertir de radianes a grados.

Esta página proporciona el código c # para esto (tenga en cuenta que es muy diferente de las fórmulas), y también una explicación y un buen diagrama de por qué esto es correcto.

daphshez
fuente
17
-1 Esto está mal. Está asumiendo que la tierra es una esfera, mientras que WGS-84 asume un elipsoide.
starblue
42
@starblue: No estoy seguro de que estés en posición de etiquetar la respuesta dada como "correcta" o "incorrecta". La aproximación esférica (para obtener las coordenadas x, y, z de estilo ECEF) usando lat / lngs disponibles (que se refieren a WGS-84) es "adecuada" para las necesidades del póster original o "no adecuada". Para estimaciones de distancia y rumbo, apuesto a que esta simple conversión está bien. Si está lanzando satélites, tal vez no. Después de todo, el propio WGS-84 está "equivocado" ... en el sentido de que no es un modelo perfecto de la superficie terrestre; todos los modelos elipsoidales son aproximaciones. Lástima que el OP no nos dijo lo que estaba tratando de hacer.
Dan H
11
@Dan H La pregunta solicita WGS-84, y si responde algo más, al menos debería discutir las diferencias / errores, lo que esta respuesta no hace.
starblue
@ daphna-shezaf no puede hacer una conversión inversa ... También hice la conversión de radianes a grados, pero el resultado no es el mismo ...
gracias, pasa una hora averiguando por qué no funciona, resulta que cambié algo de cos (lat) y sin (lat)
aeroson
6

Teoría para convertir GPS(WGS84)a coordenadas cartesianas https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_geodetic_to_ECEF_coordinates

Lo siguiente es lo que estoy usando:

  • La longitud en GPS (WGS84) y las coordenadas cartesianas son las mismas.
  • La latitud debe ser convertida por los parámetros del elipsoide WGS 84, el eje semi-mayor es 6378137 m
  • El recíproco de aplanamiento es 298,257223563.

Adjunté un código VB que escribí:

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

Tenga en cuenta que la haltitud es superior a WGS 84 ellipsoid.

Normalmente GPSnos dará una altura Hsuperior MSL. La MSLaltura debe convertirse en altura por hencima de la WGS 84 ellipsoidutilizando el modelo geopotencialEGM96 ( Lemoine et al, 1998 ).
Esto se hace interpolando una cuadrícula del archivo de altura del geoide con una resolución espacial de 15 minutos de arco.

O si tiene algún nivel profesional GPS tiene Altitud H( msl, altura por encima del nivel medio del mar ) y UNDULATION, la relación entre el geoidy el ellipsoid (m)de la salida de referencia elegida de la tabla interna. puedes obtenerh = H(msl) + undulation

A XYZ por coordenadas cartesianas:

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

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

z = R *sin(lat)
Howie
fuente
¿Cuál es el valor de R?
eych
4
Supongo que es el radio de la esfera, que son 6371 km para la Tierra.
Matthias
5

El software proj.4 proporciona un programa de línea de comandos que puede realizar la conversión, por ejemplo

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

También proporciona una API C . En particular, la función pj_geodetic_to_geocentrichará la conversión sin tener que configurar primero un objeto de proyección.

Brian Hawkins
fuente
5

En python3.x se puede hacer usando:

# 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
fuente
3

Si le interesa obtener coordenadas basadas en un elipsoide en lugar de una esfera, eche un vistazo a http://en.wikipedia.org/wiki/Geodetic_system#From_geodetic_to_ECEF : proporciona las fórmulas y las constantes WGS84 que necesita para la conversión .

Las fórmulas allí también tienen en cuenta la altitud relativa a la superficie del elipsoide de referencia (útil si está obteniendo datos de altitud de un dispositivo GPS).

Stjepan Rajko
fuente
Votar a favor aunque no haya publicado el contenido del enlace aquí.
Mad Physicist
2

¿Por qué implementar algo que ya ha sido implementado y probado?

C #, por ejemplo , tiene NetTopologySuite, que es el puerto .NET de JTS Topology Suite.

Específicamente, tiene una falla grave en su cálculo. La Tierra no es una esfera perfecta, y la aproximación del radio de la Tierra podría no ser suficiente para obtener medidas precisas.

Si en algunos casos es aceptable usar funciones caseras, GIS es un buen ejemplo de un campo en el que se prefiere usar una biblioteca confiable y probada.

Yuval Adam
fuente
1
+1. Usar una biblioteca confiable es más preciso que una función casera y también más fácil .
MarkJ
5
¿Cómo se convierte NetTopologySuite de largo / tardío a cartesion?
vinayan
1
NTS no incluye capacidades de conversión de coordenadas, tal vez necesite Proj.NET projnet.codeplex.com
D_Guidi
6
Ridículo, la respuesta ni siquiera proporciona capacidad de conversión.
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
fuente
¿Podrías elaborarlo? Yo creé una sencilla aplicación que niveles para transformar una sola coordenada usando su enfoque. Sin embargo, siempre falla ya que las dimensiones de la fuente (2) y las dimensiones del objetivo (3) difieren, lo que resulta en una excepciónjava.lang.IllegalArgumentException: dimension must be <= 3
oschrenk
Hmmm ... He mirado a JTS por un rato. Las líneas hasta el nuevo LineString () inclusive se parecen a JTS. Pero no veo las cosas de CRS y Transform en JTS. Entonces: ¿están ahí y los extraño? ¿Estaban allí y se eliminaron en 1.12? O: ¿es una biblioteca diferente?
Dan H
0

Puede hacerlo de esta manera en 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
fuente
¿Qué es el parámetro alt?
baliman
altitud, ¿qué estás haciendo aquí si no sabes cómo funciona el GPS?)
MushyPeas