Геопространственные координаты и расстояние в километрах
https://stackoverflow.com/questions/389211
-
23-08-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianВопрос
Это продолжение этот вопрос.
Кажется, я застрял на этом.По сути, мне нужно иметь возможность конвертировать взад и вперед в ссылку на координаты либо в стандартной системе градусов, либо путем измерения расстояния к северу от южного полюса вдоль международной линии даты, а затем расстояния на восток, начиная с этой точки на дате. линия.Для этого (а также для некоторых более общих задач по измерению расстояний) у меня есть один метод для определения расстояния между двумя точками широты и долготы, а также другой метод, который принимает точку широты и долготы, курс и расстояние и возвращает точка широты и долготы в конце этого курса.
Вот два статических метода, которые я определил:
/* 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;
}
Значения «конец в» должны составлять ок.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;
Решение
У вас серьезный случай магических чисел в коде.Выражение:
(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)
должно быть одинаковым, но я почти уверен, что ваш код дает разные ответы.
Итак, я думаю, вам нужно вернуться к справочным материалам по сферической тригонометрии и посмотреть, что вы делаете неправильно.(Мне потребуется некоторое время, чтобы найти свою книгу по этой теме — ее нужно будет распаковать из той коробки, в которой она находится.)
[...время идет...распаковка завершена...]
Дан сферический треугольник с углами А, Б, С по вершинам и сторонам а, б, с противоположные этим вершинам (то есть сторона а из Б к С, и т. д.), формула косинуса:
cos a = cos b . cos c + sin b . sin c . cos A
Применяя это к проблеме, мы можем назвать две данные точки Б и С, и создадим прямоугольный сферический треугольник с прямым углом при А.
ASCII-искусство в худшем виде:
+ C
/|
/ |
a / | b
/ |
/ |
/ |
B +------+ A
c
Сторона с равен разнице долгот;сторона б равен разнице широт;угол А равен 90°, так что А = 0.Поэтому я считаю, что уравнение для а является:
cos a = cos Δλ . cos Δφ + sin Δλ . sin Δφ . cos 90°
a = arccos (cos Δλ . cos Δφ)
Угол а затем в радианах преобразуется в расстояние путем умножения на радиус Земли.Альтернативно, учитывая а в градусах (и долях градуса), то в одном градусе 60 морских миль, следовательно, 60*1,1515 статутных миль, и 60*1,1515*1,609344 километров в одном градусе.Если вам не нужно расстояние в метрах, я не вижу необходимости в коэффициенте 1000.
Пол Томблин указывает на Авиационный формуляр v1.44 в качестве источника уравнения - и действительно, оно существует вместе с более устойчивой в цифровом отношении версией для случаев, когда разница в положении невелика.
Переходя к базовой тригонометрии, мы также знаем, что:
cos (A - B) = cos A . cos B + sin A . sin B
Применив это дважды в приведенном мною уравнении, вполне можно получить формулу из Авиационного формуляра.
(Моя ссылка: «Астрономия:Принципы и практика, четвертое издание» А. Э. Роя и Д. Кларка (2003);мой экземпляр — первое издание 1977 года, Адам Хилгер, ISBN 0-85274-346-7.)
Примечание: Проверьте (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
Вам дана отправная точка Б, угол Икс = 90º - B, длина (угол) а, и угол А = 90°.То, что вам нужно, это б (дельта широты) и с (дельта долготы).
Итак, у нас есть:
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
Это означает, что вы конвертируете пройденное расстояние в угол. а, возьмите его синус, умножьте на косинус направления курса и возьмите арксинус результата.
Данный а, б (только что рассчитано) и А и Б, мы можем применить формулу косинуса, чтобы получить с.Обратите внимание, что мы не можем просто повторно применить формулу синуса, чтобы получить с поскольку у нас нет значения С и, поскольку мы играем со сферической тригонометрией, не существует удобного правила, согласно которому 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];
Я выяснил большую проблему с этими формулами, не считая ошибок реализации, упомянутых в других ответах и обновлениях.
Большая проблема заключалась в следующем:Метод «Расстояние» (для вычисления расстояния между двумя точками) заключался в вычислении расстояний по большому кругу.Что, конечно, имеет смысл — это кратчайший путь между двумя точками. Однако, расстояние по большому кругу между двумя точками, лежащими на одной параллели (линии широты), НЕ совпадает с расстоянием между этими двумя точками при путешествии прямо вдоль линии широты, если только вы не находитесь на экваторе.
Так:функции работают правильно;однако альтернативная система координат, которую я предложил в исходном вопросе, требует, чтобы мы смотрели только на расстояние на север вдоль IDL, за которым следует расстояние на восток вдоль параллели на полученной широте.А вычисление расстояния по определенной параллели сильно отличается от расчета расстояния по большому кругу!
В любом случае, вот оно.
