Задачи на сфере: прямая геодезическая задача
Прямая геодезическая задача — это нахождение положения точки по координатам исходного пункта и значениям начального направления и расстояния.
Общие положения
В качестве модели Земли принимается сфера с радиусом R, равным среднему радиусу земного эллипсоида. Введём следующие обозначения:
- φ — географическая широта,
- λ — географическая долгота,
- α — азимут дуги большого круга,
- σ — сферическое расстояние (длина дуги большого круга, выраженная в долях радиуса шара).
Линейное расстояние по дуге большого круга s связано со сферическим расстоянием σ формулой s = R σ.
Прямая и обратная геодезические задачи являются важными элементами более сложных геодезических задач.
Постановка задачи
- Исходные данные
- координаты пункта Q₁, начальное направление и расстояние на сфере — φ₁, λ₁, α₁, σ.
- Определяемые величины
- координаты пункта Q₂: φ₂, λ₂.
На рисунке синим цветом выделены заданные элементы сферического треугольника, красным цветом неизвестные.
Алгоритм
Существует огромное множество подходов к решению поставленной задачи. Рассмотрим простой и надёжный векторный метод.
Последовательность решения:
- преобразовать углы (90° − σ) и (180° − α₁) в декартовы координаты,
- развернуть координатные оси вокруг оси Y на угол (φ₁ − 90°),
- развернуть координатные оси вокруг оси Z на угол −λ₁,
- преобразовать декартовы координаты в сферические.
Если третий пункт пропустить, на выходе вместо долготы λ₂ получится разность долгот (λ₂ − λ₁), что упростит алгоритм. Останется только прибавить долготу первого пункта. С другой строны, благодаря третьему пункту долгота λ₂ всегда будет в диапазоне [−180°, +180°].
Пример реализации алгоритма в виде функции языка Си:
/*
* Решение прямой геодезической задачи
*
* Аргументы исходные:
* pt1 - {широта, долгота} точки Q₁
* azi - азимут начального направления
* dist - расстояние (сферическое)
*
* Аргументы определяемые:
* pt2 - {широта, долгота} точки Q₂
*/
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 = sqrt(x[0] * x[0] + x[1] * x[1]);
y[1] = atan2(x[1], x[0]);
y[0] = atan2(x[2], p);
return sqrt(p * p + x[2] * x[2]);
}
Пример программной реализации
Исходники вышеприведённых функций можно найти в архиве Sph.tar.gz в файле 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.tar.gz этот код находится в файле dir.c. Создадим исполняемый модуль dir компилятором gcc:
$ gcc -o dir dir.c sph.c -lm
Впрочем, в архиве есть Makefile. Для MS Windows готовую программу dir.exe можно найти в архиве Sph-win32.zip.
Программа читает данные со стандартного ввода консоли и отправляет результаты на стандартный вывод. Для чтения и записи файлов используются символы перенаправления потока «>» и «<» соответственно. Из каждой строки ввода программа считывает координаты двух точек φ₁, λ₁, φ₂, λ₂, которые должны быть в градусах, решает обратную задачу и записывает в строку вывода α₁, α₂, s (азимуты прямого и обратного направлений в градусах и расстояние между пунктами в единицах, определённых константой A_E).
Создадим файл inv.dat, содержащий одну строку данных:
30 0 52 54
После запуска программы
$ inv < inv.dat
получим α₁, α₂, s:
44.804060 262.415109 5001.1309
Решение обратной задачи средствами PROJ.4
В пакет PROJ.4 входит программа geod, предназначенная для решения прямых и обратных геодезических задач на сфере. Так выглядит команда обработки файла inv.dat:
$ geod +a=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°.
Альтернативные методы
Некоторые элементы альтернативных методов решения обратной задачи представлены в статье Вычисление расстояния и начального азимута между двумя точками на сфере.