Konwersja z longitude \ latitude na współrzędne kartezjańskie

104

Mam kilka punktów współrzędnych ze środkiem Ziemi, podanych jako szerokość i długość geograficzna ( WGS-84 ).

Jak przekonwertować je na współrzędne kartezjańskie (x, y, z), których początek znajduje się w środku Ziemi?

daphshez
źródło
1
Czy udało Ci się przeliczyć długość i szerokość geograficzną WGS-84 na współrzędne kartezjańskie? Mam też elewację. Wypróbowałem tutaj zaakceptowaną odpowiedź, ale nie daje mi ona właściwej odpowiedzi. Porównałem swoje wyniki z tą witryną: apsalin.com/convert-geodetic-to-cartesian.aspx .
Yasmin

Odpowiedzi:

47

Niedawno zrobiłem coś podobnego do tego, używając „Formuły Haversine'a” na danych WGS-84, która jest pochodną „Prawa Hawersinów” z bardzo zadowalającymi wynikami.

Tak, WGS-84 zakłada, że ​​Ziemia jest elipsoidą, ale uważam, że przy zastosowaniu podejścia takiego jak „Formuła Haversine'a” uzyskuje się tylko około 0,5% średniego błędu, co może być dopuszczalną ilością błędu w twoim przypadku. Zawsze będziesz miał pewną ilość błędów, chyba że mówisz o odległości kilku stóp, a nawet wtedy istnieje teoretyczna krzywizna Ziemi ... Jeśli potrzebujesz bardziej sztywnego podejścia zgodnego z WGS-84, sprawdź „Formułę Vincenty'ego”.

Rozumiem, skąd się bierze starblue , ale dobra inżynieria oprogramowania często wymaga kompromisów, więc wszystko zależy od dokładności, jakiej potrzebujesz do tego, co robisz. Na przykład wynik obliczony na podstawie „Formuły odległości Manhattan” w porównaniu z wynikiem z „Formuły odległości” może być lepszy w pewnych sytuacjach, ponieważ jest obliczeniowo tańszy. Pomyśl „który punkt jest najbliżej?” scenariusze, w których nie potrzebujesz dokładnego pomiaru odległości.

Jeśli chodzi o „Formułę Haversine'a”, jest łatwa do wdrożenia i przyjemna, ponieważ wykorzystuje „Trygonometrię sferyczną” zamiast podejścia opartego na „Prawie Cosinusów”, które jest oparte na dwuwymiarowej trygonometrii, dzięki czemu uzyskuje się niezłą równowagę dokładności ponad złożonością.

Dżentelmen o imieniu Chris Veness ma świetną stronę internetową pod adresem http://www.movable-type.co.uk/scripts/latlong.html, która wyjaśnia niektóre interesujące Cię koncepcje i demonstruje różne implementacje programistyczne; Powinno to również odpowiedzieć na pytanie dotyczące konwersji X / Y.

bn.
źródło
1
0,5% błędu - 0,5% czego? W kontekście tego pytania mógłby to być promień Ziemi, więc 0,5% mogłoby wynosić 30 km :)
MarkJ
2
Sprawdziłeś swój link. Cytat 0,5% dotyczy błędu w odległości ortodromy między dwoma punktami, więc nie jest ściśle związany z tym pytaniem. Myślę, że przy przeliczaniu współrzędnych szerokości geograficznej szerokości geograficznej na współrzędne kartezjańskie z początkiem w środku Ziemi błędy wynikające z przyjęcia kulistej ziemi mogą być znaczące. Nie jest jasne, co pytający chce zrobić ze współrzędnymi kartezjańskimi. Albo po prostu wygodniej jest w nich pracować z jakiegoś dziwnego powodu, czy może jest to wymóg eksportu danych? Jeśli to drugie, dokładność byłaby ważna.
MarkJ
132

Oto odpowiedź, którą znalazłem:

Aby uzupełnić definicję, w kartezjańskim układzie współrzędnych:

  • oś x przechodzi przez długość, szerokość (0,0), więc długość geograficzna 0 styka się z równikiem;
  • oś y przechodzi przez (0,90);
  • a oś z przechodzi przez bieguny.

Konwersja to:

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

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

z = R *sin(lat)

Gdzie R jest przybliżonym promieniem Ziemi (np. 6371 km).

Jeśli twoje funkcje trygonometryczne oczekują radianów (co prawdopodobnie robią), będziesz musiał najpierw przeliczyć swoją długość i szerokość geograficzną na radiany. Oczywiście potrzebujesz reprezentacji dziesiętnej, a nie stopni \ minut \ sekund (zobacz np. Tutaj o konwersji).

Wzór na konwersję wsteczną:

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

asin jest oczywiście sinusoidą. przeczytaj o atan2 w Wikipedii . Nie zapomnij przekonwertować z radianów na stopnie.

Ta strona zawiera kod C # do tego (zwróć uwagę, że jest on bardzo różny od formuł), a także wyjaśnienia i ładny diagram, dlaczego jest to poprawne,

daphshez
źródło
19
-1 To jest złe. Zakładasz, że Ziemia jest kulą, podczas gdy WGS-84 przyjmuje elipsoidę.
starblue
43
@starblue: Nie jestem pewien, czy jesteś w stanie oznaczyć daną odpowiedź „właściwą” lub „złą”. Przybliżenie sferyczne (w celu uzyskania współrzędnych x, y, z w stylu ECEF) przy użyciu dostępnych lat / lng (odnoszących się do WGS-84) jest albo „odpowiednie” dla potrzeb oryginalnego plakatu, albo „nieodpowiednie”. Założę się, że w przypadku szacunków odległości i namiaru ta prosta konwersja jest w porządku. Jeśli wystrzeliwuje satelity, może nie. W końcu sam WGS-84 jest „zły” ... w tym sensie, że nie jest idealnym modelem powierzchni Ziemi; wszystkie modele elipsoidalne są przybliżeniami. Szkoda, że ​​OP nie powiedział nam, co próbuje zrobić.
Dan H,
11
@Dan H Pytanie dotyczy WGS-84, a jeśli odpowiesz na coś innego, powinieneś przynajmniej omówić różnice / błąd, czego ta odpowiedź nie zawiera.
starblue
@ daphna-shezaf nie może dokonać konwersji wstecznej ... Zrobiłem również powrót z radianów do stopni, ale wynik nie jest taki sam ...
dziękuję, spędź godzinę na zastanawianiu się, dlaczego to nie działa, okazuje się, że zamieniłem trochę cos (łac.) i sin (łac.)
aeroson
6

Teoria konwersji GPS(WGS84)na współrzędne kartezjańskie https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_geodetic_to_ECEF_coordinates

Oto, czego używam:

  • Długość geograficzna w GPS (WGS84) i współrzędne kartezjańskie są takie same.
  • Szerokość geograficzna musi być przeliczona przez parametry elipsoidy WGS 84, półoś wielka wynosi 6378137 m, a
  • Odwrotność spłaszczenia wynosi 298,257223563.

I załączeniu kod VB I napisał:

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

Zwróć uwagę, że hjest wysokość powyżej WGS 84 ellipsoid.

Zwykle GPSdaje nam wzrost Hpowyżej MSL. MSLWysokość musi być przekształcany do wysokości hwyższej od WGS 84 ellipsoidstosując geopotencjadlna modelu EGM96( Lemoine i wsp, 1998 ).
Odbywa się to poprzez interpolację siatki pliku wysokości geoidy z rozdzielczością przestrzenną 15 minut łuku.

Lub jeśli masz jakiegoś profesjonalistę na poziomie, GPSma Wysokość H( npm, wysokość nad średnim poziomem morza ) i UNDULATIONzwiązek między geoidi ellipsoid (m)a wybranego wyjścia odniesienia z wewnętrznej tabeli. możesz dostaćh = H(msl) + undulation

Do XYZ według współrzędnych kartezjańskich:

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

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

z = R *sin(lat)
Howie
źródło
Jaka jest wartość R?
eych
4
Myślę, że jest to promień kuli, który wynosi 6371 km dla Ziemi.
Matthias
6

W python3.x można to zrobić za pomocą:

# 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
źródło
5

Proj.4 oprogramowanie zapewnia program wiersza poleceń, które można zrobić konwersję, np

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

Zapewnia również C API . W szczególności funkcja pj_geodetic_to_geocentricdokona konwersji bez konieczności uprzedniego konfigurowania obiektu projekcji.

Brian Hawkins
źródło
3

Jeśli zależy Ci na uzyskaniu współrzędnych w oparciu o elipsoidę, a nie na kuli, spójrz na http://en.wikipedia.org/wiki/Geodetic_system#From_geodetic_to_ECEF - zawiera wzory oraz stałe WGS84 potrzebne do konwersji .

Formuły tam również uwzględniają wysokość względem referencyjnej powierzchni elipsoidy (przydatne, jeśli otrzymujesz dane o wysokości z urządzenia GPS).

Stjepan Rajko
źródło
Głosowanie za zgodą, mimo że nie zamieściłeś tutaj treści linku.
Mad Physicist
2

Po co wdrażać coś, co zostało już wdrożone i sprawdzone w testach?

Na przykład język C # ma NetTopologySuite, który jest portem .NET pakietu JTS Topology Suite.

W szczególności masz poważny błąd w obliczeniach. Ziemia nie jest idealną kulą, a przybliżenie promienia Ziemi może jej nie uciąć do precyzyjnych pomiarów.

Jeśli w niektórych przypadkach użycie funkcji homebrew jest dopuszczalne, GIS jest dobrym przykładem dziedziny, w której o wiele lepiej jest używać niezawodnej, sprawdzonej w testach biblioteki.

Yuval Adam
źródło
1
+1. Korzystanie z niezawodnej biblioteki jest dokładniejsze niż funkcja homebrew, a także łatwiejsze .
MarkJ
5
W jaki sposób NetTopologySuite konwertuje z długiego / późnego na kartezja?
vinayan,
1
NTS nie zawiera funkcji konwersji współrzędnych, może potrzebujesz Proj.NET projnet.codeplex.com
D_Guidi
6
Śmieszne, odpowiedź nie zapewnia nawet możliwości konwersji.
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
źródło
Czy byłbyś w stanie to rozwinąć? I stworzył prosty wniosek, że poziomy, aby przekształcić jedno współrzędnych przy użyciu podejścia. Zawsze jednak zawodzi, ponieważ wymiary źródła (2) i wymiary celu (3) różnią się, co skutkuje wyjątkiemjava.lang.IllegalArgumentException: dimension must be <= 3
oschrenk
Hmmm ... Patrzyłem trochę na JTS. Linie do góry, włącznie z nową linią LineString (), wyglądają jak JTS. Ale nie widzę rzeczy CRS i Transform w JTS. Więc: są tam i tęsknię za nimi? Czy były i zostały usunięte w 1.12? Albo: czy to inna biblioteka?
Dan H,
0

Możesz to zrobić w ten sposób na Javie.

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
źródło
co to jest parametr alt?
baliman
wysokość, co tu w ogóle robisz, jeśli nie wiesz, jak działa GPS;)
MushyPeas