Понимание терминов в формуле длины степени?

13

Онлайн-калькуляторы, такие как http://www.csgnetwork.com/degreelenllavcalc.html (просмотр страницы источника), используют приведенные ниже формулы для получения метров за градус. Я в целом понимаю, как расстояние на градус варьируется в зависимости от местоположения широты, но я не понимаю, как это приводит к приведенному ниже. Более конкретно, откуда берутся константы, 3 члена «cos» в каждой формуле и коэффициенты (2, 4, 6; 3 и 5) для «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));
казарка
источник
3
На окружности члены вида cos (m * x) для m = 0, 1, 2, ... играют ту же роль, что и мономы 1, x, x ^ 2, x ^ 3, ... do для Тейлора серия на линии. Когда вы видите расширение такого рода, вы можете думать об этом одинаково: каждый член дает приближение высшего порядка для функции. Обычно такие тригонометрические ряды бесконечны; но при практическом использовании они могут быть усечены, как только погрешность аппроксимации приемлема. Некоторые такие технологии лежат в основе каждой ГИС, потому что многие сфероидальные проекции вычисляются с использованием таких рядов.
whuber
Это очень полезно для расчета расстояний, где расстояние между линиями широты меняется, также полезно, чтобы помочь определить, где наносить точки на карту меркатора, если у вас есть сетка x, y в качестве наложения
Совет: не забывайте использовать радианы для lat(хотя получающиеся переменные latlenи longlenуказаны в метрах на градус, а не в метрах на радиан). Если вы используете градусы для lat, вы можете даже получить отрицательное значение для longlen.
Люк Хатчисон

Ответы:

23

Главный радиус сфероида WGS84 составляет a = 6378137 метров, а его обратное уплощение равно f = 298.257223563, откуда квадрат эксцентриситета равен

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

Меридиональный радиус кривизны на широте фи является

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

а радиус кривизны вдоль параллели равен

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

Кроме того, радиус параллели равен

r = N cos(phi)

Это мультипликативные поправки к сферическим значениям M и N , которые равны сферическому радиусу a , к которому они уменьшаются, когда e2 = 0.

фигура

В желтой точке в 45 градусах северной широты синий диск радиуса M - это осциллирующий круг («поцелуй») в направлении меридиана, а красный диск радиуса N - это осциллирующий круг в направлении параллели: оба диски содержат направление «вниз» в этой точке. Эта цифра преувеличивает уплощение Земли на два порядка.

Радиусы кривизны определяют длины градусов: когда круг имеет радиус R , его периметр длины 2 pi R охватывает 360 градусов, откуда длина одного градуса равна pi * R / 180. Подставляя M и r для R - то есть умножение M и r на pi / 180 - дает простые точные формулы для длин градусов.

Эти формулы, которые основаны исключительно на заданных значениях a и f (которые можно найти во многих местах ) и описании сфероида как эллипсоида вращения, согласуются с расчетами в вопросе с точностью до 0,6 частей на миллион (несколько сантиметров), что примерно соответствует порядку величины самых маленьких коэффициентов в вопросе, что указывает на их согласие. (Приближение всегда немного низкое.) На графике относительная погрешность длины градуса широты черная, а долготы - пунктирная красная:

фигура

Соответственно, мы можем понимать, что вычисления в вопросе являются приближениями (через усеченные тригонометрические ряды) к формулам, приведенным выше.


Коэффициенты могут быть вычислены из ряда косинусов Фурье для M и r как функции широты. Они даны в терминах эллиптических функций e2, которые были бы слишком беспорядочными для воспроизведения здесь. Для сфероида WGS84 мои расчеты дают

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

(Вы можете догадаться, как p4входит в формулу. :) Близость этих значений к параметрам в коде свидетельствует о правильности этой интерпретации. Эта улучшенная аппроксимация точнее, чем одна часть на миллиард, везде.


Чтобы проверить этот ответ, я выполнил Rкод для выполнения обоих расчетов:

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

Точный расчет с radiiможет быть использован для печати таблиц длин градусов, как в

zapsmall(radii(phi) * pi / 180)

Вывод в метрах и выглядит следующим образом (с некоторыми удаленными линиями):

          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

Ссылки

Л.М. Бугаевский и Ю.П. Снайдер. Картографические проекции - справочное руководство. Taylor & Francis, 1995. (Приложение 2 и Приложение 4)

JP Снайдер, Карта проекций - рабочее руководство. USGS Professional Paper 1395, 1987. (Глава 3)

Whuber
источник
Я не знаю, почему когда-либо использовалось бы такое сложное приближение к простой паре формул ...
whuber
Какой полный, отличный ответ! Это кажется правильным; теперь мне просто нужно освежить в этом математику, чтобы понять это. :)
Брент
@ Брент Я добавил цифру, чтобы помочь тебе понять математику.
whuber
0

Это формула Haversine , хотя и выражена странным образом.

tmcw
источник
Это явно не формула Haversine! Это (связано с) его возмущение, используемое для сфероида. Он даже не находит расстояния между произвольными парами точек, для чего используется формула Хаверсайна (на сфере).
whuber
1
Другими словами, формула Хаверсайна вычисляет расстояние по большому кругу, и эта формула является его возмущением, которое вычисляет более точное расстояние эллипсоида?
Брент