Comprendre les termes de la formule Longueur du diplôme?

13

Les calculateurs en ligne tels que http://www.csgnetwork.com/degreelenllavcalc.html (voir la source de la page) utilisent les formules ci-dessous pour obtenir des mètres par degré. Je comprends en général comment la distance par degré varie en fonction de la latitude, mais je ne comprends pas comment cela se traduit ci-dessous. Plus précisément, d'où viennent les constantes, les 3 termes "cos" dans chaque formule et les coefficients (2, 4, 6; 3 et 5) pour "lat"?

    // Set up "Constants"
    m1 = 111132.92;     // latitude calculation term 1
    m2 = -559.82;       // latitude calculation term 2
    m3 = 1.175;         // latitude calculation term 3
    m4 = -0.0023;       // latitude calculation term 4
    p1 = 111412.84;     // longitude calculation term 1
    p2 = -93.5;         // longitude calculation term 2
    p3 = 0.118;         // longitude calculation term 3

    // Calculate the length of a degree of latitude and longitude in meters
    latlen = m1 + (m2 * Math.cos(2 * lat)) + (m3 * Math.cos(4 * lat)) +
            (m4 * Math.cos(6 * lat));
    longlen = (p1 * Math.cos(lat)) + (p2 * Math.cos(3 * lat)) +
                (p3 * Math.cos(5 * lat));
Brent
la source
3
Sur un cercle, les termes de la forme cos (m * x) pour m = 0, 1, 2, ... jouent le même rôle que les monômes 1, x, x ^ 2, x ^ 3, ... font pour Taylor série sur la ligne. Lorsque vous voyez une expansion de ce type, vous pouvez penser de la même manière: chaque terme donne une approximation d'ordre supérieur à une fonction. Habituellement, ces séries trigonométriques sont infinies; mais en pratique, ils peuvent être tronqués dès que l'erreur d'approximation est acceptable. Une telle technologie se trouve sous le capot de chaque SIG, car de nombreuses projections sphéroïdales sont calculées à l'aide de ces séries.
whuber
Ceci est très utile pour calculer les distances où la distance entre les lignes de latitude varie, également utile pour aider à déterminer où tracer des points sur une carte Mercator si vous avez une grille x, y en superposition
Astuce: n'oubliez pas d'utiliser des radians pour lat(même si les variables résultantes latlenet longlensont en mètres par degré, pas en mètres par radian). Si vous utilisez des degrés pour lat, vous pouvez même vous retrouver avec une valeur négative pour longlen.
Luke Hutchison

Réponses:

23

Le rayon principal du sphéroïde WGS84 est a = 6378137 mètres et son aplatissement inverse est f = 298,257223563, d'où l'excentricité au carré est

e2 = (2 - 1/f)/f = 0.0066943799901413165.

Le rayon de courbure méridien à la latitude phi est

M = a(1 - e2) / (1 - e2 sin(phi)^2)^(3/2)

et le rayon de courbure le long du parallèle est

N = a / (1 - e2 sin(phi)^2)^(1/2)

De plus, le rayon du parallèle est

r = N cos(phi)

Ce sont des corrections multiplicatives des valeurs sphériques de M et N , toutes deux égales au rayon sphérique a , ce à quoi elles se réduisent lorsque e2 = 0.

Figure

Au point jaune à 45 degrés de latitude nord, le disque bleu de rayon M est le cercle osculateur ("baiser") en direction du méridien et le disque rouge de rayon N est le cercle osculateur en direction du parallèle: les deux les disques contiennent la direction "vers le bas" à ce stade. Ce chiffre exagère l'aplatissement de la terre de deux ordres de grandeur.

Les rayons de courbure déterminent les longueurs des degrés: lorsqu'un cercle a un rayon de R , son périmètre de longueur 2 pi R couvre 360 ​​degrés, d'où la longueur d'un degré est pi * R / 180. En substituant M et r pour R - - c'est-à-dire en multipliant M et r par pi / 180 - donne des formules exactes simples pour les longueurs en degrés.

Ces formules - qui sont basées uniquement sur les valeurs données de a et f (qui peuvent être trouvées à de nombreux endroits ) et sur la description du sphéroïde comme ellipsoïde de rotation - correspondent aux calculs de la question à 0,6 partie près. millions (quelques centimètres), ce qui correspond approximativement au même ordre de grandeur que les plus petits coefficients de la question, ce qui indique qu'ils sont d'accord. (L'approximation est toujours un peu faible.) Dans le graphique, l'erreur relative de longueur d'un degré de latitude est noire et celle de longitude est tiretée en rouge:

Figure

En conséquence, nous pouvons comprendre que les calculs de la question sont des approximations (via des séries trigonométriques tronquées) des formules données ci-dessus.


Les coefficients peuvent être calculés à partir de la série des cosinus de Fourier pour M et r en fonction de la latitude. Ils sont donnés en termes de fonctions elliptiques de e2, qui seraient trop compliquées à reproduire ici. Pour le sphéroïde WGS84, mes calculs donnent

  m1 = 111132.95255
  m2 = -559.84957
  m3 = 1.17514
  m4 = -0.00230
  p1 = 111412.87733
  p2 = -93.50412
  p3 = 0.11774
  p4 = -0.000165

(Vous pouvez deviner comment p4entre la formule. :) La proximité de ces valeurs avec les paramètres du code atteste de l'exactitude de cette interprétation. Cette approximation améliorée est précise à bien meilleure qu'une partie par milliard partout.


Pour tester cette réponse, j'ai exécuté du Rcode pour effectuer les deux calculs:

#
# Radii of meridians and parallels on a spheroid.  Defaults to WGS84 meters.
# Input is latitude (in degrees).
#
radii <- function(phi, a=6378137, e2=0.0066943799901413165) {
  u <- 1 - e2 * sin(phi)^2
  return(cbind(M=(1-e2)/u, r=cos(phi)) * (a / sqrt(u))) 
}
#
# Approximate calculation.  Same interface (but no options).
#
m.per.deg <- function(lat) {
  m1 = 111132.92;     # latitude calculation term 1
  m2 = -559.82;       # latitude calculation term 2
  m3 = 1.175;         # latitude calculation term 3
  m4 = -0.0023;       # latitude calculation term 4
  p1 = 111412.84;     # longitude calculation term 1
  p2 = -93.5;         # longitude calculation term 2
  p3 = 0.118;         # longitude calculation term 3

  latlen = m1 + m2 * cos(2 * lat) + m3 * cos(4 * lat) + m4 * cos(6 * lat);
  longlen = p1 * cos(lat) + p2 * cos(3 * lat) + p3 * cos(5 * lat);
  return(cbind(M.approx=latlen, r.approx=longlen))
}
#
# Compute the error of the approximation `m.per.deg` compared to the 
# correct formula and plot it as a function of latitude.
#
phi <- pi / 180 * seq(0, 90, 10)
names(phi) <- phi * 180 / pi
matplot(phi * 180 / pi, 10^6 * ((m.per.deg(phi) - radii(phi) * pi / 180) / 
       (radii(phi) * pi / 180)),
        xlab="Latitude (degrees)", ylab="Relative error * 10^6",lwd=2, type="l")

Le calcul exact avec radiipeut être utilisé pour imprimer des tableaux des longueurs de degrés, comme dans

zapsmall(radii(phi) * pi / 180)

La sortie est en mètres et ressemble à ceci (avec quelques lignes supprimées):

          M         r
0  110574.3 111319.49
10 110607.8 109639.36
20 110704.3 104647.09
...
80 111659.9  19393.49
90 111694.0      0.00

Les références

LM Bugayevskiy et JP Snyder, Projections cartographiques - Un manuel de référence. Taylor & Francis, 1995. (Annexe 2 et Annexe 4)

JP Snyder, Projections cartographiques - Un manuel de travail. USGS Professional Paper 1395, 1987. (Chapitre 3)

whuber
la source
Je ne sais pas pourquoi une approximation aussi compliquée d 'une simple paire de formules serait jamais utilisée ....
whuber
Quelle réponse approfondie et excellente! Cela semble correct; maintenant j'ai juste besoin de rafraîchir ce calcul pour le comprendre. :)
Brent
@Brent J'ai ajouté un chiffre pour vous aider à comprendre les mathématiques.
whuber
0

C'est la formule Haversine , bien qu'exprimée de manière étrange.

tmcw
la source
Ce n'est clairement pas la formule Haversine! C'est (lié à) une perturbation de celui-ci utilisé pour le sphéroïde. Il ne trouve même pas de distances entre des paires de points arbitraires, ce à quoi la formule de Haversine est utilisée (sur la sphère).
whuber
1
En d'autres termes, la formule de Haversine calcule la distance du grand cercle, et cette formule en est-elle une perturbation qui calcule la distance ellipsoïde plus précise?
Brent