Обсудить в форуме Комментариев 29Редактировать в вики
Прямая геодезическая задача — это нахождение положения точки по координатам исходного пункта и значениям начального направления и расстояния.
Содержание |
В качестве модели Земли принимается сфера с радиусом R, равным среднему радиусу земного эллипсоида. Аналогом прямой линии на плоскости является геодезическая линия на поверхности. На сфере геодезическая линия — дуга большого круга.
Введём следующие обозначения:
Линейное расстояние по дуге большого круга s связано со сферическим расстоянием σ формулой s = R σ.
Прямая и обратная геодезические задачи являются важными элементами более сложных геодезических задач.
На рисунке синим цветом выделены заданные элементы сферического треугольника, красным цветом неизвестные.
Существует великое множество подходов к решению поставленной задачи. Рассмотрим простой и надёжный векторный метод.
Последовательность решения:
Если третий пункт пропустить, на выходе вместо долготы λ₂ получится разность долгот (λ₂ − λ₁), что упростит алгоритм. Останется только прибавить долготу первого пункта. С другой строны, благодаря третьему пункту долгота λ₂ всегда будет в диапазоне [−180°, +180°].
Пример реализации алгоритма в виде функции языка Си:
/* * Решение прямой геодезической задачи * * Аргументы исходные: * pt1 - {широта, долгота} точки Q1 * azi - азимут начального направления * dist - расстояние (сферическое) * * Аргументы определяемые: * pt2 - {широта, долгота} точки Q2 */ void SphereDirect(double pt1[], double azi, double dist, double pt2[]) { double pt[2], x[3]; pt[0] = M_PI_2 - dist; pt[1] = M_PI - azi; SpherToCart(pt, x); // сферические -> декартовы Rotate(x, pt1[0] - M_PI_2, 1); // первое вращение Rotate(x, -pt1[1], 2); // второе вращение CartToSpher(x, pt2); // декартовы -> сферические return; }
Следует заметить, что прямая и обратная задача математически идентичны, и алгоритмы их решения зеркально отражают друг друга.
В данном случае в качестве сферических координат φ, λ подставим углы (90° − σ), (180° − α₁).
Реализация на Си:
/* * Преобразование сферических координат в вектор * * Аргументы исходные: * 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; }
В данном случае в роли сферических координат φ, λ окажутся φ₂, λ₂.
Реализация на Си:
/* * Преобразование вектора в сферические координаты * * Аргументы исходные: * 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) // градусы -> радианы
Теперь напишем программу, которая обращается к функции SphereDirect для решения прямой задачи:
#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, azi1, dist, azi2; while (fgets(buf, 1024, stdin) != NULL) { sscanf(buf, "%lf %lf %lf %lf", &lat1, &lon1, &azi1, &dist); pt1[0] = Radians(lat1); pt1[1] = Radians(lon1); SphereDirect(pt1, Radians(azi1), dist / A_E, pt2); // Решение прямой задачи SphereInverse(pt2, pt1, &azi2, &dist); // Вычисление обратного азимута printf("%f\t%f\t%f\n", Degrees(pt2[0]), Degrees(pt2[1]), Degrees(azi2)); } return 0; }
В архиве Sph.zip этот код находится в файле dir.c. Создадим исполняемый модуль dir компилятором gcc:
$ gcc -o dir dir.c sph.c -lm
Впрочем, в архиве есть Makefile. Для MS Windows готовую программу dir.exe можно найти в архиве Sph-win32.zip.
Программа читает данные из стандартного ввода консоли и отправляет результаты на стандартный вывод. Для чтения и записи файлов используются символы перенаправления потока «>» и «<» соответственно. Из каждой строки ввода программа считывает координаты первого пункта φ₁, λ₁, начальный азимут α₁ в градусах и расстояние s в километрах; решает прямую задачу; записывает в строку вывода координаты второго пункта φ₂, λ₂ и обратный азимут α₂ в градусах.
Создадим файл dir.dat, содержащий одну строку данных:
30 0 44.804060 5001.1309
После запуска программы
$ dir < dir.dat
получим φ₂, λ₂, α₂:
52.000000 54.000001 262.415109
В архиве Sph-py.zip находятся скрипты на языке Питон. Выполнение скрипта в командной консоли:
$ python dir.py dir.dat
В пакет PROJ входит программа geod, предназначенная для решения прямых и обратных геодезических задач на сфере. Так выглядит команда обработки файла dir.dat:
$ geod +a=6371000 +b=6371000 -f "%f" +units=km dir.dat
Параметр +a определяет радиус сферы, -f — формат вывода угловых величин, +units — единица измерения расстояний. В итоге получим идентичный результат:
52.000000 54.000001 -97.584891
Различие значений α₂ на 360° объясняется тем, что dir выводит азимуты в диапазоне от 0° до 360°, а geod от −180° до +180°.
С помощью geod можно также расставить промежуточные точки вдоль геодезической линии либо по дуге малого круга на заданном расстоянии от исходного пункта. В обоих случаях нужно задать положение начальной точки параметрами +lat_1, +lon_1 и либо координаты второй точки +lat_2, +lon_1, либо расстояние и азимут ко второй точке +S, +A. За подробностями обращайтесь к документации.
Большая часть других методов основана на сферической тригонометрии. Многие из них используют вычисление φ₂ или (λ₂ − λ₁) по таким функциям, как синус, косинус или гаверсинус. Это приводит к неоднозначности результатов вблизи особых значений, когда производная функции равна нулю. Такие методы не могут считаться универсальными.
К наиболее надёжным относится следующий способ:
В сферической тригонометрии углы и стороны должны быть в диапазоне [0, 180°]. Алгоритмизация формул требует анализа и обработки случаев, когда входные величины не попадают в эти рамки.
Обсудить в форуме Комментариев 29Редактировать в вики
Последнее обновление: 2020-05-10 07:46
Дата создания: 11.03.2014
Автор(ы): ErnieBoyd
© GIS-Lab и авторы, 2002-2021. При использовании материалов сайта, ссылка на GIS-Lab и авторов обязательна. Содержание материалов - ответственность авторов. (подробнее).