Я пытаюсь выяснить, как математически вывести общие точки двух пересекающихся кругов на земной поверхности с учетом центра широты и долготы и радиуса для каждой точки.
Например, учитывая:
- Широта / долгота (37,673442, -90,234036) радиус 107,5 нм
- Шир. / Долг. (36.109997, -90.953669) Радиус 145 НМ
Я должен найти две точки пересечения с одной из них (36.948, -088.158).
Было бы тривиально легко решить это на плоской плоскости, но у меня нет опыта решения уравнений на несовершенной сфере, такой как поверхность Земли.
Ответы:
Это не намного сложнее в сфере, чем в самолете, как только вы поймете, что
Рассматриваемые точки - это взаимные пересечения трех сфер: сфера с центром под местоположением x1 (на поверхности земли) данного радиуса, сфера с центром под местоположением x2 (на поверхности земли) данного радиуса и сама земля , которая представляет собой сферу с центром в O = (0,0,0) данного радиуса.
Пересечение каждой из первых двух сфер с поверхностью земли представляет собой круг, который определяет две плоскости. Таким образом, взаимные пересечения всех трех сфер лежат на пересечении этих двух плоскостей: линии .
Следовательно, задача сводится к пересечению линии со сферой, что легко.
Вот подробности. Входными данными являются точки P1 = (lat1, lon1) и P2 = (lat2, lon2) на поверхности земли, рассматриваемые как сфера, и два соответствующих радиуса r1 и r2.
Преобразовать (широта, долгота) в (x, y, z) геоцентрические координаты. Как обычно, потому что мы можем выбрать единицы измерения, в которых Земля имеет единичный радиус,
В этом примере P1 = (-90,234036 градусов, 37,673442 градусов) имеет геоцентрические координаты x1 = (-0,00323306, -0,7915, 0,61116) и P2 = (-90,953669 градусов, 36,109997 градусов) имеет геоцентрические координаты x2 = (-0,0134475, -0. 0,589337).
Преобразуйте радиусы r1 и r2 (которые измеряются вдоль сферы) в углы вдоль сферы. По определению, одна морская миля (NM) составляет 1/60 градуса дуги (что равно pi / 180 * 1/60 = 0,0002908888 радиан). Поэтому, как углы,
Геодезический круг радиуса r1 вокруг x1 является пересечением поверхности Земли с евклидовой сферы радиуса греха (r1) с центром в Cos (r1) * x1.
Плоскость, определяемая пересечением сферы радиуса sin (r1) вокруг cos (r1) * x1 и земной поверхности, перпендикулярна x1 и проходит через точку cos (r1) x1, откуда ее уравнение равно x.x1 = cos (r1) («.» представляет обычное скалярное произведение ); аналогично для другого самолета. На пересечении этих двух плоскостей будет уникальная точка x0, представляющая собой линейную комбинацию x1 и x2. Запись x0 = a x1 + b * x2, два плоских уравнения
Используя тот факт, что x2.x1 = x1.x2, который я напишу как q, решение (если оно существует) дается
В текущем примере я вычисляю a = 0,973503 и b = 0,0260194.
Очевидно, нам нужно q ^ 2! = 1. Это означает, что x1 и x2 не могут быть ни одной, ни антиподальной точкой.
Теперь все остальные точки на линии пересечения двух плоскостей отличаются от x0 некоторым кратным вектора n, взаимно перпендикулярного обеим плоскостям. Крестовый продукт
выполняет задание, если n отлично от нуля: еще раз, это означает, что x1 и x2 не совпадают и не диаметрально противоположны. (Мы должны позаботиться о том, чтобы вычислить перекрестное произведение с высокой точностью, поскольку оно включает в себя вычитания с большим количеством отмен, когда x1 и x2 находятся близко друг к другу.) В этом примере n = (0.0272194, -0.00631254, -0.00803124) ,
Поэтому мы ищем до двух точек вида x0 + t * n, которые лежат на поверхности земли: их длина равна 1. Эквивалентно, их квадратная длина равна 1:
Член с x0.n исчезает, потому что x0 (будучи линейной комбинацией x1 и x2) перпендикулярен n. Два решения легко
и его негатив. Еще раз требуется высокая точность, потому что когда x1 и x2 близки, x0.x0 очень близко к 1, что приводит к некоторой потере точности с плавающей запятой. В этом примере t = 1,07509 или t = -1,07509. Поэтому две точки пересечения равны
Наконец, мы можем преобразовать эти решения обратно в (широта, долгота), преобразовав геоцентрические (x, y, z) в географические координаты:
Для долготы используйте обобщенный арктангенс, возвращающий значения в диапазоне от -180 до 180 градусов (в вычислительных приложениях эта функция принимает в качестве аргументов как x, так и y, а не только отношение y / x; иногда ее называют «ATan2»).
Я получаю два решения (-88.151426, 36.989311) и (-92.390485, 38.238380), показанные на рисунке желтыми точками.
На осях отображаются геоцентрические (x, y, z) координаты. Серый участок - это часть земной поверхности от -95 до -87 градусов по долготе, от 33 до 40 градусов по широте (отмечена одной градусной сеткой). Земная поверхность была сделана частично прозрачной, чтобы показать все три сферы. Правильность вычисленных решений очевидна тем, как желтые точки находятся на пересечениях сфер.
источник
Эллипсоидальной случай:
Эта проблема является обобщением проблемы нахождения морских границ, определенных как «срединные линии», и по этой теме имеется обширная литература. Мое решение этой проблемы заключается в использовании эквидистантной азимутальной проекции:
Этот алгоритм сходится квадратично и дает точное решение на эллипсоиде. (Точность требуется в случае морских границ, поскольку она определяет права на рыболовство, добычу нефти и полезных ископаемых.)
Формулы приведены в разделе 14 геодезических на эллипсоиде вращения . Эллипсоидальная эквидистантная азимутальная проекция предоставлена GeographicLib . Версия MATLAB доступна в геодезических проекциях для эллипсоида .
источник
Вот некоторый код R, чтобы сделать это:
источник
Исходя из ответа @ whuber , вот код Java, который полезен по двум причинам:
Он не оптимизирован и не завершен (я пропустил такие очевидные классы, как
Point
), но должен сделать свое дело.Также, что важно, обратите внимание на использование
atan2
- это обратное тому, что вы ожидаете от ответа @ whuber (я не знаю почему, но это работает):источник
Рабочий код 'R' для ответа @wuhber.
источник
Если один из кругов - это Нортстар, то с юнит-сферой есть самый простой способ.
Вы можете измерить свою широту с Nortstar. Тогда у вас есть относительная позиция в этой сфере. v1 (0, sin (la), cos (la)) Вы знаете положение (угол) другой звезды (star2) из альманаха. v2 (sin (lo2) * cos (la2), sin (la2), cos (lo2) * cos (la2)) Его векторы. Из уравнения сферы.
lo2 - относительная долгота. Это неизвестно .
Угол между вами и звездой2 вы тоже можете измерить, (м) И вы знаете, что внутренним продуктом двух единичного вектора является cos (угол) между ними. cos (m) = точка (v1, v2) теперь вы можете рассчитать относительную долготу (lo2). lo2 = экос ((соз (м) -sin (ли) * sin (ла2)) / (соз (ли) * соз (ла2)))
В конце концов вы добавляете настоящую долготу звезды2 к lo2. (или sub, зависит от его на западной стороне от вас, или на востоке.) lo2 теперь ваша долгота.
Извините за мой английский, я никогда не изучаю этот язык.
2 вещи: Северная звезда означает Полярную звезду.
Другая. Поскольку угол измеряется относительно горизонтали, всегда требуется коррекция на 90 градусов. Это справедливо и для угла м.
PS: реальный угол означает: положение звезды - время коррекции.
источник