Геопространственные координаты и расстояние в километрах
Это продолжение этого вопроса.
Кажется, я застрял на этом. По сути, мне нужно иметь возможность конвертировать туда и обратно для обращения к координатам либо в стандартной системе градусов, либо путем измерения расстояния к северу от южного полюса вдоль международной линии даты, а затем расстояния на восток, начиная с этой точки даты линия. Чтобы сделать это (а также некоторые более общие вещи для измерения расстояния), у меня есть один метод для определения расстояния между двумя точками широты и долготы, и другой метод, который берет точку широты и долготы, курс и расстояние и возвращает широта / долгота в конце этого курса.
Вот два статических метода, которые я определил:
/* Takes two lon/lat pairs and returns the distance between them in kilometers.
*/
public static double distance (double lat1, double lon1, double lat2, double lon2) {
double theta = toRadians(lon1-lon2);
lat1 = toRadians(lat1);
lon1 = toRadians(lon1);
lat2 = toRadians(lat2);
lon2 = toRadians(lon2);
double dist = sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2)*cos(theta);
dist = toDegrees(acos(dist)) * 60 * 1.1515 * 1.609344 * 1000;
return dist;
}
/* endOfCourse takes a lat/lon pair, a heading (in degrees clockwise from north), and a distance (in kilometers), and returns
* the lat/lon pair that would be reached by traveling that distance in that direction from the given point.
*/
public static double[] endOfCourse (double lat1, double lon1, double tc, double dist) {
double pi = Math.PI;
lat1 = toRadians(lat1);
lon1 = toRadians(lon1);
tc = toRadians(tc);
double dist_radians = toRadians(dist / (60 * 1.1515 * 1.609344 * 1000));
double lat = asin(sin(lat1) * cos(dist_radians) + cos(lat1) * sin(dist_radians) * cos(tc));
double dlon = atan2(sin(tc) * sin(dist_radians) * cos(lat1), cos(dist_radians) - sin(lat1) * sin(lat));
double lon = ((lon1-dlon + pi) % (2*pi)) - pi;
double[] endPoint = new double[2];
endPoint[0] = lat; endPoint[1] = lon;
return endPoint;
}
И вот функция, которую я использую, чтобы проверить это:
public static void main(String args[]) throws java.io.IOException, java.io.FileNotFoundException {
double distNorth = distance(0.0, 0.0, 72.0, 0.0);
double distEast = distance(72.0, 0.0, 72.0, 31.5);
double lat1 = endOfCourse(0.0, 0.0, 0.0, distNorth)[0];
double lon1 = endOfCourse(lat1, 0.0, 90.0, distEast)[1];
System.out.println("end at: " + lat1 + " / " + lon1);
return;
}
Значения "end at" должны быть ок. 72,0 / 31,5. Но вместо этого я получаю примерно 1,25 / 0,021.
Я предполагаю, что я, должно быть, упускаю что-то глупое, забывая конвертировать единицы или что-то в этом роде... Любая помощь будет с благодарностью
ОБНОВЛЕНИЕ 1:
Я (правильно) написал функцию расстояния, чтобы вернуть метры, но по ошибке написал километры в комментариях... что, конечно, смутило меня, когда я вернулся к нему сегодня. В любом случае, теперь это исправлено, и я исправил ошибку факторинга в методе endOfCourse, и я также понял, что забыл преобразовать обратно в градусы из радиан в этом методе. В любом случае: хотя кажется, что теперь я получаю правильное значение широты (71,99...), число долготы далеко (я получаю 3,54 вместо 11,5).
ОБНОВЛЕНИЕ 2: у меня была опечатка в тесте, как упомянуто ниже. Это сейчас исправлено в коде. Однако значение долготы по-прежнему неверно: теперь я получаю -11,34 вместо 11,5. Я думаю, что должно быть что-то не так с этими строками:
double dlon = atan2(sin(tc) * sin(dist_radians) * cos(lat1), cos(dist_radians) - sin(lat1) * sin(lat));
double lon = ((lon1-dlon + pi) % (2*pi)) - pi;
5 ответов
У вас серьезный случай магических чисел в коде. Выражение:
(60 * 1.1515 * 1.609344 * 1000)
появляется дважды, но нет особых объяснений этому. С некоторой помощью: 1.609344 - это количество километров в миле; 60 - количество минут в градусе; 1000 - это количество метров в километре; и 1.1515 - это количество уставных миль в морской миле (спасибо, DanM). Одна морская миля - это одна минута широты на экваторе.
Я предполагаю, что вы используете модель сферической земли, а не сфероидальную землю? Алгебра не достаточно сложна, чтобы быть сфероидальной.
Первая формула - преобразование между двумя парами широты и долготы - нечетна. Вам нужны и дельта-лат (Δλ), и дельта-лон (Δφ), чтобы разобраться в ответе. Далее расстояние между парами:
(60° N, 30° W), (60° N, 60° W)
(60° N, 60° W), (60° N, 90° W)
должно быть одинаковым - но я уверен, что ваш код дает разные ответы.
Итак, я думаю, вам нужно вернуться к справочным материалам по сферической тригонометрии и посмотреть, что вы делаете неправильно. (Мне понадобится некоторое время, чтобы найти мою книгу на эту тему - ее нужно будет распаковать из любой коробки, в которой она находится.)
[ ... проходит время... распаковка завершена... ]
Учитывая сферический треугольник с углами A, B, C в вершинах и сторонах a, b, c, противоположных этим вершинам (то есть сторона a находится от B до C и т. Д.), Формула косинуса имеет вид:
cos a = cos b . cos c + sin b . sin c . cos A
Применяя это к задаче, мы можем назвать две точки, заданные B и C, и создадим прямоугольный сферический треугольник с прямым углом к A.
Искусство ASCII в худшем:
+ C
/|
/ |
a / | b
/ |
/ |
/ |
B +------+ A
c
Сторона с равна разнице в долготе; сторона b равна разнице в широте; угол A равен 90°, поэтому cos A = 0. Поэтому я считаю, что уравнение для a:
cos a = cos Δλ . cos Δφ + sin Δλ . sin Δφ . cos 90°
a = arccos (cos Δλ . cos Δφ)
Угол a в радианах затем преобразуется в расстояние путем умножения на радиус Земли. В качестве альтернативы, если задано a в градусах (и долях градуса), то есть 60 морских миль на один градус, следовательно, 60 * 1.1515 статутных миль и 60 * 1.1515 * 1.609344 километра на один градус. Если вы не хотите расстояние в метрах, я не вижу необходимости в коэффициент 1000.
Пол Томблин указывает на Aviation Formulary v1.44 в качестве источника уравнения - и, действительно, он существует вместе с более численно устойчивой версией для случаев, когда разница в положении мала.
Переходя к базовой тригонометрии, мы также знаем, что:
cos (A - B) = cos A . cos B + sin A . sin B
Применяя это дважды в уравнении, которое я дал, вполне может оказаться в формуле в Авиационном формуляре.
(Моя ссылка: "Астрономия: принципы и практика, четвертое издание" А. Э. Роя и Д. Кларка (2003); моя копия - первое издание 1977 года, Адам Хилгер, ISBN 0-85274-346-7.)
NB Check out (Google) 'define:"морская миля"'; Похоже, что морская миля в настоящее время составляет 1852 м (1,852 км) по определению. Множитель 1.1515 соответствует старому определению морской мили как приблизительно 6080 футов. Использование bc со шкалой 10 я получаю:
(1852/(3*0.3048))/1760
1.1507794480
Какой фактор работает для вас, зависит от того, какова ваша основа.
Рассматривая вторую проблему из первых принципов, мы имеем немного иную схему, и нам нужно "другое" сферическое тригонометрическое уравнение, формула синуса:
sin A sin B sin C
----- = ----- = -----
sin a sin b sin c
Адаптация предыдущей диаграммы:
+ C
/|
/ |
a / | b
| / |
|X/ |
|/ |
B +------+ A
c
Вам задается начальная точка B, угол X = 90º - B, длина (угол) a и угол A = 90 °. Что вам нужно, так это b (дельта по широте) и c (дельта по долготе).
Итак, имеем:
sin a sin b
----- = ----
sin A sin B
Или же
sin a . sin B
sin b = -------------
sin A
Или, поскольку A = 90°, sin A = 1, а sin B = sin (90° - X) = cos X:
sin b = sin a . cos X
Это означает, что вы преобразуете пройденное расстояние в угол a, берете синус этого, умножаете на косинус направления курса и берете арксинус результата.
Учитывая a, b (только что вычисленные) и A и B, мы можем применить формулу косинуса, чтобы получить c. Обратите внимание, что мы не можем просто повторно применить формулу синуса, чтобы получить c, поскольку у нас нет значения C, и, поскольку мы играем со сферической тригонометрией, нет удобного правила, что C = 90° - B (сумма из углов в сферическом треугольнике может быть больше 180° (рассмотрим равносторонний сферический треугольник со всеми углами, равными 90°, что вполне выполнимо).
Проверьте http://www.movable-type.co.uk/scripts/latlong.html
На этом сайте есть много разных формул и кода Javascript, которые должны вам помочь. Я успешно перевел его как на C#, так и на UDF SQL Server, и я использую их повсеместно.
Например, для Javascript для расчета расстояния:
var R = 6371; // km
var φ1 = lat1.toRadians();
var φ2 = lat2.toRadians();
var Δφ = (lat2-lat1).toRadians();
var Δλ = (lon2-lon1).toRadians();
var a = Math.sin(Δφ/2) * Math.sin(Δφ/2) +
Math.cos(φ1) * Math.cos(φ2) *
Math.sin(Δλ/2) * Math.sin(Δλ/2);
var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
var d = R * c;
Наслаждайтесь!
Ваше преобразование между км и радианами неверно. Морская миля составляет 1/60 градуса, поэтому при условии, что 1,15... это ваше преобразование из миль в морские мили, а 1,6... это ваше преобразование из км в статутные мили,
nm = km / (1.1515 * 1.609344);
deg = nm / 60;
rad = toRadians(deg);
Другими словами, я думаю, что вы в 1000 раз.
Относительно вашего обновленного вопроса: не должен
double lon1 = endOfCourse(lat1, 0.0, 90.0, distEast)[0];
быть
double lon1 = endOfCourse(lat1, 0.0, 90.0, distEast)[1];
Я понял большую проблему с этими формулами, кроме ошибок реализации, упомянутых в других ответах и обновлениях.
Большая проблема заключалась в следующем: метод Distance (для вычисления расстояния между двумя точками) вычислял расстояния большого круга. Что, конечно, имеет смысл - это кратчайший путь между двумя точками. Однако расстояние по большому кругу между двумя точками, которые лежат на одной параллели (линия широты), НЕ совпадает с расстоянием между этими двумя точками при движении прямо вдоль линии широты, если только вы не на экваторе.
Итак: функции работают правильно; однако альтернативная система координат, которую я предложил в первоначальном вопросе, требует, чтобы мы смотрели только на расстояние север вдоль IDL, а затем расстояние восток вдоль параллели на результирующей широте. И вычисление расстояния по определенной параллели сильно отличается от вычисления расстояния по большому кругу!
Во всяком случае, у вас есть это.