Conversão de longitude \ latitude em coordenadas cartesianas

103

Eu tenho alguns pontos de coordenadas centrados na Terra dados como latitude e longitude ( WGS-84 ).

Como posso convertê-los em coordenadas cartesianas (x, y, z) com a origem no centro da Terra?

Daphshez
fonte
1
Você conseguiu converter longitude e latitude WGS-84 em coordenadas cartesianas? Eu também tenho elevação. Tentei a resposta aceita aqui, mas não me deu a resposta certa. Eu comparei meus resultados com este site: apsalin.com/convert-geodetic-to-cartesian.aspx .
Yasmin

Respostas:

46

Recentemente, fiz algo semelhante a isso usando a "Fórmula Haversine" nos dados WGS-84, que é um derivado da "Lei de Haversine" com resultados muito satisfatórios.

Sim, o WGS-84 assume que a Terra é um elipsóide, mas acredito que você obtenha apenas cerca de 0,5% de erro médio usando uma abordagem como a "Fórmula Haversine", que pode ser uma quantidade aceitável de erro no seu caso. Você sempre terá algum erro, a menos que esteja falando sobre uma distância de alguns metros e mesmo assim, teoricamente, há uma curvatura da Terra ... Se você precisar de uma abordagem compatível com WGS-84 mais rígida, verifique a "Fórmula de Vincenty".

Eu entendo de onde vem o starblue , mas uma boa engenharia de software geralmente trata de compensações, então tudo depende da precisão que você precisa para o que está fazendo. Por exemplo, o resultado calculado da "Fórmula da distância de Manhattan" em comparação com o resultado da "Fórmula da distância" pode ser melhor para certas situações, pois é computacionalmente mais barato. Pense "qual é o ponto mais próximo?" cenários onde você não precisa de uma medição de distância precisa.

Em relação à "Fórmula Haversine", é fácil de implementar e é agradável porque usa "Trigonometria Esférica" ​​em vez de uma abordagem baseada na "Lei dos Cossenos" que é baseada na trigonometria bidimensional, portanto, você obtém um bom equilíbrio de precisão sobre a complexidade.

Um cavalheiro chamado Chris Veness tem um ótimo site em http://www.movable-type.co.uk/scripts/latlong.html que explica alguns dos conceitos nos quais você está interessado e demonstra várias implementações programáticas; isso deve responder também à sua pergunta de conversão x / y.

bn.
fonte
1
0,5% de erro - 0,5% de quê? No contexto desta questão, poderia ser o raio da Terra, então 0,5% poderia ser 30 km :)
MarkJ
2
Verifiquei seu link. A cotação de 0,5% é para o erro na distância do grande círculo entre dois pontos, não sendo estritamente relevante para esta questão. Eu pensaria que ao converter lat-long em coordenadas cartesianas com a origem no centro da Terra, os erros de presumir uma Terra esférica poderiam ser significativos. Não está claro o que o questionador deseja fazer com as coordenadas cartesianas. Ou é apenas mais conveniente trabalhar neles por algum motivo bizarro ou talvez seja algum requisito para exportação de dados? Neste último caso, a precisão seria importante.
MarkJ
130

Aqui está a resposta que encontrei:

Só para completar a definição, no sistema de coordenadas cartesianas:

  • o eixo x passa por longo, lat (0,0), então a longitude 0 encontra o equador;
  • o eixo y passa por (0,90);
  • e o eixo z passa pelos pólos.

A conversão é:

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

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

z = R *sin(lat)

Onde R está o raio aproximado da terra (por exemplo, 6371 km).

Se suas funções trigonométricas esperam radianos (o que provavelmente acontece), você precisará primeiro converter sua longitude e latitude em radianos. Você obviamente precisa de uma representação decimal, não graus \ minutos \ segundos (veja por exemplo, aqui sobre conversão).

A fórmula para conversão reversa:

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

asin é, naturalmente, arco seno. leia sobre atan2 na wikipedia . Não se esqueça de converter de volta de radianos para graus.

Esta página fornece o código c # para isso (observe que é muito diferente das fórmulas), e também algumas explicações e um belo diagrama de por que isso está correto,

Daphshez
fonte
17
-1 Isso está errado. Você está assumindo que a Terra é uma esfera, enquanto o WGS-84 assume um elipsóide.
starblue
42
@starblue: Não tenho certeza se você está em posição de rotular a resposta dada como "certa" ou "errada". A aproximação esférica (para obter as coordenadas x, y, z do estilo ECEF) usando lat / lngs disponíveis (que são referenciados a WGS-84) é "adequada" para as necessidades do autor original ou "inadequada". Para estimativas de distância e direção, aposto que essa conversão simples é adequada. Se ele está lançando satélites, talvez não. Afinal, o próprio WGS-84 está "errado" ... por não ser um modelo perfeito da superfície da Terra; todos os modelos elipsoidais são aproximações. Pena que o OP não nos disse o que estava tentando fazer.
Dan H
11
@Dan H A questão pede WGS-84, e se você responder alguma outra coisa, você deve pelo menos discutir as diferenças / erros, o que esta resposta não faz.
starblue
@ daphna-shezaf não pode fazer uma conversão reversa ... Também fiz a conversão reversa de radianos para graus, mas o resultado não é o mesmo ...
obrigado, passe horas descobrindo por que isso não funciona, acabei descobrindo que troquei alguns cos (lat) e sin (lat)
aeroson
6

Teoria para converter GPS(WGS84)em coordenadas cartesianas https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_geodetic_to_ECEF_coordinates

O seguinte é o que estou usando:

  • A longitude no GPS (WGS84) e as coordenadas cartesianas são as mesmas.
  • A latitude precisa ser convertida pelos parâmetros do semi-eixo maior do elipsóide WGS 84 é 6378137 m, e
  • O recíproco de achatamento é 298,257223563.

Anexei um código VB que escrevi:

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

Observe que a haltitude está acima de WGS 84 ellipsoid.

Normalmente GPSnos dará altura Hacima MSL. A MSLaltura deve ser convertida para a altura hacima do WGS 84 ellipsoidusando o modelo geopotencialEGM96 ( Lemoine et al, 1998 ).
Isso é feito interpolando uma grade do arquivo de altura do geóide com uma resolução espacial de 15 minutos de arco.

Ou se você tiver algum nível profissional GPS tem Altitude H( msl, altura acima do nível médio do mar ) e UNDULATION, a relação entre o geoide o ellipsoid (m)da saída do datum escolhido da tabela interna. você pode terh = H(msl) + undulation

Para XYZ por coordenadas cartesianas:

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

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

z = R *sin(lat)
Howie
fonte
Qual é o valor de R?
eych
4
Eu acho que é o raio da esfera, que é 6.371 km para a Terra.
Matthias
5

O software proj.4 fornece um programa de linha de comando que pode fazer a conversão, por exemplo

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

Ele também fornece um API C . Em particular, a função pj_geodetic_to_geocentricfará a conversão sem ter que configurar um objeto de projeção primeiro.

Brian Hawkins
fonte
5

Em python3.x, isso pode ser feito 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
fonte
3

Se você se preocupa em obter coordenadas com base em um elipsóide em vez de uma esfera, dê uma olhada em http://en.wikipedia.org/wiki/Geodetic_system#From_geodetic_to_ECEF - ele fornece as fórmulas, bem como as constantes WGS84 de que você precisa para a conversão .

As fórmulas ali também levam em consideração a altitude relativa à superfície do elipsóide de referência (útil se você estiver obtendo dados de altitude de um dispositivo GPS).

Stjepan Rajko
fonte
Votar positivamente mesmo que você não tenha postado o conteúdo do link aqui.
Mad Physicist
2

Por que implementar algo que já foi implementado e testado?

C #, por exemplo, tem o NetTopologySuite, que é a porta .NET do JTS Topology Suite.

Especificamente, você tem uma falha grave em seu cálculo. A Terra não é uma esfera perfeita, e a aproximação do raio da Terra pode não cortá-la para medições precisas.

Se em alguns casos for aceitável usar funções homebrew, o GIS é um bom exemplo de um campo no qual é preferível usar uma biblioteca confiável e comprovada por teste.

Yuval Adam
fonte
1
+1. Usar uma biblioteca confiável é mais preciso do que uma função homebrew e também mais fácil .
MarkJ
5
Como o NetTopologySuite converte de longo / atrasado em cartesion?
vinayan de
1
NTS não inclui recursos de conversão de coordenadas, talvez você precise do Proj.NET projnet.codeplex.com
D_Guidi
6
Ridículo, a resposta nem oferece capacidade de conversão.
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
fonte
Você seria capaz de elaborar? Eu criei um aplicativo simples com camadas para transformar uma única coordenada usando sua abordagem. No entanto, ele sempre falha, pois as dimensões da fonte (2) e as dimensões do alvo (3) diferem, resultando em uma exceçãojava.lang.IllegalArgumentException: dimension must be <= 3
oschrenk
Hmmm ... Eu olhei para o JTS um pouco. As linhas até e incluindo o novo LineString () parecem JTS. Mas eu não vejo o CRS e as coisas de Transform no JTS. Então: eles estão aí e eu estou sentindo falta deles? Estiveram lá e removidos em 1.12? Ou: é uma biblioteca diferente?
Dan H
0

Você pode fazer isso em 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
fonte
o que é o parâmetro alt?
baliman
altitude, o que você está fazendo aqui se não sabe como funciona o GPS;)
MushyPeas