Обсудить в форуме Комментариев 29Редактировать в вики
Обратная геодезическая задача — это нахождение начального направления и расстояния между двумя точками с известными координатами.
Содержание |
В качестве модели Земли принимается сфера с радиусом R, равным среднему радиусу земного эллипсоида. Аналогом прямой линии на плоскости является геодезическая линия на поверхности. На сфере геодезическая линия — дуга большого круга.
Введём следующие обозначения:
Линейное расстояние по дуге большого круга s связано со сферическим расстоянием σ формулой s = R σ.
Прямая и обратная геодезические задачи являются важными элементами более сложных геодезических задач.
На рисунке синим цветом выделены заданные элементы сферического треугольника, красным цветом неизвестные.
Существует великое множество подходов к решению поставленной задачи. Рассмотрим простой и надёжный векторный метод.
Последовательность решения:
Можно устранить второй пункт, если в первом заменить долготу λ₂ на разность долгот (λ₂ − λ₁).
Пример реализации алгоритма в виде функции языка Си:
/* * Решение обратной геодезической задачи * * Аргументы исходные: * pt1 - {широта, долгота} точки Q1 * pt2 - {широта, долгота} точки Q2 * * Аргументы определяемые: * azi - азимут начального направления * dist - расстояние (сферическое) */ void SphereInverse(double pt1[], double pt2[], double *azi, double *dist) { double x[3], pt[2]; SpherToCart(pt2, x); // сферические -> декартовы Rotate(x, pt1[1], 2); // первое вращение Rotate(x, M_PI_2 - pt1[0], 1); // второе вращение CartToSpher(x, pt); // декартовы -> сферические *azi = M_PI - pt[1]; *dist = M_PI_2 - pt[0]; return; }
Следует заметить, что прямая и обратная задача математически идентичны, и алгоритмы их решения зеркально отражают друг друга.
В данном случае в качестве сферических координат φ, λ подставим φ₂, λ₂.
Реализация на Си:
/* * Преобразование сферических координат в вектор * * Аргументы исходные: * y - {широта, долгота} * * Аргументы определяемые: * x - вектор {x, y, z} */ void SpherToCart(double y[], double x[]) { double p; p = cos(y[0]); x[2] = sin(y[0]); x[1] = p * sin(y[1]); x[0] = p * cos(y[1]); return; }
Представим оператор вращения вокруг оси X на угол θ в следующем виде:
Операторы вращения вокруг осей Y и Z получаются перестановкой символов.
Реализация вращения вокруг i-ой координатной оси на Си:
/* * Вращение вокруг координатной оси * * Аргументы: * x - входной/выходной 3-вектор * a - угол вращения * i - номер координатной оси (0..2) */ void Rotate(double x[], double a, int i) { double c, s, xj; int j, k; j = (i + 1) % 3; k = (i - 1) % 3; c = cos(a); s = sin(a); xj = x[j] * c + x[k] * s; x[k] = -x[j] * s + x[k] * c; x[j] = xj; return; }
В данном случае в роли сферических координат φ, λ окажутся углы (90° − σ), (180° − α₁).
Реализация на Си:
/* * Преобразование вектора в сферические координаты * * Аргументы исходные: * x - {x, y, z} * * Аргументы определяемые: * y - {широта, долгота} * * Возвращает: * длину вектора */ double CartToSpher(double x[], double y[]) { double p; p = hypot(x[0], x[1]); y[1] = atan2(x[1], x[0]); y[0] = atan2(x[2], p); return hypot(p, x[2]); }
Исходники вышеприведённых функций можно найти в архиве Sph.zip в файле sph.c. Кроме того, в файл sph.h включены следующие определения:
#define A_E 6371.0 // радиус Земли в километрах #define Degrees(x) (x * 57.29577951308232) // радианы -> градусы #define Radians(x) (x / 57.29577951308232) // градусы -> радианы
Теперь напишем программу, которая обращается к функции SphereInverse для решения обратной задачи:
#include <stdio.h> #include <stdlib.h> #include "sph.h" int main(int argc, char *argv[]) { char buf[1024]; double pt1[2], pt2[2]; double lat1, lon1, lat2, lon2, azi1, azi2, dist; while (fgets(buf, 1024, stdin) != NULL) { sscanf(buf, "%lf %lf %lf %lf", &lat1, &lon1, &lat2, &lon2); pt1[0] = Radians(lat1); pt1[1] = Radians(lon1); pt2[0] = Radians(lat2); pt2[1] = Radians(lon2); SphereInverse(pt2, pt1, &azi2, &dist); // Решение обратной задачи SphereInverse(pt1, pt2, &azi1, &dist); // Вычисление обратного азимута printf("%f\t%f\t%.4f\n", Degrees(azi1), Degrees(azi2), dist * A_E); } return 0; }
В архиве Sph.zip этот код находится в файле inv.c. Создадим исполняемый модуль inv компилятором gcc:
$ gcc -o inv inv.c sph.c -lm
Впрочем, в архиве есть Makefile. Для MS Windows готовую программу inv.exe можно найти в архиве Sph-win32.zip.
Программа читает данные из стандартного ввода консоли и отправляет результаты на стандартный вывод. Для чтения и записи файлов используются символы перенаправления потока «>» и «<» соответственно. Из каждой строки ввода программа считывает координаты двух точек φ₁, λ₁, φ₂, λ₂, которые должны быть в градусах, решает обратную задачу и записывает в строку вывода α₁, α₂, s (азимуты прямого и обратного направлений в градусах; расстояние между пунктами в километрах, а точнее, в единицах, определённых константой A_E).
Создадим файл inv.dat, содержащий одну строку данных:
30 0 52 54
После запуска программы
$ inv < inv.dat
получим α₁, α₂, s:
44.804060 262.415109 5001.1309
В архиве Sph-py.zip находятся скрипты на языке Питон. Выполнение скрипта в командной консоли:
$ python inv.py inv.dat
В пакет PROJ входит программа geod, предназначенная для решения прямых и обратных геодезических задач на сфере. Так выглядит команда обработки файла inv.dat:
$ geod +a=6371000 +b=6371000 -I -f "%f" -F "%.4f" +units=km inv.dat
Параметр +a определяет радиус сферы, -I — решение обратных задач, -f — формат вывода угловых величин, -F — формат вывода длин линий, +units — единица измерения расстояний. В результате получим идентичный вывод:
44.804060 -97.584891 5001.1309
Различие значений α₂ на 360° объясняется тем, что inv выводит азимуты в диапазоне от 0° до 360°, а geod от −180° до +180°.
Некоторые элементы альтернативных методов решения обратной задачи представлены в статье Вычисление расстояния и начального азимута между двумя точками на сфере.
В большинстве своём другие методы основаны на сферической тригонометрии. Многие из них используют вычисление σ или α₁ по таким функциям, как синус, косинус или гаверсинус. Это приводит к неоднозначности результатов вблизи особых значений, когда производная функции равна нулю. Такие методы не могут считаться универсальными.
К наиболее надёжным относится следующий способ:
В сферической тригонометрии углы и стороны должны быть в диапазоне [0, 180°]. Алгоритмизация формул требует анализа и обработки случаев, когда входные величины не попадают в эти рамки.
Обсудить в форуме Комментариев 29Редактировать в вики
Последнее обновление: 2020-05-10 07:46
Дата создания: 11.03.2014
Автор(ы): ErnieBoyd
© GIS-Lab и авторы, 2002-2021. При использовании материалов сайта, ссылка на GIS-Lab и авторов обязательна. Содержание материалов - ответственность авторов. (подробнее).