Преобразование вручную повернутого широты / долготы в обычные широты / долготы?

24

Сначала я должен уточнить, что у меня нет предыдущего опыта работы в этой области, поэтому я не знаю технической терминологии. Мой вопрос заключается в следующем:

У меня есть два набора данных погоды:

  • Первый имеет регулярную систему координат (я не знаю, имеет ли она конкретное имя), в диапазоне от -90 до 90 и от -180 до 180, а полюса находятся на широтах -90 и 90.

  • Во втором, хотя она должна соответствовать той же области, я заметил что - то другое: широта и долгота не была таким же, как они имеют другую точку отсчета (в описании называются повернутые сетки ). Вместе с парами широта / долгота приходит следующая информация: широта южного полюса: -35,00, длина южного полюса: -15,00, угол: 0,0.

Мне нужно преобразовать вторую пару lon / lat в первую. Это может быть так же просто, как добавить 35 к широтам и 15 к долготам, поскольку угол равен 0, и это кажется простым смещением, но я не уверен.

Редактировать: у меня есть информация о координатах следующая

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

По-видимому, вторая система координат определяется общим вращением сферы

«Один из вариантов для этих параметров:

  • Географическая широта в градусах южного полюса системы координат, например тетап;

  • Географическая долгота в градусах южного полюса системы координат, например, лямбдап;

  • Угол поворота в градусах относительно новой полярной оси (измеренный по часовой стрелке при взгляде с южного на северный полюс) системы координат, предполагая, что новая ось была получена при первом повороте сферы на угол лямбдапа вокруг географической полярной оси и затем поворачивается на (90 + тетап) градусов так, чтобы южный полюс двигался вдоль (ранее повернутого) гринвичского меридиана. "

но все же я не знаю, как преобразовать это в первый.

СКД
источник
2
Так это данные GRIB ? Если так, возможно, нам нужен тег grib.
Кирк Куйкендалл
@skd ссылки на ECMWF не действительны. Вы можете редактировать?
gansub
@gansub Я редактировал ссылки. Я не знаю, является ли информация точно такой же, так как это было давно, но я полагаю, что новая ссылка может обеспечить некоторый контекст для будущей ссылки.
скд
@skd, когда вы говорите angle=0.0, вы имеете в виду подшипник ? У меня есть файл netcdf с координатами повернутого полюса, но ни одного угла не упоминается.
FaCoffee
@ CF84 Я на самом деле не уверен. Я предполагаю, что если нет никакого упоминания об угле, то это то же самое, что и угол = 0
скд

Ответы:

24

Ручное изменение направления вращения должно сработать; где-то должна быть формула для вращения сферических систем координат, но так как я не могу ее найти, вот деривация ( ' отмечает повернутую систему координат; нормальные географические координаты используют простые символы):

Сначала преобразуйте данные во втором наборе данных из сферических (lon ', lat') в (x ', y', z '), используя:

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

Затем используйте две матрицы вращения, чтобы повернуть вторую систему координат так, чтобы она совпадала с первой «нормальной». Мы будем вращать оси координат, чтобы мы могли использовать матрицы вращения оси . Нам нужно поменять знак в матрице на match, чтобы соответствовать смыслу вращения, используемому в определении ECMWF, который, кажется, отличается от стандартного положительного направления.

Поскольку мы отменяем вращение, описанное в определении системы координат, мы сначала поворачиваем на by = - (90 + lat0) = -55 градусов вокруг оси y (вдоль повернутого гринвичского меридиана), а затем на φ = - lon0 = +15 градусов вокруг оси 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')

В расширенном виде это становится:

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

Затем преобразовать обратно в «нормальный» (широта, долгота), используя

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

Если у вас нет atan2, вы можете реализовать его самостоятельно, используя atan (y / x) и исследуя знаки x и y

Убедитесь, что вы преобразовали все углы в радианы, прежде чем использовать тригонометрические функции, иначе вы получите странные результаты; в конце переведите обратно в градусы, если это то, что вы предпочитаете ...

Пример (координаты повернутой сферы ==> стандартные географические координаты):

  • Южный полюс повернутой CS есть (lat0, lon0)

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

  • главный меридиан повернутого CS - это меридиан -15 ° в географическом направлении (повернут на 55 ° к северу)

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

  • симметрия требует, чтобы оба экватора пересекались под углом 90 ° / -90 ° в новой системе координат или под углом 75 ° / -105 ° в географических координатах

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

РЕДАКТИРОВАТЬ: Переписать ответ благодаря очень конструктивному комментарию от Whuber: матрицы и расширения теперь синхронизированы, используя соответствующие знаки для параметров вращения; добавлена ​​ссылка на определение матриц; убрал atan (y / x) из ответа; добавлены примеры конвертации.

РЕДАКТИРОВАТЬ 2: можно получить выражения для того же результата без явного преобразования в декартово пространство. x, y, zВ результате могут быть замещены их соответствующими выражениями, и то же самое можно повторить для x', y'и z'. После применения некоторых тригонометрических тождеств появляются следующие одношаговые выражения:

lat = arcsin(cos(ϑ) sin(lat') - cos(lon') sin(ϑ) cos(lat'))
lon = atan2(sin(lon'), tan(lat') sin(ϑ) + cos(lon') cos(ϑ)) - φ
mkadunc
источник
1
Идея хороша, но некоторые детали нуждаются в исправлении. lon0 = -15, а не +15. Все три строки в разложении матрицы произведения некорректны. ATan2 (или его эквивалент) должен быть использован, модифицированный для возврата любой разумной долготы, когда x = y = 0. Обратите внимание: поскольку x ^ 2 + y ^ 2 + z ^ 2 = 1, в конце вы получите просто lat = Arcsin (z).
whuber
1
Спасибо. Я исправил ответ, чтобы хотя бы сделать математику правильной. Вращения теперь должны соответствовать описанию в определении CS, но трудно быть уверенным в их знаке без примера (кроме положения южного полюса).
mkadunc
Отлично сработано! Я удивлен, что этот ответ не получает больше голосов, потому что он предоставляет полезный и труднодоступный материал.
whuber
Это действительно очень трудно найти материал, большое спасибо за ответ. В итоге я использовал это программное обеспечение code.zmaw.de/projects/cdo для преобразования из повернутой сетки в обычную сетку. Я предполагаю, что сначала он преобразует координаты, как в этом ответе, а затем интерполирует их, чтобы получить результаты в точках прямоугольной сетки. Хотя и немного поздно, я оставлю это ей для дальнейшего использования.
скд
1
@alfe Я не эксперт по сферам Блоха, но принцип выглядит очень похоже на то, что я сделал, но вместо преобразования в декартово пространство с 3 действительными координатами подсказка предлагает преобразовать в пространство с 2 воображаемыми координатами (что означает 4 реальных компонента) и выполнение вращения там. Приведенный в действие вашим комментарием, я собрал все выражения и добавил результат, в котором промежуточный декартовой шаг больше не очевиден.
mkadunc
6

В случае, если кому-то интересно, я поделился сценарием MATLAB на обмене файлами, преобразуя обычный лат / лон в повернутый лат / лон и наоборот: преобразование повернутой сетки

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
источник
Если ссылка не работает, пожалуйста, вставьте код для будущих читателей. Спасибо.
Майкл Стимсон
1
Конечно - код вставлен.
simondk
2

Это преобразование также можно вычислить с помощью программного обеспечения proj (с использованием командной строки или программно), используя то, что proj называет наклонным переводом ( ob_tran), применяемым к преобразованию latlon. Настраиваемые параметры проекции:

  • o_lat_p = широта северного полюса => 35 ° в примере
  • lon_0 = долгота южного полюса => -15 ° в примере
  • o_lon_p = 0

кроме того, -m 57.2957795130823(180 / pi) требуется для того, чтобы учитывать прогнозируемые значения в градусах.

Воспроизведение примеров, предложенных mkadunc, дает тот же результат (обратите внимание, что здесь порядок lon latне (lat,lon)указан, coodinates набираются при стандартном вводе, вывод помечается =>):

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

invprojКоманда используется для преобразования из «спроецированных» (то есть повернутых) координат в географические, а projдля обратного.

Davide
источник
1

Я разработал страницу asp.net для преобразования координат из повернутых в не повернутые на основе доменов CORDEX.

Это основано на вышеуказанных методах. Вы можете свободно использовать его по этой ссылке:

Преобразование вручную повернутого лата / долготы в обычный широту / долготу

Sohrab kolsoomi аяск
источник
Cordex Data Extractor - это программное обеспечение для Windows, предназначенное для извлечения данных из файла CORDEX NetCDF. Cordex Data Extractor не нуждается в файле справки, потому что все процессы были выполнены за кодами, и пользователь просто вводит даты, координаты и имя переменной. Пожалуйста, посмотрите это видео: youtu.be/RmpZblZPXjI agrimetsoft.com/cordexDataExtractor.aspx
Sohrab kolsoomi ayask
1

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

ПИТОН:

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
источник
0

Какое программное обеспечение вы используете? Каждое ГИС-программное обеспечение будет иметь возможность показывать вам текущую информацию о системе координат / проекции. , который может помочь вам в получении названия вашей текущей системы координат.

Кроме того, если вы используете ArcGIS, вы можете использовать инструмент Project для повторного проецирования второго набора данных, импортируя настройки из первого.

ujjwalesri
источник
2
К сожалению, я не использую никакого программного обеспечения. Это просто наборы данных сетки, и они содержат следующую информацию: - Для первого: ecmwf.int/publications/manuals/d/gribapi/fm92/grib1/detail/… - Для второго: ecmwf.int/publications/ руководства / d / gribapi / fm92 / grib1 / detail /…
skd
Поскольку угол поворота равен 0, я думаю, что простой перевод должен выровнять второй набор данных к первому, как вы сказали добавив 15 к X и 35 к Y
ujjwalesri