Grundlegendes zu Begriffen in der Formel für die Länge des Grades?

Online-Taschenrechner wie http://www.csgnetwork.com/degreelenllavcalc.html ( Seitenquelle anzeigen) verwenden die folgenden Formeln, um Meter pro Grad zu erhalten. Ich verstehe im Allgemeinen, wie die Entfernung pro Grad je nach geografischer Breite variiert, aber ich verstehe nicht, wie sich dies auf das Folgende auswirkt. Woher kommen die Konstanten, die 3 "cos" -Terme in jeder Formel und die Koeffizienten (2, 4, 6; 3 und 5) für "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));

Der Hauptradius des WGS84-Sphäroids beträgt a = 6378137 Meter und seine inverse Abflachung beträgt f = 298.257223563, wobei die quadratische Exzentrizität gleich ist

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

Der meridionale Krümmungsradius am Breitengrad phi beträgt

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

und der Krümmungsradius entlang der Parallele ist

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

Weiterhin ist der Radius der Parallele

r = N cos(phi)

Dies sind multiplikative Korrekturen der sphärischen Werte von M und N , die beide gleich dem sphärischen Radius a sind , worauf sie sich reduzieren, wenn e2 = 0 ist.

Am gelben Punkt bei 45 Grad nördlicher Breite ist die blaue Scheibe mit dem Radius M der Oszillationskreis ("Kusskreis") in Richtung des Meridians und die rote Scheibe mit dem Radius N der Oszillationskreis in Richtung der Parallele: beides Discs enthalten an dieser Stelle die Abwärtsrichtung. Diese Zahl übertrifft die Abflachung der Erde um zwei Größenordnungen.

Die Krümmungsradien bestimmen die Längen der Grad: wenn ein Kreis mit einem Radius von HAS R , dessen Umfang mit einer Länge von 2 pi R Abdeckungen 360 Grad, von wo aus die Länge von einem Grad ist pi * R / 180. Substituieren M und R für R - - das heißt, multiplizieren M und r werden mit pi / 180 - ergibt einfache exakte Formeln für die Gradlängen.

Diese Formeln - die ausschließlich auf den vorgegebenen Werten von a und f (die an vielen Stellen zu finden sind ) und der Beschreibung des Sphäroids als Rotationsellipsoid beruhen - stimmen mit den Berechnungen in der Frage auf 0,6 Teile pro Sekunde überein Millionen (einige Zentimeter), was ungefähr der Größenordnung der kleinsten Koeffizienten in der Frage entspricht, was darauf hinweist, dass sie übereinstimmen. (Die Annäherung ist immer etwas niedrig.) In der Darstellung ist der relative Längenfehler eines Breitengrads schwarz und der Längengrad rot gestrichelt:

Dementsprechend können wir die Berechnungen in der Frage als Annäherungen (über abgeschnittene trigonometrische Reihen) an die oben angegebenen Formeln verstehen.

Die Koeffizienten können aus der Fourier-Cosinus-Reihe für M und r als Funktionen des Breitengrads berechnet werden . Sie sind in Bezug auf gegeben elliptischen Funktionen von e2 angegeben, die zu chaotisch wären, um hier reproduziert zu werden. Für den WGS84-Sphäroid geben meine Berechnungen

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

(Sie können sich vorstellen, wie p4 die Formel eingegeben wird. :) Die Nähe dieser Werte zu den Parametern im Code bestätigt die Richtigkeit dieser Interpretation. Diese verbesserte Annäherung ist überall viel genauer als ein Teil pro Milliarde.

Um diese Antwort zu testen, habe ich ausgeführt R Code , um beide Berechnungen durchzuführen:

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

Die genaue Rechnung mit radii kann verwendet werden, um Tabellen der Längen von Graden wie in zu drucken

zapsmall(radii(phi) * pi / 180)

Die Ausgabe erfolgt in Metern und sieht wie folgt aus (einige Zeilen wurden entfernt):

          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

Verweise

LM Bugayevskiy und JP Snyder, Kartenprojektionen - Ein Referenzhandbuch. Taylor & Francis, 1995. (Anhang 2 und Anhang 4)

JP Snyder, Kartenprojektionen - Ein Arbeitshandbuch. USGS Professional Paper 1395, 1987. (Kapitel 3)

Das ist die Haversine-Formel , obwohl sie auf seltsame Weise ausgedrückt wird.