Например, у меня есть координаты для трех базовых точек на береговой линии, и мне нужно найти координаты точки у берега, которая равноудалена от всех трех. Это простое упражнение по геометрии, но все измерения должны учитывать геодезию.
Если бы я подходил к этому евклидовым образом, я мог бы измерить геодезические пути, соединяющие базовые точки, найти середины сторон полученного треугольника и создать перпендикулярные ортодромы для каждого из этих путей. Предполагается, что три локсодрома сходятся в равноудаленной точке. Если это правильный метод, должен быть более простой способ сделать это в Arc.
Ответы:
Этот ответ состоит из нескольких разделов:
Анализ и уменьшение проблемы , показывая, как найти нужную точку с помощью «стандартных» процедур.
Иллюстрация: рабочий прототип , дающий рабочий код.
Пример , показывающий примеры решений.
Подводные камни , обсуждение потенциальных проблем и способы их решения.
Реализация ArcGIS , комментарии о создании специального инструмента ArcGIS и о том, где можно получить необходимые подпрограммы.
Анализ и уменьшение проблемы
Давайте начнем с наблюдения, что в (идеально круглой) сферической модели всегда найдется решение - фактически два решения. Для заданных базовых точек A, B и C каждая пара определяет свой «перпендикулярный биссектрис», который представляет собой набор точек, равноудаленных от двух данных точек. Этот биссектриса является геодезической (большой круг). Сферическая геометрия эллиптическая : любые две геодезические пересекаются (в двух уникальных точках). Таким образом, точки пересечения биссектрисы AB и биссектрисы BC по определению равноудалены от A, B и C, что решает проблему. (См. Первый рисунок ниже.)
На эллипсоиде все выглядит сложнее, но поскольку это небольшое возмущение сферы, мы можем ожидать аналогичного поведения. (Анализ этого приведет нас слишком далеко.) Сложные формулы, используемые (внутри ГИС) для вычисления точных расстояний на эллипсоиде, не являются концептуальным осложнением: проблема в основном та же. Чтобы увидеть, насколько на самом деле проста проблема, давайте сформулируем ее несколько абстрактно. В этом утверждении «d (U, V)» относится к истинному, полностью точному расстоянию между точками U и V.
Эти три расстояния все зависит от неизвестного X . Таким образом, различия в расстояниях u (X) = d (X, A) - d (X, B) и v (X) = d (X, B) - d (X, C) являются действительными функциями X. Опять же, несколько абстрактно, мы можем собрать эти различия в упорядоченную пару. Мы также будем использовать (lat, lon) в качестве координат для X, что позволит нам рассматривать его как упорядоченную пару, скажем, X = (phi, lambda). В этой настройке функция
является функцией от части двумерного пространства, принимающей значения в двумерном пространстве, и наша задача сводится к
Вот где абстракция окупается: существует множество отличного программного обеспечения для решения этой (чисто численной многомерной задачи поиска корней). Это работает так, что вы пишете подпрограмму для вычисления F , а затем передаете ее программному обеспечению вместе с любой информацией об ограничениях на ее вход ( phi должно лежать в диапазоне от -90 до 90 градусов, а лямбда должна лежать в диапазоне от -180 до 180 градусов). Он запускается за доли секунды и возвращает (обычно) только одно значение ( фи , лямбда ), если он может найти его.
Есть детали для обработки, потому что в этом есть искусство: есть разные методы решения на выбор, в зависимости от того, как F «ведет себя»; это помогает «управлять» программным обеспечением, предоставляя ему разумную отправную точку для его поиска (это один из способов получить ближайшее решение, а не любое другое); и вам обычно нужно указать, насколько точным вы хотите, чтобы решение было (чтобы он знал, когда прекратить поиск). (Для получения дополнительной информации о том, что аналитики ГИС должны знать о таких деталях, которые часто встречаются в проблемах ГИС, посетите разделы «Рекомендации», которые необходимо включить в курс «Компьютерные науки для геопространственных технологий») и посмотрите в конце раздела «Разное». )
Иллюстрация: рабочий прототип
Анализ показывает нам нужно запрограммировать две вещи: неочищенную первоначальную оценку решения и вычисление F сам.
Первоначальная оценка может быть сделана «сферическим средним» трех базовых точек. Это получается путем представления их в геоцентрических декартовых (x, y, z) координатах, усреднения этих координат и проецирования этого среднего обратно в сферу и повторного ее выражения по широте и долготе. Размер сферы не имеет значения, и поэтому расчеты производятся просто: поскольку это только отправная точка, нам не нужны эллипсоидальные вычисления.
Для этого рабочего прототипа я использовал Mathematica 8.
(Конечное
Ifусловие проверяет, может ли среднее значение явно не указывать долготу; если это так, оно возвращается к прямому арифметическому среднему значению широт и долгот своего входа - возможно, это не лучший выбор, но, по крайней мере, правильный. Для тех, кто использует этот код для руководства по реализации, обратите внимание, что аргументы MathematicaArcTanобращены вспять по сравнению с большинством других реализаций: его первый аргумент - это координата x, второй - координата y, и он возвращает угол, заданный вектором ( х, у).)Что касается второй части, поскольку Mathematica - как ArcGIS и почти все другие ГИС - содержит код для вычисления точных расстояний на эллипсоиде, писать практически нечего. Мы просто вызываем процедуру поиска корня:
Наиболее примечательным аспектом этой реализации является то, как она избегает необходимости ограничивать широту (
f) и долготу (q), всегда вычисляя их по модулю 180 и 360 градусов соответственно. Это позволяет избежать необходимости ограничивать проблему (что часто создает осложнения). Параметры управленияMaxIterationsи т. Д. Настроены так, чтобы этот код обеспечивал максимально возможную точность.Чтобы увидеть это в действии, давайте применим его к трем базовым пунктам, приведенным в связанном вопросе :
Вычисленные расстояния между этим решением и тремя точками
(это метры). Они соглашаются через одиннадцатую значащую цифру (которая на самом деле является слишком точной, поскольку расстояния редко бывают точными с точностью до миллиметра или около того). Вот изображение этих трех точек (черный), их три взаимных биссектрисы и решение (красный):
пример
Чтобы протестировать эту реализацию и лучше понять, как ведет себя проблема, вот контурный график среднеквадратичного расхождения расстояний для трех широко расположенных базовых точек. (Среднеквадратичное расхождение получается путем вычисления всех трех разностей d (X, A) -d (X, B), d (X, B) -d (X, C) и d (X, C) -d (X). , A), усредняя их квадраты и беря квадратный корень.Это равно нулю, когда X решает проблему, и в противном случае возрастает, когда X удаляется от решения, и таким образом измеряет, насколько «мы близки» к тому, чтобы быть решением в любом месте. )
Базовые точки (60, -120), (10, -40) и (45,10) показаны красным цветом в этой проекции Плато Карри; решение (49.2644488, -49.9052992), для вычисления которого потребовалось 0,03 секунды, выделено желтым цветом. Его среднеквадратичное отклонение составляет менее трех нанометров , несмотря на то, что все соответствующие расстояния составляют тысячи километров. Темные области показывают небольшие значения RMS, а светлые области показывают высокие значения.
Эта карта ясно показывает, что рядом находится другое решение (-49.2018206, 130.0297177) (вычислено для среднеквадратичного значения в два нанометра путем установки начального значения поиска, диаметрально противоположного первому решению.)
Ловушки
Численная нестабильность
Когда базовые точки почти коллинеарны и близки друг к другу, все решения будут располагаться почти на полмира, и их будет очень сложно точно определить. Причина в том, что небольшие изменения в местоположении по всему земному шару - перемещение его к базовым точкам или от них - вызовет только невероятно малые изменения в разнице расстояний. Просто нет точности и точности, встроенной в обычный расчет геодезических расстояний, чтобы определить результат.
Например, начиная с базовых точек в (45.001, 0), (45, 0) и (44.999,0), которые разделены вдоль простого меридиана только на 111 метров между каждой парой,
triполучают решение (11.8213, 77.745) ). Расстояния от него до базовых точек составляют 8 127 964,998; 77; 8 127 964,998; 41; и 8 127 964,998 65 метров соответственно. Они согласны с точностью до миллиметра! Я не уверен, насколько точным может быть этот результат, но не удивлюсь, если другие реализации вернут местоположения далеко от этого, демонстрируя почти такое же хорошее равенство трех расстояний.Время вычислений
Эти вычисления, поскольку они включают в себя значительный поиск с использованием сложных вычислений расстояния, не являются быстрыми и обычно требуют заметной доли секунды. Приложения реального времени должны знать об этом.
Реализация ArcGIS
Python является предпочтительной средой сценариев для ArcGIS (начиная с версии 9). В пакете scipy.optimize есть многомерный rootfinder,
rootкоторый должен делать то, чтоFindRootделает в коде Mathematica . Конечно, сама ArcGIS предлагает точные расчеты эллипсоидального расстояния. Остальное - это все детали реализации: решить, как будут получены координаты базовой точки (из слоя, набранного пользователем, из текстового файла, из мыши) и как будет представлен результат (в виде координат). отображается на экране (в виде графической точки, как новый точечный объект в слое), напишите этот интерфейс, перенесите код Mathematica, показанный здесь (прямо), и все будет готово.источник
Как вы заметили, эта проблема возникает при определении морских границ; ее часто называют «трехточечной» проблемой, и вы можете найти ее в Google и найти несколько статей, посвященных ее решению. Одна из этих статей написана мной (!), И я предлагаю точное и быстро сходящееся решение. См. Раздел 14 http://arxiv.org/abs/1102.1215.
Метод состоит из следующих шагов:
Необходимая формула для трехточечного решения в проекции приведена в статье. Пока вы используете точную азимутальную равноудаленную проекцию, ответ будет точным. Сходимость является квадратичной, что означает, что требуется всего несколько итераций. Это почти наверняка превзойдет общие методы поиска корней, предложенные @whuber.
Я не могу напрямую помочь вам с ArcGIS. Вы можете взять мой пакет python для выполнения геодезических расчетов с https://pypi.python.org/pypi/geographiclib, и написание кода на основе этого просто.
редактировать
Проблема нахождения три-точки в вырожденном случае @uuber (45 + eps, 0) (45,0) (45-eps, 0) была рассмотрена Кэли в статье « О геодезических линиях на сплюснутом сфероиде» , Phil. Магнето (1870), http://books.google.com/books?id=4XGIOoCMYYAC&pg=PA15
В этом случае тройную точку получают, следуя геодезической из (45,0) с азимутом 90 и находя точку, в которой исчезает геодезическая шкала. Для эллипсоида WGS84 эта точка равна (-0,10690908732248, 89,89291072793167). Расстояние от этой точки до каждого из (45,001,0), (45,0), (44,999) составляет 10010287,665788943 м (в пределах нанометра или около того). Это примерно на 1882 км больше, чем оценка Whuber (которая показывает, насколько нестабильным является этот случай). Для сферической земли трипункт будет, конечно, (0,90) или (0, -90).
ДОБАВЛЕНИЕ: Вот реализация азимутального эквидистантного метода с использованием Matlab
Тестирование это с помощью Octave я получаю
как три точки для Нью-Йорка, Токио и Веллингтона.
Этот метод является неточным для соседних коллинеарных точек, например, [45,001,0], [45,0], [44,999,0]. В этом случае вы должны решить для M 12 = 0 для геодезической, исходящей из [45,0] при азимуте 90. Следующая функция добивается цели (используя метод Ньютона):
Для примера это дает:
источник
Мне было любопытно посмотреть, как быстро подход @ cffk сходится к решению, поэтому я написал тест с использованием arcobjects, который вывел этот вывод. Расстояния в метрах:
Вот исходный код. (Правка) Изменен FindCircleCenter для обработки пересечений (центральных точек), которые падают с края азимутальной проекции:
Есть также альтернативный подход в выпуске журнала MSDN Magazine в июне 2013 года « Оптимизация метода Amoeba с использованием C #» .
редактировать
Ранее опубликованный код в некоторых случаях сходился к антиподу. Я изменил код так, чтобы он выводил этот результат для контрольных точек @ cffk.
Вот результат, который он теперь производит:
Вот пересмотренный код:
редактировать
Вот результаты, которые я получаю с esriSRProjCS_WGS1984N_PoleAziEqui
источник