Czy ręcznie przekształcasz obrócone lat / lon w zwykłe lat / lon?

Najpierw powinienem wyjaśnić, że nie mam wcześniejszego doświadczenia w tej dziedzinie, więc nie znam terminologii technicznej. Moje pytanie jest następujące:

Mam dwa zestawy danych pogodowych:

• Pierwszy ma regularny układ współrzędnych (nie wiem, czy ma określoną nazwę), w zakresie od -90 do 90 i -180 do 180, a bieguny znajdują się na szerokościach -90 i 90.

• W drugim, chociaż powinien on odpowiadać temu samemu regionowi, zauważyłem coś innego: szerokość i długość geograficzna nie były takie same, ponieważ mają inny punkt odniesienia (w opisie nazywany jest siatką obróconą ). Wraz z parami lat / lon otrzymujemy następujące informacje: południowy biegun lat: -35,00, południowy biegun lon: -15,00, kąt: 0,0.

Muszę przekształcić drugą parę lon / lat na pierwszą. Może to być tak proste, jak dodanie 35 do szerokości i 15 do długości, ponieważ kąt wynosi 0 i wydaje się prostą zmianą, ale nie jestem pewien.

Edycja: Mam informacje o współrzędnych są następujące

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

Najwyraźniej drugi układ współrzędnych jest określony przez ogólny obrót kuli

„Jeden z tych parametrów to:

• Szerokość geograficzna w stopniach południowego bieguna układu współrzędnych, na przykład dotknięcie;

• Długość geograficzna w stopniach południowego bieguna układu współrzędnych, na przykład lambdap;

• Kąt obrotu w stopniach wokół nowej osi biegunowej (mierzony zgodnie z ruchem wskazówek zegara, patrząc od bieguna południowego do północnego) układu współrzędnych, przy założeniu, że nową oś uzyskano najpierw poprzez obrócenie kuli o stopnie lambda wokół geograficznej osi biegunowej , a następnie obracanie o (90 + dotknięcie) stopni, tak aby biegun południowy poruszał się wzdłuż (wcześniej obróconego) południka Greenwich. ”

ale wciąż nie wiem, jak przekonwertować to na pierwsze.

Ręczne odwrócenie obrotu powinno załatwić sprawę; powinien gdzieś istnieć wzór na obracanie sferycznych układów współrzędnych, ale ponieważ nie mogę go znaleźć, oto pochodna ( „ oznacza obrócony układ współrzędnych; normalne współrzędne geograficzne używają prostych symboli):

Najpierw przekonwertuj dane w drugim zestawie danych ze sferycznego (lon ', lat') na (x ', y', z '), używając:

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

Następnie użyj dwóch macierzy obrotu, aby obrócić drugi układ współrzędnych, tak aby pokrywał się z pierwszym „normalnym”. Będziemy obracać osie współrzędnych, abyśmy mogli użyć macierzy obrotu osi . Musimy odwrócić znak w macierzy ϑ, aby pasował do sensu obrotu stosowanego w definicji ECMWF, który wydaje się różnić od standardowego kierunku dodatniego.

Ponieważ cofamy obrót opisany w definicji układu współrzędnych, najpierw obracamy o ϑ = - (90 + lat0) = -55 stopni wokół osi y '(wzdłuż obróconego południka Greenwich), a następnie o φ = - lon0 = +15 stopni wokół osi Z):

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

Po rozwinięciu staje się:

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

Następnie przekonwertuj z powrotem na „normalny” (lat, lon) za pomocą

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

Jeśli nie masz atan2, możesz go zaimplementować samodzielnie, używając atan (t / x) i sprawdzając znaki xiy

Przed użyciem funkcji trygonometrycznych przekonwertuj wszystkie kąty na radiany, aby uzyskać dziwne wyniki; przelicz z powrotem na stopnie, jeśli tak wolisz ...

Przykład (obrócone współrzędne kuli ==> standardowe współrzędne geograficzne):

• biegun południowy obróconego CS to (lat0, lon0)

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

• głównym południkiem obróconego CS jest południk -15 ° geograficznie (obrócony o 55 ° w kierunku północnym)

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

• symetria wymaga, aby oba równiki przecinały się pod kątem 90 ° / -90 ° w nowym CS lub 75 ° / -105 ° we współrzędnych geograficznych

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

EDYCJA: Przepisałem odpowiedź dzięki bardzo konstruktywnemu komentarzowi whuber'a: macierze i rozwinięcie są teraz zsynchronizowane, używając odpowiednich znaków dla parametrów obrotu; dodano odniesienie do definicji macierzy; usunięto atan (t / x) z odpowiedzi; dodane przykłady konwersji.

EDYCJA 2: Możliwe jest uzyskanie wyrażeń dla tego samego wyniku bez wyraźnej transformacji w przestrzeń kartezjańską. x, y, zW rezultacie może być podstawiona przez odpowiadające im wyrażenia, a tym samym mogą być powtarzane x', y'a z'. Po zastosowaniu niektórych tożsamości trygonometrycznych pojawiają się następujące wyrażenia jednoetapowe:

lat = arcsin(cos(ϑ) sin(lat') - cos(lon') sin(ϑ) cos(lat'))
lon = atan2(sin(lon'), tan(lat') sin(ϑ) + cos(lon') cos(ϑ)) - φ

Jeśli ktoś jest zainteresowany, udostępniłem skrypt MATLAB na wymianie plików, przekształcając zwykłe lat / lon w obrócone lat / lon i odwrotnie: Obrócona transformacja siatki

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

Transformację tę można również obliczyć za pomocą oprogramowania proj (przy użyciu wiersza polecenia lub programowo), wykorzystując to, co proj nazywa tłumaczeniem ukośnym ( ob_tran) zastosowanym do transformacji latlon. Parametry projekcji do ustawienia to:

• o_lat_p = szerokość geograficzna bieguna północnego => 35 ° w przykładzie
• lon_0 = długość bieguna południowego => -15 ° w przykładzie
• o_lon_p = 0

dodatkowo -m 57.2957795130823wymagane jest (180 / pi), aby uwzględnić przewidywane wartości w stopniach.

Replikacja przykładów zaproponowanych przez mkadunc daje ten sam wynik (zauważ, że w tym przypadku kolejność lon latnie występuje (lat,lon), znaki standardowe są wpisywane na standardowym wejściu, dane wyjściowe są oznaczone przez =>):

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

invprojpolecenie służy do konwersji współrzędnych „rzutowanych” (tj. obróconych) na geograficzne, podczas gdy projrobi się odwrotnie.

Opracowałem stronę asp.net do konwersji współrzędnych z obróconych na nieobrócone na podstawie domen CORDEX.

Opiera się na powyższych metodach. Możesz go swobodnie używać w tym linku:

Ręczne przekształcanie obróconych lat / lon na zwykłe lat / lon

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

PYTON:

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

Z jakiego oprogramowania korzystasz? Każde oprogramowanie GIS będzie miało funkcję pokazującą aktualne informacje o systemie / projekcji. , co może pomóc w uzyskaniu nazwy bieżącego układu współrzędnych.

Ponadto, jeśli korzystasz z ArcGIS, możesz użyć narzędzia Projekt , aby ponownie rzutować drugi zestaw danych, importując ustawienia z pierwszego.

Moduł kartopii dla pytona ma funkcję obracania biegunów. Rozumie, jak się wydaje transformacje.

https://scitools.org.uk/cartopy/docs/latest/gallery/rotated_pole.html#sphx-glr-gallery-rotated-pole-py

Chyba może być pomocna.