Mengubah lat / lon yang diputar secara manual menjadi lat / lon biasa?

Pertama saya harus mengklarifikasi saya tidak memiliki pengalaman sebelumnya dengan bidang ini, jadi saya tidak tahu terminologi teknisnya. Pertanyaan saya adalah sebagai berikut:

Saya memiliki dua set data cuaca:

• Yang pertama memiliki sistem koordinat reguler (saya tidak tahu apakah itu memiliki nama tertentu), mulai dari -90 hingga 90 dan -180 hingga 180, dan kutub berada di lintang -90 dan 90.

• Dalam yang kedua, meskipun harus sesuai dengan wilayah yang sama, saya perhatikan sesuatu yang berbeda: lintang dan bujur tidak sama, karena mereka memiliki titik referensi lain (dalam uraian disebut grid yang diputar ). Bersama-sama dengan pasangan lat / lon, muncul informasi berikut: kutub selatan lat: -35,00, kutub selatan lon: -15,00, sudut: 0,0.

Saya perlu mengubah pasangan kedua lon / lat ke yang pertama. Ini bisa sesederhana menambahkan 35 ke garis lintang dan 15 ke garis bujur, karena sudutnya 0 dan sepertinya bergeser sederhana, tapi saya tidak yakin.

Sunting: Informasi yang saya miliki tentang koordinat adalah sebagai berikut

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

Rupanya, sistem koordinat kedua didefinisikan oleh rotasi bola secara umum

"Satu pilihan untuk parameter ini adalah:

• Lintang geografis dalam derajat kutub selatan sistem koordinat, contohnya;

• Garis bujur geografis dalam derajat kutub selatan sistem koordinat, lambdap misalnya;

• Sudut rotasi dalam derajat tentang sumbu kutub baru (diukur searah jarum jam ketika melihat dari selatan ke kutub utara) dari sistem koordinat, dengan asumsi sumbu baru telah diperoleh dengan terlebih dahulu memutar bola melalui derajat lambdap tentang sumbu kutub geografis tentang sumbu kutub geografis , dan kemudian berputar melalui (90 + thetap) derajat sehingga kutub selatan bergerak di sepanjang meridian Greenwich (yang sebelumnya diputar). "

tapi tetap saja saya tidak tahu bagaimana mengubahnya menjadi yang pertama.

Membalik rotasi secara manual harus melakukan trik; harus ada rumus untuk memutar sistem koordinat bola di suatu tempat, tetapi karena saya tidak dapat menemukannya, inilah derivasi ( ' tandai sistem koordinat yang diputar; koordinat geografis normal menggunakan simbol polos):

Pertama-tama, konversikan data dalam dataset kedua dari spherical (lon ', lat') menjadi (x ', y', z ') menggunakan:

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

Kemudian gunakan dua matriks rotasi untuk memutar sistem koordinat kedua sehingga bertepatan dengan yang 'normal' pertama. Kami akan memutar sumbu koordinat, sehingga kami dapat menggunakan matriks rotasi sumbu . Kita perlu membalikkan tanda dalam matriks ϑ untuk mencocokkan pengertian rotasi yang digunakan dalam definisi ECMWF, yang tampaknya berbeda dari arah positif standar.

Karena kita membatalkan rotasi yang dijelaskan dalam definisi sistem koordinat, pertama-tama kita memutar dengan ϑ = - (90 + lat0) = -55 derajat di sekitar sumbu y '(sepanjang meridian Greenwich yang dirotasi) dan kemudian dengan φ = - lon0 = +15 derajat di sekitar sumbu 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')

Diperluas, ini menjadi:

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

Kemudian konversikan kembali ke 'normal' (lat, lon) menggunakan

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

Jika Anda tidak memiliki atan2, Anda dapat menerapkannya sendiri dengan menggunakan atan (y / x) dan memeriksa tanda-tanda x dan y

Pastikan Anda mengubah semua sudut menjadi radian sebelum menggunakan fungsi trigonometri, atau Anda akan mendapatkan hasil yang aneh; konversi kembali ke derajat pada akhirnya jika itu yang Anda inginkan ...

Contoh (koordinat bola yang diputar ==> koordinat geografis standar):

• kutub selatan dari CS yang diputar adalah (lat0, l0)

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

• meridian utama dari CS yang diputar adalah meridian -15 ° secara geografis (diputar 55 ° ke arah utara)

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

• simetri mensyaratkan bahwa kedua ekuator berpotongan pada 90 ° / -90 ° dalam CS baru, atau 75 ° / -105 ° dalam koordinat geografis

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

EDIT: Menulis ulang jawaban berkat komentar yang sangat konstruktif oleh whuber: matriks dan ekspansi sekarang disinkronkan, menggunakan tanda-tanda yang tepat untuk parameter rotasi; menambahkan referensi ke definisi matriks; dihapus atan (y / x) dari jawabannya; menambahkan contoh konversi.

EDIT 2: Dimungkinkan untuk mendapatkan ekspresi untuk hasil yang sama tanpa tranformasi eksplisit ke ruang kartesius. The x, y, zdalam hasil dapat diganti dengan ekspresi yang sesuai mereka, dan sama dapat diulang untuk x', y'dan z'. Setelah menerapkan beberapa identitas trigonometri, ekspresi langkah-tunggal berikut ini muncul:

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

Jika ada yang tertarik, saya telah membagikan skrip MATLAB pada pertukaran file yang mengubah lat / lon biasa menjadi lat / lon yang diputar dan sebaliknya: Rotasi kisi-kisi yang diputar

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

Transformasi ini juga dapat dihitung dengan perangkat lunak proj (baik menggunakan command-line atau secara terprogram) dengan menggunakan apa yang disebut proj sebagai terjemahan miring ( ob_tran) yang diterapkan pada transformasi latlon. Parameter proyeksi yang akan ditetapkan adalah:

• o_lat_p = lintang kutub utara => 35 ° dalam contoh
• lon_0 = bujur kutub selatan => -15 ° dalam contoh
• o_lon_p = 0

tambahan, -m 57.2957795130823(180 / pi) diperlukan untuk mempertimbangkan nilai yang diproyeksikan dalam derajat.

Meniru contoh yang diusulkan oleh mkadunc memberikan hasil yang sama (perhatikan bahwa di sini pesanan lon lattidak (lat,lon), coodinate diketik dalam input standar, output ditandai dengan =>):

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

invprojperintah digunakan untuk mengubah dari "proyeksi" (yaitu diputar) koordinat ke geografis, sedangkan projuntuk melakukan yang sebaliknya.

Saya telah mengembangkan halaman asp.net untuk mengonversi koordinat dari rotated ke non-rotate berdasarkan domain CORDEX.

Itu berdasarkan metode di atas. Anda dapat menggunakannya secara bebas di tautan ini:

Mengubah lat / lon yang diputar secara manual menjadi lat / lon biasa

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

Perangkat lunak apa yang Anda gunakan? Setiap perangkat lunak SIG akan memiliki fasilitas untuk menunjukkan kepada Anda sistem informasi proyeksi / proyeksi saat ini. , yang dapat membantu Anda mendapatkan nama sistem koordinat Anda saat ini.

Selain itu, jika Anda menggunakan ArcGIS, Anda dapat menggunakan alat Proyek untuk memproyeksikan ulang dataset kedua, mengimpor pengaturan dari yang pertama.