¿Transformar manualmente lat / lon girado a lat / lon regular?

Invertir manualmente la rotación debería hacer el truco; debería haber una fórmula para rotar los sistemas de coordenadas esféricas en algún lugar, pero como no puedo encontrarla, aquí está la derivación ( ' marca el sistema de coordenadas rotado; las coordenadas geográficas normales usan símbolos planos):

Primero convierta los datos en el segundo conjunto de datos desde esférico (lon ', lat') a (x ', y', z ') usando:

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

Luego use dos matrices de rotación para rotar el segundo sistema de coordenadas para que coincida con el primero "normal". Giraremos los ejes de coordenadas, por lo que podemos utilizar las matrices de rotación del eje . Necesitamos invertir el signo en la matriz para que coincida con el sentido de rotación utilizado en la definición de ECMWF, que parece ser diferente de la dirección positiva estándar.

Ya que estamos deshaciendo la rotación descrita en la definición del sistema de coordenadas, primero giramos por ϑ = - (90 + lat0) = -55 grados alrededor del eje y '(a lo largo del meridiano de Greenwich girado) y luego por φ = -lon0 = +15 grados alrededor del eje 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, esto se convierte en:

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

Luego convierta de nuevo a 'normal' (lat, lon) usando

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

Si no tiene atan2, puede implementarlo usted mismo utilizando atan (y / x) y examinando los signos de x y y

Asegúrate de convertir todos los ángulos a radianes antes de usar las funciones trigonométricas, o obtendrás resultados extraños; convertir de nuevo a grados al final si eso es lo que prefieres ...

Ejemplo (coordenadas de esfera giradas == > coordenadas geográficas estándar):

• el polo sur de la CS girada es (lat0, lon0)

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

• el meridiano principal de la CS girada es el meridiano -15 ° en geográfico (rotado 55 ° hacia el norte)

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

• La simetría

  requiere que ambos ecuadores se intersecten a 90 ° / -90 ° en la nueva CS, o 75 ° / -105 ° en las coordenadas geográficas

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

EDITAR: Reescribió la respuesta gracias a un comentario muy constructivo de whuber: las matrices y la expansión están ahora sincronizadas, utilizando signos adecuados para los parámetros de rotación; Se agregó referencia a la definición de las matrices; eliminado atan (y / x) de la respuesta; Se agregaron ejemplos de conversión.

EDIT 2: es posible derivar expresiones para el mismo resultado sin la transformación explícita en el espacio cartesiano. El x , y , z en el resultado se puede sustituir con sus expresiones correspondientes, y lo mismo se puede repetir para x' , y' y z' . Después de aplicar algunas identidades trigonométricas, emergen las siguientes expresiones de un solo paso:

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

En caso de que alguien esté interesado, he compartido un script de MATLAB en el intercambio de archivos que transforma lat / lon regular a lat / rotado y viceversa: Transformada de grid girada

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

He desarrollado una página asp.net para convertir coordenadas de rotadas a no rotadas según los dominios CORDEX.

Se basa en los métodos anteriores. Puedes usarlo libremente en este enlace:

Transformación manual de lat / lon girado a lat / lon regular

enlace

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

¿Qué software estás usando? Cada software GIS tendrá la facilidad de mostrarle la información actual del sistema / proyección coordinada. , que puede ayudarlo a obtener el nombre de su sistema de coordenadas actual.

Además, si está utilizando ArcGIS, puede usar Proyecto herramienta para volver a proyectar el segundo conjunto de datos, importando la configuración del primero.