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

24

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.

skd
źródło
2
Czy to są dane GRIB ? Jeśli tak, może potrzebujemy tagu grib.
Kirk Kuykendall
@skd wydaje się, że linki ECMWF są nieprawidłowe. Umiesz edytować
gansub,
@gansub Zredagowałem linki. Nie wiem, czy informacje są dokładnie takie same, odkąd minęło sporo czasu, ale wierzę, że nowy link może zapewnić kontekst do wykorzystania w przyszłości.
skd
@skd, kiedy mówisz angle=0.0, masz na myśli łożysko ? Mam plik netcdf ze współrzędnymi obróconego bieguna, ale nie ma wzmianki o żadnym kącie.
FaCoffee,
@ CF84 Właściwie nie jestem pewien. Myślę, że jeśli nie ma wzmianki o kącie, to jest taki sam jak kąt = 0
skd

Odpowiedzi:

24

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(ϑ)) - φ
mkadunc
źródło
1
Pomysł jest dobry, ale niektóre szczegóły wymagają naprawy. l0 = -15, a nie +15. Wszystkie trzy linie w rozwinięciu produktu macierzowego są niepoprawne. ATan2 (lub jego odpowiednik) musi zostać użyty, zmodyfikowany, aby zwracał dowolną rozsądną długość geograficzną, gdy x = y = 0. Zauważ, że ponieważ x ^ 2 + y ^ 2 + z ^ 2 = 1, na końcu dostajesz po prostu lat = Arcsin (z).
whuber
1
Dzięki. Naprawiłem odpowiedź, aby przynajmniej poprawić matematykę. Obroty powinny teraz pasować do opisu w definicji CS, ale trudno jest mieć pewność co do ich znaku bez przykładu (innego niż położenie bieguna południowego).
mkadunc,
Dobra robota! Dziwię się, że ta odpowiedź nie otrzymuje więcej głosów, ponieważ zawiera użyteczny i trudny do znalezienia materiał.
whuber
Naprawdę trudno jest znaleźć materiał, dziękuję bardzo za odpowiedź. W końcu użyłem tego oprogramowania code.zmaw.de/projects/cdo do konwersji z obracanej siatki na zwykłą siatkę. Domyślam się, że najpierw przekształca współrzędne jak w tej odpowiedzi, a następnie interpoluje je, aby dać wyniki w punktach prostokątnej siatki. Chociaż trochę się spóźniam, zostawiam ją jej na przyszłość.
skd
1
@alfe Nie jestem ekspertem od sfer Blocha, ale zasada wygląda bardzo podobnie do tego, co zrobiłem, ale zamiast konwertować do przestrzeni kartezjańskiej z 3 rzeczywistymi współrzędnymi, podpowiedź sugeruje konwersję do przestrzeni z 2 wyobrażonymi współrzędnymi (co oznacza 4 prawdziwe komponenty) i wykonanie tam obrotu. Wywołane przez twój komentarz, połączyłem wszystkie wyrażenia razem i dodałem wynik, w którym pośredni krok kartezjański nie jest już widoczny.
mkadunc,
6

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];
simondk
źródło
Na wypadek gdyby link zniknął, możesz wstawić kod dla przyszłych czytelników. Dzięki.
Michael Stimson
1
Pewnie - wstawiono kod.
simondk
2

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.

Davide
źródło
1

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

Sohrab kolsoomi ayask
źródło
Cordex Data Extractor to program komputerowy dla systemu Windows do wyodrębniania danych z pliku CORDEX NetCDF. Cordex Data Extractor nie potrzebuje pliku pomocy, ponieważ wszystkie procesy zostały wykonane za kodami, a użytkownik po prostu wprowadza daty, współrzędne i nazwę zmiennej. Proszę obejrzeć ten film: youtu.be/RmpZblZPXjI agrimetsoft.com/cordexDataExtractor.aspx
Sohrab kolsoomi ayask
1

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))
użytkownik126158
źródło
0

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.

ujjwalesri
źródło
2
Niestety nie używam żadnego oprogramowania. Są to tylko zestawy danych siatki i zawierają następujące informacje: - Dla pierwszego: ecmwf.int/publications/manuals/d/gribapi/fm92/grib1/detail/… - Dla drugiego: ecmwf.int/publications/ instrukcje / d / gribapi / fm92 / grib1 / detail /…
skd
Ponieważ kąt obrotu wynosi 0, myślę, że proste tłumaczenie powinno wyrównać drugi zestaw danych do pierwszego, tak jak powiedziałeś, dodając 15 do X i 35 do Y
ujjwalesri