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

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?

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.

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,

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)

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

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.

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

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.

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

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;


}