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

24

Primero debo aclarar que no tengo experiencia previa en el campo, por lo que no conozco la terminología técnica. Mi pregunta es la siguiente:

Tengo dos conjuntos de datos meteorológicos:

  • El primero tiene el sistema de coordenadas regular (no sé si tiene un nombre específico), que va desde -90 a 90 y -180 a 180, y los polos están en las latitudes -90 y 90.

  • En el segundo, aunque debería corresponder a la misma región, noté algo diferente: la latitud y la longitud no eran las mismas, ya que tienen otro punto de referencia (en la descripción se llama una cuadrícula girada ). Junto con los pares lat / lon, viene la siguiente información: lat del polo sur: -35.00, polo sur lon: -15.00, ángulo: 0.0.

Necesito transformar el segundo par de lon / lat en el primero. Podría ser tan simple como agregar 35 a las latitudes y 15 a las longitudes, ya que el ángulo es 0 y parece un cambio simple, pero no estoy seguro.

Editar: la información que tengo sobre las coordenadas es la siguiente

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

Aparentemente, el segundo sistema de coordenadas está definido por una rotación general de la esfera.

"Una opción para estos parámetros es:

  • La latitud geográfica en grados del polo sur del sistema de coordenadas, por ejemplo, el toque;

  • La longitud geográfica en grados del polo sur del sistema de coordenadas, lambdap por ejemplo;

  • El ángulo de rotación en grados sobre el nuevo eje polar (medido en el sentido de las agujas del reloj cuando se mira desde el polo sur al polo norte) del sistema de coordenadas, suponiendo que el nuevo eje se haya obtenido primero girando la esfera a través de grados lambdap sobre el eje polar geográfico , y luego girando (90 + thetap) grados para que el polo sur se moviera a lo largo del meridiano de Greenwich (previamente girado) ".

pero aún no sé cómo convertir esto al primero.

skd
fuente
2
Entonces, ¿son estos datos GRIB ? Si es así, tal vez necesitamos una etiqueta grib.
Kirk Kuykendall
@skd los enlaces ECMWF no parecen ser válidos. ¿Puedes editar?
gansub
@gansub He editado los enlaces. No sé si la información es exactamente la misma ya que ha pasado mucho tiempo, pero creo que el nuevo enlace puede proporcionar algún contexto para referencia futura.
skd
@skd cuando dices angle=0.0, ¿te refieres al rodamiento ? Tengo un archivo netcdf con las coordenadas del polo girado, pero no se menciona ningún ángulo.
FaCoffee
@ CF84 En realidad no estoy seguro. Supongo que si no se menciona el ángulo, entonces es lo mismo que angle = 0
skd

Respuestas:

24

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 alguna parte, pero como no puedo encontrarla, aquí está la derivación ( ' marca el sistema de coordenadas rotado; las coordenadas geográficas normales usan símbolos simples):

Primero convierta los datos en el segundo conjunto de datos de 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'. Vamos a girar los ejes de coordenadas, por lo que podemos usar las matrices de rotación de ejes . Necesitamos invertir el signo en la matriz ϑ para que coincida con el sentido de rotación utilizado en la definición ECMWF, que parece ser diferente de la dirección positiva estándar.

Como estamos deshaciendo la rotación descrita en la definición del sistema de coordenadas, primero rotamos por ϑ = - (90 + lat0) = -55 grados alrededor del eje y '(a lo largo del meridiano de Greenwich rotado) 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')

Ampliado, 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 nuevamente a 'normal' (lat, lon)

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

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

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

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

  • el polo sur del CS rotado es (lat0, lon0)

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

  • el primer meridiano del CS rotado es el meridiano -15 ° en posición geográfica (girado 55 ° hacia el norte)

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

  • la simetría requiere que ambos ecuadores se crucen a 90 ° / -90 ° en el nuevo CS, o 75 ° / -105 ° en 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 ahora están sincronizadas, utilizando signos adecuados para los parámetros de rotación; referencia agregada a la definición de las matrices; eliminó atan (y / x) de la respuesta; ejemplos agregados de conversión.

EDIT 2: es posible derivar expresiones para el mismo resultado sin una transformación explícita en el espacio cartesiano. El x, y, zen el resultado puede estar sustituido con sus correspondientes expresiones, y el mismo se puede repetir para x', y'y z'. Después de aplicar algunas identidades trigonométricas, surgen 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(ϑ)) - φ
mkadunc
fuente
1
La idea es buena, pero algunos detalles necesitan ser arreglados. lon0 = -15, no +15. Las tres líneas en la expansión del producto matriz son incorrectas. Se debe usar ATan2 (o su equivalente) , modificado para devolver cualquier longitud razonable cuando x = y = 0. Tenga en cuenta que debido a que x ^ 2 + y ^ 2 + z ^ 2 = 1, al final obtendrá simplemente lat = Arcsin (z).
whuber
1
Gracias. Arregle la respuesta para al menos corregir las matemáticas. Las rotaciones ahora deberían coincidir con la descripción en la definición CS, pero es difícil estar seguro acerca de su signo sin un ejemplo (aparte de la posición del polo sur).
mkadunc
¡Bien hecho! Me sorprende que esta respuesta no obtenga más votos, porque proporciona material útil y difícil de encontrar.
whuber
De hecho, es muy difícil encontrar material, muchas gracias por la respuesta. Terminé usando este software code.zmaw.de/projects/cdo para convertir de una grilla rotada a una grilla regular. Supongo que primero transforma las coordenadas como en esta respuesta y luego las interpola para dar los resultados en los puntos de una cuadrícula rectangular. Aunque un poco tarde, la dejo para referencia futura.
skd
1
@alfe No soy un experto en esferas de Bloch, pero el principio se parece mucho a lo que he hecho, pero en lugar de convertir al espacio cartesiano con 3 coordenadas reales, la sugerencia sugiere convertir a un espacio con 2 coordenadas imaginarias (lo que significa 4 componentes reales) y ejecutando la rotación allí. Activado por su comentario, puse todas las expresiones juntas y agregué un resultado en el que el paso cartesiano intermedio ya no es aparente.
mkadunc
6

En caso de que alguien esté interesado, he compartido un script MATLAB en el intercambio de archivos transformando lat / lon regular a lat / lon girado y viceversa: transformación de cuadrícula 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];
simondk
fuente
En caso de que el enlace se cierre, puede insertar el código para futuros lectores. Gracias.
Michael Stimson
1
Claro - código insertado.
simondk
2

Esta transformación también se puede calcular con proj de software (ya sea usando la línea de comandos o mediante programación) mediante el empleo de lo que llama proj traducción oblicua ( ob_tran) aplicada a una transformación LatLon. Los parámetros de proyección a configurar son:

  • o_lat_p = latitud del polo norte => 35 ° en el ejemplo
  • lon_0 = longitud del polo sur => -15 ° en el ejemplo
  • o_lon_p = 0

además, -m 57.2957795130823se requiere (180 / pi) para considerar los valores proyectados en grados.

Al replicar los ejemplos propuestos por mkadunc se obtiene el mismo resultado (observe que aquí el orden lon latno es (lat,lon), las coordenadas se escriben en la entrada estándar, la salida se marca con =>):

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

invprojEl comando se utiliza para convertir de coordenadas "proyectadas" (es decir, rotadas) a geográficas, mientras que projpara hacer lo contrario.

Davide
fuente
1

He desarrollado una página asp.net para convertir coordenadas de rotadas a no rotadas en función de 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

Sohrab kolsoomi ayask
fuente
Cordex Data Extractor es un software de escritorio de Windows para extraer datos del archivo CORDEX NetCDF. Cordex Data Extractor no necesita el archivo de ayuda debido a que todos los procesos se han realizado en códigos subyacentes y el usuario solo ingresa fechas, coordenadas y nombre de variable. Mire este video: youtu.be/RmpZblZPXjI agrimetsoft.com/cordexDataExtractor.aspx
Sohrab kolsoomi ayask
1

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

PITÓN:

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))
usuario126158
fuente
0

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

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

ujjwalesri
fuente
2
Desafortunadamente no estoy usando ningún software. Esos son solo conjuntos de datos de cuadrícula y vienen con la siguiente información: - Para el primero: ecmwf.int/publications/manuals/d/gribapi/fm92/grib1/detail/… - Para el segundo: ecmwf.int/publications/ manuales / d / gribapi / fm92 / grib1 / detail / ...
skd
Dado que el ángulo de rotación es 0, creo que una simple traducción debe alinear el segundo conjunto de datos a la primera, como usted ha dicho la adición de 15 a X y 35 a Y
ujjwalesri