Comment calculez-vous le rayon de la Terre à une latitude géodésique donnée?

19

(Je vois qu'il y a une équation sur wikipedia qui fait exactement ce que je demande mais il n'y a pas de références. Je n'ai aucun moyen de confirmer la validité de cette équation!)

Je comprends déjà la différence entre Latitude géocentrique et Latitude géodésique.

En supposant que les rayons semi-majeurs aet semi-mineurs connus bsont donnés. Comment calculez-vous le rayon à une latitude géodésique donnée?

J'ai besoin d'une sorte de confirmation d'expert (dérivation, lien vers la dérivation, confirmation de l'expert, explication, etc.).

Trevor Boyd Smith
la source

Réponses:

25

Cette question suppose un modèle ellipsoïdal de la terre. Sa surface de référence est obtenue en faisant tourner une ellipse autour de son petit axe (tracée verticalement par convention). Une telle ellipse n'est qu'un cercle qui a été étiré horizontalement d'un facteur a et verticalement d'un facteur b . En utilisant le paramétrage standard du cercle unitaire,

t --> (cos(t), sin(t))

(qui définit le cosinus et le sinus), on obtient un paramétrage

t --> (a cos(t), b sin(t)).

(Les deux composantes de cette paramétrisation décrivent un parcours autour de la courbe: elles précisent, en coordonnées cartésiennes, notre position au "temps" t .)

La latitude géodésique , f , de tout point est l'angle que "up" fait avec le plan équatorial. Lorsque a diffère de b , la valeur de f diffère de celle de t (sauf le long de l'équateur et aux pôles).

Figure

Dans cette image, la courbe bleue est un quadrant d'une telle ellipse (considérablement exagérée par rapport à l'excentricité de la Terre). Le point rouge dans le coin inférieur gauche est son centre. La ligne en pointillés désigne le rayon à un point de la surface. Sa direction "vers le haut" y est représentée par un segment noir: elle est, par définition, perpendiculaire à l'ellipse en ce point. En raison de l'excentricité exagérée, il est facile de voir que "up" n'est pas parallèle au rayon.

Dans notre terminologie, t est lié à l'angle fait par le rayon à l'horizontale et f est l'angle fait par ce segment noir. (Notez que n'importe quel point de la surface peut être vu de cette perspective. Cela nous permet de limiter à la fois t et f entre 0 et 90 degrés; leurs cosinus et sinus seront positifs, donc nous n'avons pas à nous soucier des négatifs racines carrées dans les formules.)

L'astuce consiste à convertir la paramétrisation t en une unité en termes de f , car en termes de t le rayon R est facile à calculer (via le théorème de Pythagore). Son carré est la somme des carrés des composantes du point,

R(t)^2 = a^2 cos(t)^2 + b^2 sin(t)^2.

Pour effectuer cette conversion, nous devons relier la direction "ascendante" f au paramètre t . Cette direction est perpendiculaire à la tangente de l'ellipse. Par définition, une tangente à une courbe (exprimée en vecteur) est obtenue en différenciant sa paramétrisation:

Tangent(t) = d/dt (a cos(t), b sin(t)) = (-a sin(t), b cos(t)).

(La différenciation calcule le taux de changement. Le taux de changement de notre position lorsque nous parcourons la courbe est, bien sûr, notre vitesse , et cela pointe toujours le long de la courbe.)

Faites-le pivoter de 90 degrés dans le sens des aiguilles d'une montre pour obtenir la perpendiculaire, appelée vecteur "normal":

Normal(t) = (b cos(t), a sin(t)).

La pente de ce vecteur normal, égale à (a sin (t)) / (b cos (t)) ("montée sur course"), est aussi la tangente de l'angle qu'il fait à l'horizontale, d'où

tan(f) = (a sin(t)) / (b cos(t)).

De manière équivalente,

(b/a) tan(f) = sin(t) / cos(t) = tan(t).

(Si vous avez une bonne idée de la géométrie euclidienne, vous pouvez obtenir cette relation directement à partir de la définition d'une ellipse sans passer par un trig ou un calcul, simplement en reconnaissant que les expansions horizontales et verticales combinées de a et b respectivement ont pour effet de changer toutes les pentes par ce facteur b / a .)

Regardez à nouveau la formule pour R (t) ^ 2: nous connaissons a et b - ils déterminent la forme et la taille de l'ellipse - nous n'avons donc qu'à trouver cos (t) ^ 2 et sin (t) ^ 2 en termes de f , ce que l'équation précédente nous permet de faire facilement:

cos(t)^2 = 1/(1 + tan(t)^2) 
         = 1 / (1 + (b/a)^2 tan(f)^2) 
         = a^2 / (a^2 + b^2 tan(f)^2);
sin(t)^2 = 1 - cos(t)^2 
         = b^2 tan(f)^2 / (a^2 + b^2 tan(f)^2).

(Lorsque tan (f) est infini, nous sommes au pôle, il suffit donc de définir f = t dans ce cas.)

C'est la connexion dont nous avons besoin. Substituez ces valeurs pour cos (t) ^ 2 et sin (t) ^ 2 dans l'expression de R (t) ^ 2 et simplifiez pour obtenir

R(f)^2 = ( a^4 cos(f)^2 + b^4 sin(f)^2 ) / ( a^2 cos(f)^2 + b^2 sin(f)^2 ).

Une simple transformation montre que cette équation est la même que celle trouvée sur Wikipédia. Parce que a ^ 2 b ^ 2 = (ab) ^ 2 et (a ^ 2) ^ 2 = a ^ 4,

R(f)^2 = ( (a^2 cos(f))^2 + (b^2 sin(f))^2 ) / ( (a cos(f))^2 + (b sin(f))^2 )
whuber
la source
+1 .. sauf que je pense que la formule finale a un paren hors de propos ... ne devrait pas (b^4 sin(f))^2être changé en (b^4 sin(f)^2)?
Kirk Kuykendall
vraiment content qu'il y ait des experts autour de ce sujet =).
Trevor Boyd Smith,
Un fichier Geogebra (html) peut-il être publié sur ce site? J'ai un rayon de la verticale principale qui pourrait démontrer visuellement ce qui se passe.
Vous pouvez exporter l'original au format .png, @Dan: utilisez la boîte de dialogue Fichier | Exporter. Je recommande d'utiliser de grandes polices (16 ou 18 points semblent bien fonctionner) et de zoomer aussi loin que possible sur l'image.
whuber
Je suppose que l'interactivité sera alors perdue. La démo montre comment la variation des rayons et de la latitude d'intérêt change les propriétés.
3

Il est intéressant de constater que ma solution analphabète en mathématiques a fait le travail avec 5 minutes de réflexion et de codage, le facteur d'aplatissement ne devrait-il pas être pris en compte plutôt qu'un modèle elliptique parfait?

        double pRad = 6356.7523142;
        double EqRad = 6378.137;                      
        return pRad + (90 - Math.Abs(siteLatitude)) / 90 * (EqRad - pRad); 
Howard Grover
la source
1
Où pRad est le rayon polaire et EqRad est le rayon équatorial.
Stefan Steiger
c'est la seule réponse que j'ai pu lire. Il semble fonctionner pour moi.
Sean Bradley
1
Je vois que vous faites une interpolation linéaire du rayon entre le pôle et l'équateur. Bien qu'il n'y ait aucune raison de croire que l' interpolation linéaire est exacte , je vais l'utiliser comme "assez bonne" pour la Terre, compte tenu de son léger facteur d'aplatissement. BTW Je pense que c'est un peu plus facile à lire l'équivalent:, return E + (P - E) * Abs(Lat) / 90donc pas besoin d'avoir 90 - ...dans la formule.
ToolmakerSteve
2

entrez la description de l'image ici

C'est du moins la formule que j'ai trouvée au Centre d'analyse et d'évaluation des données des États-Unis (DAAC) pour le wiki du programme de modernisation du calcul haute performance (HPCMP) du ministère de la Défense (DoD) . Cela dit qu'ils ont beaucoup emprunté à l' entrée de Wikipédia . Pourtant, le fait qu'ils aient conservé cette formule devrait compter pour quelque chose.

RK
la source
Pouvez-vous fournir un lien vers le contenu?
Trevor Boyd Smith,
où φ est la latitude géodésique, et a (axe semi-majeur) et b (axe semi-mineur) sont, respectivement, le rayon équatorial et le rayon polaire. var a = 6378137; // m var b = 6356752.3142; // m en.wikipedia.org/wiki/Semi-major_and_semi-minor_axes en.wikipedia.org/wiki/World_Geodetic_System
Stefan Steiger