Transformando manualmente lat / lon girado em lat / lon regular?

23

Primeiro, devo esclarecer que não tenho experiência anterior com o campo, portanto não conheço a terminologia técnica. Minha pergunta é a seguinte:

Eu tenho dois conjuntos de dados climáticos:

  • O primeiro tem o sistema de coordenadas regular (não sei se tem um nome específico), variando de -90 a 90 e -180 a 180, e os polos estão nas latitudes -90 e 90.

  • No segundo, embora devesse corresponder à mesma região, notei algo diferente: latitude e longitude não eram iguais, pois têm outro ponto de referência (na descrição é chamado de grade rotativa ). Juntamente com os pares lat / lon, vêm as seguintes informações: pólo sul lat: -35,00, pólo sul lon: -15,00, ângulo: 0,0.

Eu preciso transformar o segundo par de lon / lat no primeiro. Pode ser tão simples quanto adicionar 35 às latitudes e 15 às longitudes, já que o ângulo é 0 e parece uma mudança simples, mas não tenho certeza.

Editar: As informações que tenho sobre as coordenadas são as seguintes

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

Aparentemente, o segundo sistema de coordenadas é definido por uma rotação geral da esfera

"Uma opção para esses parâmetros é:

  • A latitude geográfica em graus do polo sul do sistema de coordenadas, por exemplo;

  • A longitude geográfica em graus do polo sul do sistema de coordenadas, lambdap por exemplo;

  • O ângulo de rotação em graus sobre o novo eixo polar (medido no sentido horário ao olhar do polo sul para o polo norte) do sistema de coordenadas, assumindo que o novo eixo tenha sido obtido girando primeiro a esfera através de graus lambdap em torno do eixo polar geográfico e, em seguida, girando através de (90 + batida) graus para que o pólo sul se movesse ao longo do meridiano de Greenwich (anteriormente girado). "

mas ainda não sei como converter isso para o primeiro.

skd
fonte
2
Então, são esses dados da GRIB ? Se assim for, talvez precisemos de uma etiqueta grib.
Kirk Kuykendall
@skd os links do ECMWF não parecem ser válidos. Você pode editar?
gansub
@ gansub Eu editei os links. Não sei se as informações são exatamente as mesmas desde há muito tempo, mas acredito que o novo link possa fornecer algum contexto para referência futura.
Skd #
@skd quando você diz angle=0.0, você quer dizer o rolamento ? Eu tenho um arquivo netcdf com as coordenadas dos pólos girados, mas não há menção de nenhum ângulo.
FaCoffee
@ CF84 Na verdade, não tenho certeza. Eu acho que se não houver menção ao ângulo, então é o mesmo que angle = 0
skd

Respostas:

24

A reversão manual da rotação deve fazer o truque; deve haver uma fórmula para girar sistemas de coordenadas esféricas em algum lugar, mas como não consigo encontrá-la, aqui está a derivação ( ' marca o sistema de coordenadas giradas; as coordenadas geográficas normais usam símbolos simples):

Primeiro converta os dados no segundo conjunto de dados de esférico (lon ', lat') para (x ', y', z ') usando:

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

Em seguida, use duas matrizes de rotação para girar o segundo sistema de coordenadas para que coincida com o primeiro sistema 'normal'. Giraremos os eixos das coordenadas, para que possamos usar as matrizes de rotação dos eixos . Precisamos reverter o sinal na matriz to para corresponder ao sentido de rotação usado na definição ECMWF, que parece ser diferente da direção positiva padrão.

Como estamos desfazendo a rotação descrita na definição do sistema de coordenadas, primeiro giramos por ϑ = - (90 + lat0) = -55 graus ao redor do eixo y '(ao longo do meridiano de Greenwich girado) e depois por φ = - lon0 = +15 graus ao redor do eixo 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')

Expandido, isso se torna:

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

Em seguida, converta novamente para 'normal' (lat, lon) usando

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

Se você não possui atan2, pode implementá-lo usando atan (y / x) e examinando os sinais de x e y

Certifique-se de converter todos os ângulos em radianos antes de usar as funções trigonométricas, ou obterá resultados estranhos; converter de volta para graus no final, se é isso que você prefere ...

Exemplo (coordenadas de esfera girada ==> coordenadas geográficas padrão):

  • o pólo sul do CS girado é (lat0, lon0)

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

  • o principal meridiano do CS girado é o meridiano de -15 ° em termos geográficos (girado 55 ° em direção ao norte)

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

  • a simetria requer que ambos os equadores se cruzem a 90 ° / -90 ° no novo CS ou a 75 ° / -105 ° em coordenadas geográficas

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

EDIT: Reescreva a resposta graças a um comentário muito construtivo do whuber: as matrizes e a expansão estão agora sincronizadas, usando sinais adequados para os parâmetros de rotação; referência adicionada à definição das matrizes; removeu atan (y / x) da resposta; adicionados exemplos de conversão.

EDIT 2: É possível derivar expressões para o mesmo resultado sem transformação explícita no espaço cartesiano. A x, y, zem resultado pode ser substituído com as suas expressões correspondentes, e o mesmo pode ser repetido para x', y'e z'. Após aplicar algumas identidades trigonométricas, surgem as seguintes expressões de etapa única:

lat = arcsin(cos(ϑ) sin(lat') - cos(lon') sin(ϑ) cos(lat'))
lon = atan2(sin(lon'), tan(lat') sin(ϑ) + cos(lon') cos(ϑ)) - φ
mkadunc
fonte
1
A ideia é boa, mas alguns detalhes precisam ser corrigidos. lon0 = -15, não +15. Todas as três linhas na expansão do produto da matriz estão incorretas. ATan2 (ou equivalente) deve ser usado, modificado para retornar qualquer longitude razoável quando x = y = 0. Observe que, como x ^ 2 + y ^ 2 + z ^ 2 = 1, no final você obtém simplesmente lat = Arcsin (z).
whuber
1
Obrigado. Corrigi a resposta para pelo menos corrigir a matemática. As rotações agora devem coincidir com a descrição na definição de CS, mas é difícil ter certeza do sinal sem um exemplo (além da posição do polo sul).
Mkadunc 20/09/11
Bem feito! Estou surpreso que esta resposta não esteja recebendo mais votos, porque fornece material útil e difícil de encontrar.
whuber
É realmente muito difícil encontrar material, muito obrigado pela resposta. Acabei usando este software code.zmaw.de/projects/cdo para converter de uma grade rotacionada para uma grade regular. Meu palpite é que ele primeiro transforma as coordenadas como nesta resposta e depois as interpola para fornecer os resultados nos pontos de uma grade retangular. Embora um pouco tarde, eu a deixo para referência futura.
skd 19/03/13
1
@alfe Não sou especialista em esferas de Bloch, mas o princípio se parece muito com o que fiz, mas, em vez de converter em espaço cartesiano com 3 coordenadas reais, a dica sugere a conversão em um espaço com 2 coordenadas imaginárias (o que significa 4 componentes reais) e executando a rotação lá. Desencadeado pelo seu comentário, juntei todas as expressões e adicionei um resultado no qual a etapa cartesiana intermediária não é mais aparente.
Mkadunc
6

Caso alguém esteja interessado, compartilhei um script MATLAB na troca de arquivos, transformando lat / lon regular em lat / lon girado e vice-versa: transformação de grade rotativa

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
fonte
Caso o link fique inoperante, você pode inserir o código para futuros leitores. Obrigado.
Michael Stimson
1
Claro - código inserido.
Simdeck #
2

Essa transformação também pode ser calculada com o software proj (usando linha de comando ou programaticamente), empregando o que o proj chama de tradução oblíqua ( ob_tran) aplicada a uma transformação de latlon. Os parâmetros de projeção a serem definidos são:

  • o_lat_p = latitude do polo norte => 35 ° no exemplo
  • lon_0 = longitude do polo sul => -15 ° no exemplo
  • o_lon_p = 0

Além disso, -m 57.2957795130823é necessário (180 / pi) para considerar os valores projetados em graus.

Replicar os exemplos propostos por mkadunc dá o mesmo resultado (observe que aqui a ordem lon latnão é (lat,lon), os coordenados são digitados na entrada padrão, a saída é marcada por =>):

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

invprojO comando é usado para converter coordenadas "projetadas" (ie giradas) em geográficas, enquanto que projpara fazer o oposto.

Davide
fonte
1

Eu desenvolvi uma página asp.net para converter coordenadas de rotação para não rotação com base nos domínios CORDEX.

É baseado nos métodos acima. Você pode usá-lo livremente neste link:

Transformação manual de lat / lon girado em lat / lon regular

Sohrab kolsoomi ayask
fonte
O Cordex Data Extractor é um software de desktop do Windows para extrair dados do arquivo CORDEX NetCDF. O Cordex Data Extractor não precisa de um arquivo de ajuda, porque todos os processos foram realizados em códigos atrasados ​​e o usuário apenas digita datas, coordenadas e nome da variável. Por favor, assista a este vídeo: 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
fonte
0

Qual software você está usando? Todo software GIS terá a capacidade de mostrar as informações atuais do sistema / projeção. , o que pode ajudá-lo a obter o nome do seu sistema de coordenadas atual.

Além disso, se você estiver usando o ArcGIS, poderá usar a ferramenta Projeto para reprojetar o segundo conjunto de dados, importando configurações do primeiro.

ujjwalesri
fonte
2
Infelizmente não estou usando nenhum software. Esses são apenas conjuntos de dados da grade e vêm com as seguintes informações: - Para o primeiro: ecmwf.int/publications/manuals/d/gribapi/fm92/grib1/detail/… - Para o segundo: ecmwf.int/publications/ manuais / d / gribapi / fm92 / grib1 / detail / ...
SKD
Desde o ângulo de rotação é 0, acho que uma simples tradução deve alinhar o segundo conjunto de dados para o primeiro, como você disse acrescentando 15 a X e 35 a Y
ujjwalesri