Вычисление расстояний на сфере

Обсуждение материалов сайта: вопросы, замечания, предложения
Аватара пользователя
Максим Дубинин
MindingMyOwnBusiness
Сообщения: 8902
Зарегистрирован: 06 окт 2003, 20:20
Статьи: 231
Проекты: 12/6
Репутация: 642
Ваше звание: NextGIS
Откуда: Москва
Контактная информация:

Вычисление расстояний на сфере

Сообщение Максим Дубинин » 07 июн 2006, 00:58

Обсуждение статьи "Вычисление расстояний на сфере"

http://gis-lab.info/qa/great-circles.html
пристегивайтесь, турбулентность прямо по курсу

Аватара пользователя
Максим Дубинин
MindingMyOwnBusiness
Сообщения: 8902
Зарегистрирован: 06 окт 2003, 20:20
Статьи: 231
Проекты: 12/6
Репутация: 642
Ваше звание: NextGIS
Откуда: Москва
Контактная информация:

Сообщение Максим Дубинин » 29 дек 2006, 23:05

Добавлено вычисление начального азимута
пристегивайтесь, турбулентность прямо по курсу

DLL
Новоприбывший
Сообщения: 1
Зарегистрирован: 12 ноя 2007, 13:01
Репутация: 0

Сообщение DLL » 12 ноя 2007, 13:12

Не понял формулы начального азимута на языке Avenue.
(также не увидел ее в примере Exel)
Нельзя ли ее расписать, аналогично формуле расстояний?

Эта часть понятна:
'получение координат точек в радианах
lat1 = pnt1.getY*pi/180
lat2 = pnt2.getY*pi/180
long1 = pnt1.getX*pi/180
long2 = pnt2.getX*pi/180

'косинусы и синусы широт и разниц долгот
cl1 = lat1.cos
cl2 = lat2.cos
sl1 = lat1.sin
sl2 = lat2.sin
delta = long2 - long1
cdelta = delta.cos
sdelta = delta.sin
Далее почему-то идут x, y, z ???
Непонятно, на схеме этого нет.

Потом расчет:

Понятно, как x, y
x = (cl1*sl2) - (sl1*cl2*cdelta)
y = sdelta*cl2
Непонятно как z (в основном трудность с пониманием синтаксиса Avenue):
z = (-y/x).ATan.AsDegrees
if (x < 0) then z = z+180 end
z = -(z + 180 mod 360 - 180).AsRadians ?????
Далее ясно
anglerad2 = z - ((2*pi)*((z/(2*pi)).floor))
angledeg = (anglerad2*180)/pi

tarantyl
Новоприбывший
Сообщения: 1
Зарегистрирован: 01 июн 2009, 12:03
Репутация: 0

Re: Вычисление расстояний на сфере

Сообщение tarantyl » 01 июн 2009, 12:11

Воспользовался данной статьей и написал функционал для mysql Вычиcление расстояния между точками: гаверсинусы и mysql

Аватара пользователя
Максим Дубинин
MindingMyOwnBusiness
Сообщения: 8902
Зарегистрирован: 06 окт 2003, 20:20
Статьи: 231
Проекты: 12/6
Репутация: 642
Ваше звание: NextGIS
Откуда: Москва
Контактная информация:

Re: Вычисление расстояний на сфере

Сообщение Максим Дубинин » 01 июн 2009, 17:23

tarantyl писал(а):Воспользовался данной статьей и написал функционал для mysql
это хорошо, но раз использовали статью, то почему бы не сослаться на нее в своей?
пристегивайтесь, турбулентность прямо по курсу

neilsagitov
Новоприбывший
Сообщения: 2
Зарегистрирован: 02 июн 2009, 09:56
Репутация: 0

Re: Вычисление расстояний на сфере

Сообщение neilsagitov » 02 июн 2009, 10:03

Что-то не так у меня считает или я чего-то не понял.
Для координат 0, 0 и 0, 180 расстояние должно быть равно половине экватора?
а получается 0.

Попробовал Сидней - Нью-Йорк (примерные координаты -34, 152 и 41, -74),
тоже неверно получается, 4 084 км, вместо 16 000.

С уважением, Наиль.

neilsagitov
Новоприбывший
Сообщения: 2
Зарегистрирован: 02 июн 2009, 09:56
Репутация: 0

Re: Вычисление расстояний на сфере

Сообщение neilsagitov » 02 июн 2009, 10:05

neilsagitov писал(а):Что-то не так у меня считает или я чего-то не понял.
Забыл написать, это про 3-ю формулу.

С уважением, Наиль.

_Vovka_
Новоприбывший
Сообщения: 1
Зарегистрирован: 14 дек 2011, 15:36
Репутация: 0

Re: Вычисление расстояний на сфере

Сообщение _Vovka_ » 14 дек 2011, 15:43

Немного смущает rad = 6372795, т.к. гугл о радиусе Земли говорит 6378.1км
Да и вообще, r=40000/2/pi? потому что

"Средний диаметр планеты примерно равен 12 742 км. Это 40 000 км/π, так как метр в прошлом определялся, как 1/10 000 000 расстояния от экватора до северного полюса через Париж"
http://ru.wikipedia.org/wiki/%D0%97%D0% ... 0%BB%D1%8F

И работающая заготовка реализации на VisualBasic:

Dim distance As Long
Dim distance1 As Double
Dim u2 As Double
Dim u As Double
Dim Dg2r As Double
Dim sh2r As Double
Dim Dg1r As Double
Dim Sh1r As Double

Dim Dg1 As String
Dim Sh1 As String
Dim Dg2 As String
Dim sh2 As String

Sh1 = Inputbox$({Введите широту 1}, {Расстояние между точками Земли}, {56°30'N})
If Right$(Ucase(Sh1), 1) <> {N} Then Sh1 = {-} + Sh1
Dg1 = Inputbox$({Введите долготу 1}, {Расстояние между точками Земли}, {53°0'E})
If Right$(Ucase(Dg1), 1) <> {E} Then Dg1 = {-} + Dg1
Sh2 = Inputbox$({Введите широту 2}, {Расстояние между точками Земли}, {53°40'N})
If Right$(Ucase(Sh2), 1) <> {N} Then Sh2 = {-} + Sh2
Dg2 = Inputbox$({Введите долготу 2}, {Расстояние между точками Земли}, {53°38'E})
If Right$(Ucase(Dg2), 1) <> {E} Then Dg2 = {-} + Dg2

Sh1r=(Strleft(Sh1, {°})+Strleft(Strright(Sh1, {°}), {'})/60)*Pi/180
Dg1r=(Strleft(Dg1, {°})+Strleft(Strright(Dg1, {°}), {'})/60)*Pi/180
Sh2r=(Strleft(Sh2, {°})+Strleft(Strright(Sh2, {°}), {'})/60)*Pi/180
Dg2r=(Strleft(Dg2, {°})+Strleft(Strright(Dg2, {°}), {'})/60)*Pi/180
U=Sin(Sh1r)*Sin(Sh2r)+Cos(Sh1r)*Cos(Sh2r)*Cos(Dg1r-Dg2r)
U2=U*U
Distance1=Atn(Sqr(1/U2-1))*20000/Pi
If Distance1<0 Then
Distance1=0
End If
If Sh1r=0 Or Dg1r=0 Or Sh2r=0 Or Dg2r=0 Then
Distance=0
Else
Distance=Int(Distance1)
End If

Msgbox Distance, 64

Игорь Белов
Гуру
Сообщения: 1406
Зарегистрирован: 04 янв 2011, 22:00
Статьи: 12
Проекты: 1
Репутация: 863
Откуда: Казань

Re: Вычисление расстояний на сфере

Сообщение Игорь Белов » 19 май 2012, 17:29

Формула гаверсинусов используется, чтобы избежать проблем с небольшими расстояниями.
Довольно старое заблуждение. На самом деле теорема гаверсинусов для малых расстояний не имеет значительных преимуществ перед теоремой косинусов. С точки зрения математики гаверсинус есть перенормированный косинус:

Код: Выделить всё

hav x = (1 − cos x) / 2
Неопределённость арккосинуса при малых значениях аргумента связана не с тем, что косинус близок к единице, сколько с тем, что производная косинуса близка к нулю. Очевидно, что производная гаверсинуса также будет близка к нулю.
Предполагалось повышение точности благодаря меньшим потерям машинного округления из-за близости гаверсинуса к нулю, однако расчёты показывают, что выигрыш мизерный.
Последний раз редактировалось Игорь Белов 19 май 2012, 20:40, всего редактировалось 3 раза.

Игорь Белов
Гуру
Сообщения: 1406
Зарегистрирован: 04 янв 2011, 22:00
Статьи: 12
Проекты: 1
Репутация: 863
Откуда: Казань

Re: Вычисление расстояний на сфере

Сообщение Игорь Белов » 19 май 2012, 20:01

В догонку расчёты в Intel'овской двенадцатибайтной арифметике:

Код: Выделить всё

  x            x_c − x      x_h − x
5.96046e−08  3.52909e−23 −4.41055e−24
2.98023e−08  4.41055e−24  1.10183e−24
1.49012e−08  5.50915e−25  1.37325e−25
7.45058e−09 −7.45058e−09  1.69637e−26
3.72529e−09 −3.72529e−09  2.01948e−27
1.86265e−09 −1.86265e−09  2.01948e−28
В первой колонке аргумент в радианах; во второй — отличие арккосинуса косинуса от исходного значения; в третьей — аналогичная разница для гаверсинуса, который есть синус квадрат половинного аргумента.
Видно, что существенная разница всё-таки появляется, когда косинус превращается в единицу при значении аргумента 7.45e−09. При этом погрешность арккосинуса равна по величине аргументу, и для обоих это соответствует дальности 0.047 метра на поверхности Земного шара.
При хранении промежуточных результатов в стандартном восьмибайтном представлении ситуация принципиально не меняется.
Признаю, точность обратного гаверсинуса несколько выше по сравнению с арккосинусом. Однако можно выбирать между максимальной погрешностью в пять сантиметров при расстоянии 5 сантиметров в вычислении арккосинуса и некоторым усложнением алгоритма и потерей быстродействия (обратный гаверсинус есть удвоенный арксинус корня) при практически одинаковой точности на бОльших расстояниях.
Между тем, в библиотеке GIS-Lab есть хорошая книга Степанова Сферическая тригонометрия, в которой предлагается два иных решения задачи. Возможно, стоит упомянуть о них в статье. Первое из них лишено неопределённости при любых значениях аргумента.

Аватара пользователя
Максим Дубинин
MindingMyOwnBusiness
Сообщения: 8902
Зарегистрирован: 06 окт 2003, 20:20
Статьи: 231
Проекты: 12/6
Репутация: 642
Ваше звание: NextGIS
Откуда: Москва
Контактная информация:

Re: Вычисление расстояний на сфере

Сообщение Максим Дубинин » 19 май 2012, 20:08

ErnieBoyd, как насчет я переношу свою статью в вики, а вы вносите добавки, какие сочтете нужными?
пристегивайтесь, турбулентность прямо по курсу

Игорь Белов
Гуру
Сообщения: 1406
Зарегистрирован: 04 янв 2011, 22:00
Статьи: 12
Проекты: 1
Репутация: 863
Откуда: Казань

Re: Вычисление расстояний на сфере

Сообщение Игорь Белов » 19 май 2012, 20:14

Замечательная статья. Непременно переносите.

Аватара пользователя
Чулан
Участник
Сообщения: 54
Зарегистрирован: 13 дек 2010, 11:07
Репутация: 0
Откуда: Москва
Контактная информация:

Re: Вычисление расстояний на сфере

Сообщение Чулан » 19 май 2013, 20:16

Изображение
Почему, если найти среднее арифметическое координат (в долях градусов) начала и конца линии (на карте Яндекс проведённой), оно покажет середину прямой линии, а не середину ортодромии (дуги большого круга)? Как найти координаты середины кратчайшего (дугообразно выглядящего на плоскости) расстояния?

Аватара пользователя
Максим Дубинин
MindingMyOwnBusiness
Сообщения: 8902
Зарегистрирован: 06 окт 2003, 20:20
Статьи: 231
Проекты: 12/6
Репутация: 642
Ваше звание: NextGIS
Откуда: Москва
Контактная информация:

Re: Вычисление расстояний на сфере

Сообщение Максим Дубинин » 19 май 2013, 20:35

ErnieBoyd, в лучших традициях, год прошел, статью я перенес, теперь редактируется.
пристегивайтесь, турбулентность прямо по курсу

Аватара пользователя
Максим Дубинин
MindingMyOwnBusiness
Сообщения: 8902
Зарегистрирован: 06 окт 2003, 20:20
Статьи: 231
Проекты: 12/6
Репутация: 642
Ваше звание: NextGIS
Откуда: Москва
Контактная информация:

Re: Вычисление расстояний на сфере

Сообщение Максим Дубинин » 19 май 2013, 20:38

Чулан, показывает середину правильно, так как среднее арифметическое означает, что вы рассматриваете координаты как линейные единицы, а не угловые. Хотите, чтобы среднее арифметическое показывало середину дуги? Используйте гномоническую проекцию и ведите расчеты в ней.
пристегивайтесь, турбулентность прямо по курсу

Ответить

Вернуться в «Материалы сайта»