Transformer manuellement le lat / lon tourné en lat / lon régulier?

24

Je dois d'abord préciser que je n'ai pas d'expérience préalable dans le domaine, donc je ne connais pas la terminologie technique. Ma question est la suivante:

J'ai deux jeux de données météorologiques:

  • Le premier a le système de coordonnées régulier (je ne sais pas s'il a un nom spécifique), allant de -90 à 90 et -180 à 180, et les pôles sont aux latitudes -90 et 90.

  • Dans le second, bien que cela devrait correspondre à la même région, j'ai remarqué quelque chose de différent: la latitude et la longitude n'étaient pas les mêmes, car elles ont un autre point de référence (dans la description, on appelle une grille pivotée ). Avec les paires lat / lon, viennent les informations suivantes: lat du pôle sud: -35,00, lon du pôle sud: -15,00, angle: 0,0.

J'ai besoin de transformer la deuxième paire de lon / lat en première. Cela pourrait être aussi simple que d'ajouter 35 aux latitudes et 15 aux longitudes, car l'angle est 0 et cela semble un simple décalage, mais je ne suis pas sûr.

Edit: Les informations que j'ai sur les coordonnées sont les suivantes

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

Apparemment, le deuxième système de coordonnées est défini par une rotation générale de la sphère

"Un choix pour ces paramètres est:

  • La latitude géographique en degrés du pôle sud du système de coordonnées, le tap par exemple;

  • La longitude géographique en degrés du pôle sud du système de coordonnées, lambdap par exemple;

  • Angle de rotation en degrés autour du nouvel axe polaire (mesuré dans le sens des aiguilles d'une montre en regardant du pôle sud au nord) du système de coordonnées, en supposant que le nouvel axe a été obtenu en faisant d'abord tourner la sphère en degrés lambdap autour de l'axe polaire géographique , puis tournant de (90 + thetap) degrés de sorte que le pôle sud se déplace le long du méridien de Greenwich (précédemment tourné). "

mais je ne sais toujours pas comment convertir cela en premier.

skd
la source
2
Alors, ces données GRIB ? Si c'est le cas, nous avons peut-être besoin d'une étiquette grib.
Kirk Kuykendall
@skd les liens ECMWF ne semblent pas valides. Pouvez-vous éditer?
gansub
@gansub J'ai édité les liens. Je ne sais pas si les informations sont exactement les mêmes depuis longtemps, mais je pense que le nouveau lien peut fournir un contexte pour référence future.
skd
@skd quand vous dites angle=0.0, voulez-vous dire le roulement ? J'ai un fichier netcdf avec les coordonnées des pôles tournés, mais il n'y a aucune mention d'angle.
FaCoffee
@ CF84 Je ne suis en fait pas sûr. Je suppose que s'il n'y a aucune mention de l'angle, alors c'est la même chose que angle = 0
skd

Réponses:

24

Inverser manuellement la rotation devrait faire l'affaire; il devrait y avoir une formule pour faire tourner les systèmes de coordonnées sphériques quelque part, mais comme je ne le trouve pas, voici la dérivation ( ' marque le système de coordonnées pivoté; les coordonnées géographiques normales utilisent des symboles simples):

Convertissez d'abord les données du deuxième ensemble de données de sphérique (lon ', lat') en (x ', y', z ') en utilisant:

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

Utilisez ensuite deux matrices de rotation pour faire pivoter le deuxième système de coordonnées de sorte qu'il coïncide avec le premier «normal». Nous allons faire pivoter les axes de coordonnées, afin que nous puissions utiliser les matrices de rotation des axes . Nous devons inverser le signe dans la matrice ϑ pour correspondre au sens de rotation utilisé dans la définition du CEPMMT, qui semble être différent de la direction positive standard.

Puisque nous annulons la rotation décrite dans la définition du système de coordonnées, nous effectuons d'abord une rotation de ϑ = - (90 + lat0) = -55 degrés autour de l'axe y '(le long du méridien de Greenwich pivoté), puis de φ = - lon0 = +15 degrés autour de l'axe 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')

Développé, cela devient:

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

Puis reconvertissez en «normal» (lat, lon) en utilisant

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

Si vous n'avez pas atan2, vous pouvez l'implémenter vous-même en utilisant atan (y / x) et en examinant les signes de x et y

Assurez-vous de convertir tous les angles en radians avant d'utiliser les fonctions trigonométriques, sinon vous obtiendrez des résultats étranges; reconvertissez en degrés à la fin si c'est ce que vous préférez ...

Exemple (coordonnées sphériques tournées ==> coordonnées géographiques standard):

  • le pôle sud du CS tourné est (lat0, lon0)

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

  • le méridien principal du CS tourné est le méridien géographique de -15 ° (tourné de 55 ° vers le nord)

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

  • la symétrie nécessite que les deux équateurs se croisent à 90 ° / -90 ° dans le nouveau CS, ou 75 ° / -105 ° en coordonnées géographiques

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

EDIT: Réécriture de la réponse grâce à un commentaire très constructif de whuber: les matrices et l'expansion sont maintenant synchronisées, en utilisant des signes appropriés pour les paramètres de rotation; ajout d'une référence à la définition des matrices; supprimé atan (y / x) de la réponse; ajouté des exemples de conversion.

EDIT 2: Il est possible de dériver des expressions pour le même résultat sans transformation explicite dans l'espace cartésien. Le x, y, zdans le résultat peut être substitué par leurs expressions correspondantes, et la même chose peut être répétée pour x', y'et z'. Après avoir appliqué des identités trigonométriques, les expressions à étape unique émergent:

lat = arcsin(cos(ϑ) sin(lat') - cos(lon') sin(ϑ) cos(lat'))
lon = atan2(sin(lon'), tan(lat') sin(ϑ) + cos(lon') cos(ϑ)) - φ
mkadunc
la source
1
L'idée est bonne, mais certains détails doivent être corrigés. lon0 = -15, pas +15. Les trois lignes de l'extension du produit matriciel sont incorrectes. ATan2 (ou son équivalent) doit être utilisé, modifié pour renvoyer toute longitude raisonnable lorsque x = y = 0. Notez que parce que x ^ 2 + y ^ 2 + z ^ 2 = 1, à la fin, vous obtenez simplement lat = Arcsin (z).
whuber
1
Merci. J'ai corrigé la réponse pour au moins rendre les calculs corrects. Les rotations devraient maintenant correspondre à la description de la définition CS, mais il est difficile d'être certain de leur signe sans exemple (autre que la position du pôle sud).
mkadunc
Bien joué! Je suis surpris que cette réponse n'obtienne pas plus de votes, car elle fournit du matériel utile et difficile à trouver.
whuber
Il est en effet très difficile de trouver du matériel, merci beaucoup pour la réponse. J'ai fini par utiliser ce logiciel code.zmaw.de/projects/cdo pour convertir d'une grille pivotée en une grille régulière. Je suppose qu'il transforme d'abord les coordonnées comme dans cette réponse, puis les interpole afin de donner les résultats aux points d'une grille rectangulaire. Bien qu'un peu en retard, je la laisse pour référence future.
skd
1
@alfe Je ne suis pas un expert des sphères de Bloch, mais le principe ressemble beaucoup à ce que j'ai fait, mais au lieu de convertir en espace cartésien avec 3 coordonnées réelles, l'indice suggère de convertir en un espace avec 2 coordonnées imaginaires (ce qui signifie 4 composants réels) et y exécuter la rotation. Déclenché par votre commentaire, j'ai rassemblé toutes les expressions et ajouté un résultat dans lequel l'étape cartésienne intermédiaire n'est plus apparente.
mkadunc
6

Dans le cas où quelqu'un est intéressé, j'ai partagé un script MATLAB sur l'échange de fichiers transformant le lat / lon régulier en lat / lon tourné et vice versa: Transformation de grille pivotée

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
la source
Au cas où le lien disparaîtrait, veuillez insérer le code pour les futurs lecteurs. Merci.
Michael Stimson
1
Bien sûr - code inséré.
simondk
2

Cette transformation peut également être calculée avec un logiciel proj (en utilisant la ligne de commande ou par programme) en utilisant ce que proj appelle une traduction oblique ( ob_tran) appliquée à une transformation latlon. Les paramètres de projection à définir sont:

  • o_lat_p = latitude du pôle nord => 35 ° dans l'exemple
  • lon_0 = longitude du pôle sud => -15 ° dans l'exemple
  • o_lon_p = 0

aditionellement, -m 57.2957795130823 (180 / pi) est nécessaire pour prendre en compte les valeurs projetées en degrés.

La réplication des exemples proposés par mkadunc donne le même résultat (notez que ici l'ordre n'est lon latpas (lat,lon), les coodinates sont tapés en entrée standard, la sortie est marquée par =>):

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

invprojLa commande est utilisée pour convertir des coordonnées "projetées" (c'est-à-dire tournées) en coordonnées géographiques, tout projen faisant le contraire.

Davide
la source
1

J'ai développé une page asp.net pour convertir les coordonnées de tourné en non tourné en fonction des domaines CORDEX.

Il est basé sur les méthodes ci-dessus. Vous pouvez l'utiliser librement dans ce lien:

Transformation manuelle de lat / lon tourné en lat / lon régulier

Sohrab kolsoomi ayask
la source
Cordex Data Extractor est un logiciel de bureau Windows pour extraire des données du fichier CORDEX NetCDF. Cordex Data Extractor n'a pas besoin de fichier d'aide car tous les processus ont été effectués en arrière-plan et l'utilisateur saisit simplement les dates, les coordonnées et le nom de la variable. Veuillez regarder cette vidéo: 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
la source
0

Quel logiciel utilisez-vous? Chaque logiciel SIG aura la possibilité de vous montrer les informations système / projection actuelles coordonnées. , ce qui peut vous aider à obtenir le nom de votre système de coordonnées actuel.

De plus, si vous utilisez ArcGIS, vous pouvez utiliser l' outil Projet pour reprojeter le deuxième jeu de données, en important les paramètres du premier.

ujjwalesri
la source
2
Malheureusement, je n'utilise aucun logiciel. Ce ne sont que des jeux de données de grille et ils contiennent les informations suivantes: - Pour le premier: ecmwf.int/publications/manuals/d/gribapi/fm92/grib1/detail/… - Pour le second: ecmwf.int/publications/ manuels / d / gribapi / fm92 / grib1 / detail /…
skd
Puisque l'angle de rotation est de 0, je pense qu'une simple traduction devrait aligner le deuxième ensemble de données sur le premier, comme vous l'avez dit en ajoutant 15 à X et 35 à Y
ujjwalesri