Manuelles Umwandeln von gedrehtem Lat / Lon in reguläres Lat / Lon?

23

Zuerst sollte ich klarstellen, dass ich keine Vorkenntnisse in diesem Bereich habe, daher kenne ich die technische Terminologie nicht. Meine Frage lautet wie folgt:

Ich habe zwei Wetterdatensätze:

  • Das erste hat das reguläre Koordinatensystem (ich weiß nicht, ob es einen bestimmten Namen hat), das von -90 bis 90 und -180 bis 180 reicht, und die Pole liegen bei den Breiten -90 und 90.

  • Im zweiten Fall habe ich etwas anderes bemerkt, obwohl es der gleichen Region entsprechen sollte: Breite und Länge waren nicht gleich, da sie einen anderen Bezugspunkt haben (in der Beschreibung als gedrehtes Gitter bezeichnet ). Zusammen mit den Lat / Lon-Paaren erhalten Sie folgende Informationen: Südpol Lat: -35,00, Südpol Lon: -15,00, Winkel: 0,0.

Ich muss das zweite Paar von lon / lat in das erste verwandeln. Es könnte so einfach sein, 35 zu den Breitengraden und 15 zu den Längengraden hinzuzufügen, da der Winkel 0 ist und es eine einfache Verschiebung zu sein scheint, aber ich bin nicht sicher.

Bearbeiten: Die Informationen, die ich über die Koordinaten habe, sind die folgenden

http://rda.ucar.edu/docs/formats/grib/gribdoc/llgrid.html

Offensichtlich ist das zweite Koordinatensystem durch eine allgemeine Drehung der Kugel definiert

Eine Wahl für diese Parameter ist:

  • Die geografische Breite in Grad des Südpols des Koordinatensystems, z.

  • Die geografische Länge in Grad des Südpols des Koordinatensystems, z. B. Lambdap;

  • Der Drehwinkel in Grad um die neue Polarachse (gemessen im Uhrzeigersinn vom Südpol zum Nordpol) des Koordinatensystems, wobei angenommen wird, dass die neue Achse erhalten wurde, indem die Kugel zunächst um Lambdap-Grad um die geografische Polarachse gedreht wurde und dann um (90 + thetap) Grad drehen, so dass sich der Südpol entlang des (zuvor gedrehten) Greenwich-Meridians bewegte. "

aber ich weiß immer noch nicht, wie ich das in das erste umwandeln soll.

skd
quelle
2
Sind das also GRIB- Daten? Wenn ja, brauchen wir vielleicht ein Grib-Tag.
Kirk Kuykendall
@skd die ECMWF-Links scheinen nicht gültig zu sein. Kannst du bearbeiten
Gansub
@gansub Ich habe die Links bearbeitet. Ich weiß nicht, ob die Informationen genau die gleichen sind, da es lange her ist, aber ich glaube, dass der neue Link einen Kontext für zukünftige Verweise bieten kann.
Skd
@skd wenn du sagst angle=0.0, meinst du das Lager ? Ich habe eine netcdf-Datei mit den gedrehten Polkoordinaten, aber es wird kein Winkel erwähnt.
FaCoffee
@ CF84 Ich bin mir eigentlich nicht sicher. Ich denke, wenn es keine Erwähnung des Winkels gibt, dann ist es dasselbe wie angle = 0
skd

Antworten:

24

Manuelles Umkehren der Rotation sollte den Trick tun. Es sollte irgendwo eine Formel für rotierende sphärische Koordinatensysteme geben, aber da ich sie nicht finden kann, ist hier die Ableitung ( ' markiert das gedrehte Koordinatensystem; normale geografische Koordinaten verwenden einfache Symbole):

Konvertieren Sie zuerst die Daten im zweiten Datensatz von sphärisch (lon ', lat') nach (x ', y', z ') mit:

x' = cos(lon')*cos(lat')
y' = sin(lon')*cos(lat')
z' = sin(lat')

Verwenden Sie dann zwei Rotationsmatrizen, um das zweite Koordinatensystem so zu drehen, dass es mit dem ersten "normalen" zusammenfällt. Wir werden die Koordinatenachsen drehen, damit wir die Achsendrehungsmatrizen verwenden können . Wir müssen das Vorzeichen in der ϑ-Matrix umkehren, um es mit dem in der ECMWF-Definition verwendeten Rotationssinn abzustimmen, der sich von der standardmäßigen positiven Richtung zu unterscheiden scheint.

Da wir die in der Definition des Koordinatensystems beschriebene Drehung rückgängig machen, drehen wir zuerst um ϑ = - (90 + lat0) = -55 Grad um die y'-Achse (entlang des gedrehten Greenwich-Meridians) und dann um φ = - lon0 = +15 Grad um die z-Achse):

x   ( cos(φ), sin(φ), 0) (  cos(ϑ), 0, sin(ϑ)) (x')
y = (-sin(φ), cos(φ), 0).(  0     , 1, 0     ).(y')
z   ( 0     , 0     , 1) ( -sin(ϑ), 0, cos(ϑ)) (z')

Erweitert wird dies zu:

x = cos(ϑ) cos(φ) x' + sin(φ) y' + sin(ϑ) cos(φ) z'
y = -cos(ϑ) sin(φ) x' + cos(φ) y' - sin(ϑ) sin(φ) z'
z = -sin(ϑ) x' + cos(ϑ) z'

Dann konvertiere zurück zu 'normal' (lat, lon) mit

lat = arcsin(z)
lon = atan2(y, x)

Wenn Sie nicht über atan2 verfügen, können Sie es selbst implementieren, indem Sie atan (y / x) verwenden und die Vorzeichen von x und y untersuchen

Vergewissern Sie sich, dass Sie alle Winkel in Bogenmaß konvertieren, bevor Sie die trigonometrischen Funktionen verwenden, sonst erhalten Sie seltsame Ergebnisse. wandle es am Ende wieder in Grad um, wenn es dir lieber ist ...

Beispiel (gedrehte Kugelkoordinaten ==> geografische Standardkoordinaten):

  • Südpol des gedrehten CS ist (lat0, lon0)

    (-90 °, *) ==> (-35 °, -15 °)

  • Hauptmeridian des gedrehten CS ist der geografische Meridian von -15 ° (55 ° nach Norden gedreht)

    (0 °, 0 °) ==> (55 °, -15 °)

  • Symmetrie erfordert, dass sich beide Äquatoren im neuen CS bei 90 ° / -90 ° oder in geografischen Koordinaten bei 75 ° / -105 ° schneiden

    (0 °, 90 °) ==> (0 °, 75 °)
    (0 °, -90 °) ==> (0 °, -105 °)

BEARBEITEN: Die Antwort wurde dank eines sehr konstruktiven Kommentars von whuber umgeschrieben: Die Matrizen und die Erweiterung sind jetzt synchron, wobei die richtigen Vorzeichen für die Rotationsparameter verwendet werden. Verweis auf die Definition der Matrizen hinzugefügt; atan (y / x) aus der Antwort entfernt; Konvertierungsbeispiele hinzugefügt.

EDIT 2: Es ist möglich, Ausdrücke für dasselbe Ergebnis ohne explizite Transformation in den kartesischen Raum abzuleiten. Die x, y, zim Ergebnis kann mit ihrem entsprechenden Ausdrücke ersetzt werden, und das gleiche kann für die wiederholt werden x', y'und z'. Nach dem Anwenden einiger trigonometrischer Identitäten ergeben sich die folgenden Einzelschrittausdrücke:

lat = arcsin(cos(ϑ) sin(lat') - cos(lon') sin(ϑ) cos(lat'))
lon = atan2(sin(lon'), tan(lat') sin(ϑ) + cos(lon') cos(ϑ)) - φ
mkadunc
quelle
1
Die Idee ist gut, aber einige Details müssen noch korrigiert werden. lon0 = -15, nicht +15. Alle drei Zeilen in der Erweiterung des Matrixprodukts sind falsch. ATan2 (oder sein Äquivalent) muss verwendet werden und so modifiziert werden, dass eine angemessene Länge zurückgegeben wird, wenn x = y = 0 ist. Beachten Sie, dass Sie am Ende einfach lat = Arcsin (z) erhalten, da x ^ 2 + y ^ 2 + z ^ 2 = 1.
whuber
1
Vielen Dank. Ich habe die Antwort korrigiert, um zumindest die Mathematik zu korrigieren. Die Rotationen sollten jetzt mit der Beschreibung in der CS-Definition übereinstimmen, es ist jedoch schwierig, ohne ein Beispiel (außer der Position des Südpols) sicher zu sein, ob sie ein Vorzeichen haben.
mkadunc
Gut gemacht! Ich bin überrascht, dass diese Antwort nicht mehr Stimmen erhält, da sie nützliches und schwer zu findendes Material enthält.
whuber
Dies ist in der Tat sehr schwer zu finden, vielen Dank für die Antwort. Schließlich habe ich diese Software code.zmaw.de/projects/cdo zum Konvertieren von einem gedrehten Raster in ein reguläres Raster verwendet. Ich vermute, dass es zuerst die Koordinaten wie in dieser Antwort transformiert und sie dann interpoliert, um die Ergebnisse an den Punkten eines rechteckigen Gitters zu erhalten. Obwohl es ein bisschen spät ist, lasse ich dies für zukünftige Referenz.
Skd
1
@alfe Ich bin kein Experte für Blochkugeln, aber das Prinzip ähnelt dem, was ich gemacht habe, aber anstatt in einen kartesischen Raum mit 3 reellen Koordinaten zu konvertieren, schlägt der Hinweis vor, in einen Raum mit 2 imaginären Koordinaten zu konvertieren (was bedeutet) 4 reale Komponenten) und dort die Rotation durchführen. Ausgelöst von Ihrem Kommentar habe ich alle Ausdrücke zusammengefasst und ein Ergebnis hinzugefügt, bei dem der kartesische Zwischenschritt nicht mehr erkennbar ist.
Mkadunc
6

Falls sich jemand dafür interessiert, habe ich ein MATLAB-Skript für den Dateiaustausch freigegeben, das reguläre Lat / Lon-Transformationen in gedrehte Lat / Lon-Transformationen und umgekehrt: gedrehte Gittertransformationen durchführt

function [grid_out] = rotated_grid_transform(grid_in, option, SP_coor)

lon = grid_in(:,1);
lat = grid_in(:,2);

lon = (lon*pi)/180; % Convert degrees to radians
lat = (lat*pi)/180;

SP_lon = SP_coor(1);
SP_lat = SP_coor(2);

theta = 90+SP_lat; % Rotation around y-axis
phi = SP_lon; % Rotation around z-axis

phi = (phi*pi)/180; % Convert degrees to radians
theta = (theta*pi)/180;

x = cos(lon).*cos(lat); % Convert from spherical to cartesian coordinates
y = sin(lon).*cos(lat);
z = sin(lat);

if option == 1 % Regular -> Rotated

    x_new = cos(theta).*cos(phi).*x + cos(theta).*sin(phi).*y + sin(theta).*z;
    y_new = -sin(phi).*x + cos(phi).*y;
    z_new = -sin(theta).*cos(phi).*x - sin(theta).*sin(phi).*y + cos(theta).*z;

elseif option == 2 % Rotated -> Regular

    phi = -phi;
    theta = -theta;

    x_new = cos(theta).*cos(phi).*x + sin(phi).*y + sin(theta).*cos(phi).*z;
    y_new = -cos(theta).*sin(phi).*x + cos(phi).*y - sin(theta).*sin(phi).*z;
    z_new = -sin(theta).*x + cos(theta).*z;

end

lon_new = atan2(y_new,x_new); % Convert cartesian back to spherical coordinates
lat_new = asin(z_new);

lon_new = (lon_new*180)/pi; % Convert radians back to degrees
lat_new = (lat_new*180)/pi;

grid_out = [lon_new lat_new];
simondk
quelle
Falls der Link nicht mehr funktioniert, geben Sie bitte den Code für zukünftige Leser ein. Vielen Dank.
Michael Stimson
1
Sicher - Code eingefügt.
simondk
2

Diese Transformation kann auch mit proj- Software (entweder über die Befehlszeile oder programmgesteuert) berechnet werden, indem das verwendet wird, was proj als schräge Übersetzung ( ob_tran) bezeichnet und auf eine Latlon-Transformation angewendet wird. Die einzustellenden Projektionsparameter sind:

  • o_lat_p = Nordpol Breitengrad => 35 ° im Beispiel
  • lon_0 = Südpol-Längengrad => -15 ° im Beispiel
  • o_lon_p = 0

Zusätzlich ist -m 57.2957795130823(180 / pi) erforderlich, um die projizierten Werte in Grad zu berücksichtigen.

Wenn Sie die von mkadunc vorgeschlagenen Beispiele wiederholen, erhalten Sie dasselbe Ergebnis (beachten Sie, dass hier die Reihenfolge lon latnicht (lat,lon)zutrifft, Coodinaten in der Standardeingabe eingegeben werden und die Ausgabe durch gekennzeichnet ist =>):

invproj -f "=> %.6f" -m 57.2957795130823 +proj=ob_tran +o_proj=latlon +o_lon_p=0 +o_lat_p=35 +lon_0=-15
0 -90
=> -15.000000   => -35.000000
40 -90
=> -15.000000   => -35.000000
0 0
=> -15.000000   => 55.000000
90 0
=> 75.000000    => -0.000000
-90 0
=> -105.000000  => -0.000000

invprojDer Befehl wird zum Konvertieren von "projizierten" (dh gedrehten) Koordinaten in geografische Koordinaten verwendet, während projdas Gegenteil der Fall ist.

Davide
quelle
1

Ich habe eine asp.net-Seite zum Konvertieren von Koordinaten von gedreht auf nicht gedreht basierend auf CORDEX-Domänen entwickelt.

Es basiert auf obigen Methoden. Sie können es in diesem Link frei verwenden:

Manuelles Umwandeln von gedrehtem Lat / Lon in reguläres Lat / Lon

Sohrab kolsoomi ayask
quelle
Cordex Data Extractor ist eine Windows-Desktop-Software zum Extrahieren von Daten aus der CORDEX NetCDF-Datei. Cordex Data Extractor benötigt keine Hilfedatei, da alle Prozesse hinter Codes ausgeführt wurden und der Benutzer nur Datums-, Koordinaten- und Variablennamen eingibt. Bitte sehen Sie dieses Video: youtu.be/RmpZblZPXjI agrimetsoft.com/cordexDataExtractor.aspx
Sohrab kolsoomi ayask
1

https://www.mathworks.com/matlabcentral/fileexchange/43435-rotated-grid-transform

PYTHON:

from math import *

def rotated_grid_transform(grid_in, option, SP_coor):
    lon = grid_in[0]
    lat = grid_in[1];

    lon = (lon*pi)/180; # Convert degrees to radians
    lat = (lat*pi)/180;

    SP_lon = SP_coor[0];
    SP_lat = SP_coor[1];

    theta = 90+SP_lat; # Rotation around y-axis
    phi = SP_lon; # Rotation around z-axis

    theta = (theta*pi)/180;
    phi = (phi*pi)/180; # Convert degrees to radians

    x = cos(lon)*cos(lat); # Convert from spherical to cartesian coordinates
    y = sin(lon)*cos(lat);
    z = sin(lat);

    if option == 1: # Regular -> Rotated

        x_new = cos(theta)*cos(phi)*x + cos(theta)*sin(phi)*y + sin(theta)*z;
        y_new = -sin(phi)*x + cos(phi)*y;
        z_new = -sin(theta)*cos(phi)*x - sin(theta)*sin(phi)*y + cos(theta)*z;

    else:  # Rotated -> Regular

        phi = -phi;
        theta = -theta;

        x_new = cos(theta)*cos(phi)*x + sin(phi)*y + sin(theta)*cos(phi)*z;
        y_new = -cos(theta)*sin(phi)*x + cos(phi)*y - sin(theta)*sin(phi)*z;
        z_new = -sin(theta)*x + cos(theta)*z;



    lon_new = atan2(y_new,x_new); # Convert cartesian back to spherical coordinates
    lat_new = asin(z_new);

    lon_new = (lon_new*180)/pi; # Convert radians back to degrees
    lat_new = (lat_new*180)/pi;

    print lon_new,lat_new;

rotated_grid_transform((0,0), 1, (0,30))
user126158
quelle
0

Welche Software verwenden Sie? Jede GIS-Software kann Ihnen die aktuellen Informationen zum Koordinatensystem / zur Projektion anzeigen. , mit dessen Hilfe Sie den Namen Ihres aktuellen Koordinatensystems ermitteln können.

Wenn Sie ArcGIS verwenden, können Sie das zweite Dataset mit dem Tool " Projekt" erneut projizieren und dabei Einstellungen aus dem ersten importieren.

ujjwalesri
quelle
2
Leider benutze ich keine Software. Das sind nur Raster Datensätze und sie kommen mit den folgenden Informationen: - Für die erste: ecmwf.int/publications/manuals/d/gribapi/fm92/grib1/detail/... - Für die zweite: ecmwf.int/publications/ manual / d / gribapi / fm92 / grib1 / detail /…
skd
Da der Drehwinkel 0 ist, glaube ich , eine einfache Übersetzung die zweite Datensatz zu dem ersten ausgerichtet ist, wie gesagt , du 15 bis X 35 und Y an der Zugabe
ujjwalesri