Точность сферической тригонометрии на эллипсоиде

[Игорь Белов]

желательно что-нибудь подходящее для Delphi / C++ Builder / MS Visual Studio.

Там есть всё, чего можно желать.

[KotAlex]

версия 2015 MS Visual Studio нормально работает с Питоном (в том числе бесплатная). Язык, конечно мутноватый (хотя и не сильно мутнее Матлаба/R), но, поскольку нужны только вызовы подпрограмм, продраться можно, да и примеров много

С Питоном я уже пытался разобраться где-то пару лет назад, когда пытался скрестить QGIS с MS Visual Studio. Уже не помню, с какой такой неожиданной стороны выплыл тогда Питон, но от разбирательства с ним, насколько помню, у меня остались весьма тягостные впечатления. Результат, откровенно говоря, от всех этих разбирательств остался такой - предпочел бы с Питоном не связываться. Да и других проблем хватает по горло ...

[KotAlex]

Там есть всё, чего можно желать.

Спасибо за ссылку. По сути - это то же самое, что предложил и Эдуард Казаков. Однако, прежде чем устраивать большую разборку в дебрях исходных кодов этой библиотеки, решил разобраться с другим, что полегче.
Бегло просматривая статью
http://gis-lab.info/qa/sphere-geodesic- ... oblem.html
обнаружил в ней небольшой раздел с примером использования утилиты geod. Решил поиграть с ней.
Однако, документация по ней представляет весьма небольшую по объему страницу и осталось много непонятного.
Помогите разобраться.
Я запустил пример в этой статье :
geod +a=6371000 -f "%f" +units=km dir.dat
с таким же файлом dir.dat. Сработало прекрасно, результаты получились точь-в-точь те же самые, что и в статье.
Не могу понять значения параметров в командной строке.
1. Как я понял - параметры с предваряющим минусом - это аргументы, описанные в документации
https://trac.osgeo.org/proj/wiki/man_geod
а с плюсом - стандартные аргументы PROJ4 ? так ?
2. Чем определяется - когда исходные данные прямой задачи ( координаты точки, азимут, дистанция ) указываются в файле ( dir.dat ), а когда - в командной строке ? Или вариант с указанием данных в командной строке не работает ?
3. Как указать аргументы для двух вариантов расчета - для сферы и для эллипсоида WGS84 ?

[Ariki]

1. Да. Аргументы PROJ.4 задают параметры эллипсоида.
2. Программа читает данные из указанных файлов. Если файл не указан, читает из стандартного потока ввода. Чтобы указать данные в командной строке, можно использовать пайплайн:

Код: Выделить всё

echo "53 84 45 1000" | geod +ellps=WGS84

Под Linux точно работает; под Windows тоже должно, но не проверял. Также на странице руководства есть пример для Linux с использованием heredoc (<<).
3. Пример для эллипсоида WGS84 — см. п. 2. Пример для сферы заданного радиуса (параметр +a) из статьи Игоря Белова вы процитировали выше. Параметр -f задаёт формат вывода значений координат (как в функции printf() языка C).

[KotAlex]

...

Большое спасибо за ответы. Я тоже именно так и понял из документации. Спрашиваю потому, что мне нужна стопроцентная уверенность.
Я ключом -le :

Код: Выделить всё

 geod -le

вывел список возможных значений параметра +ellps= ... , там в самом низу присутствует значение "sphere" , Запустил :

Код: Выделить всё

geod +ellps=sphere ...

Проходит нормально. Как я понимаю, это значение параметра +ellps= эквивалентно запуску с параметром

Код: Выделить всё

geod +a=6371000 ...

( то есть, эти параметры идентичны друг другу ) ?

[Ariki]

Не совсем. У меня выводит такие параметры для sphere:
sphere a=6370997.0 b=6370997.0 Normal Sphere (r=6370997)
Лучше задать +a и +b явно.

[KotAlex]

Не совсем. У меня выводит такие параметры для sphere:
sphere a=6370997.0 b=6370997.0 Normal Sphere (r=6370997)
Лучше задать +a и +b явно.

Действительно, те значения для a и b, которые у вас ( 6370997.0 ), отличаются от того, которое у меня ( +a=6371000 ).
Хорошо, уточню вопрос - если я задам, как у вас :

Код: Выделить всё

geod +a=6370997.0 ...

то это будет полностью эквивалентно заданию :

Код: Выделить всё

geod +ellps=sphere ...

?
Параметры a и b, как я понял - большая и малая полуоси ?
Спрашиваю потому, что предпочитаю задавать параметры в виде строковых констант, а не числовых, с которыми еще напутать можно.

[Ariki]

Я так и не понял, какие параметры у вас выводил geod -le для sphere: такие же, как у меня, или те, что вы указываете. В моём представлении, они должны быть постоянными и не меняться от версии к версии, но кто его знает. Эллипсоид WGS84 стандартизирован и применяется повсеместно, так что его логично будет указывать по имени, чтобы нигде не опечататься в значениях длины полуосей. А вот сфер всяких в разных программных продуктах используется куча, и я бы предпочёл явно указывать радиус, чтобы с первого взгляда на команду было ясно, на какой сфере выполняется расчёт. А так это просто вопрос удобства: указание предопределённого эллипсоида эквивалентно ручному заданию его параметров, которые можно узнать из вывода geod -le

[KotAlex]

У меня geod -le выводит следующее :

Код: Выделить всё

  MERIT a=6378137.0      rf=298.257       MERIT 1983
    SGS85 a=6378136.0      rf=298.257       Soviet Geodetic System 85
    GRS80 a=6378137.0      rf=298.257222101 GRS 1980(IUGG, 1980)
    IAU76 a=6378140.0      rf=298.257       IAU 1976
     airy a=6377563.396    b=6356256.910    Airy 1830
   APL4.9 a=6378137.0.     rf=298.25        Appl. Physics. 1965
    NWL9D a=6378145.0.     rf=298.25        Naval Weapons Lab., 1965
 mod_airy a=6377340.189    b=6356034.446    Modified Airy
   andrae a=6377104.43     rf=300.0         Andrae 1876 (Den., Iclnd.)
  aust_SA a=6378160.0      rf=298.25        Australian Natl & S. Amer. 1969
    GRS67 a=6378160.0      rf=298.2471674270 GRS 67(IUGG 1967)
   bessel a=6377397.155    rf=299.1528128   Bessel 1841
 bess_nam a=6377483.865    rf=299.1528128   Bessel 1841 (Namibia)
   clrk66 a=6378206.4      b=6356583.8      Clarke 1866
   clrk80 a=6378249.145    rf=293.4663      Clarke 1880 mod.
      CPM a=6375738.7      rf=334.29        Comm. des Poids et Mesures 1799
   delmbr a=6376428.       rf=311.5         Delambre 1810 (Belgium)
  engelis a=6378136.05     rf=298.2566      Engelis 1985
  evrst30 a=6377276.345    rf=300.8017      Everest 1830
  evrst48 a=6377304.063    rf=300.8017      Everest 1948
  evrst56 a=6377301.243    rf=300.8017      Everest 1956
  evrst69 a=6377295.664    rf=300.8017      Everest 1969
  evrstSS a=6377298.556    rf=300.8017      Everest (Sabah & Sarawak)
  fschr60 a=6378166.       rf=298.3         Fischer (Mercury Datum) 1960
 fschr60m a=6378155.       rf=298.3         Modified Fischer 1960
  fschr68 a=6378150.       rf=298.3         Fischer 1968
  helmert a=6378200.       rf=298.3         Helmert 1906
    hough a=6378270.0      rf=297.          Hough
     intl a=6378388.0      rf=297.          International 1909 (Hayford)
    krass a=6378245.0      rf=298.3         Krassovsky, 1942
    kaula a=6378163.       rf=298.24        Kaula 1961
    lerch a=6378139.       rf=298.257       Lerch 1979
    mprts a=6397300.       rf=191.          Maupertius 1738
 new_intl a=6378157.5      b=6356772.2      New International 1967
  plessis a=6376523.       b=6355863.       Plessis 1817 (France)
   SEasia a=6378155.0      b=6356773.3205   Southeast Asia
  walbeck a=6376896.0      b=6355834.8467   Walbeck
    WGS60 a=6378165.0      rf=298.3         WGS 60
    WGS66 a=6378145.0      rf=298.25        WGS 66
    WGS72 a=6378135.0      rf=298.26        WGS 72
    WGS84 a=6378137.0      rf=298.257223563 WGS 84
   sphere a=6370997.0      b=6370997.0      Normal Sphere (r=6370997)

Как видите, здесь радиус сферы - 6379997.
У меня библиотека сферической геометрии работает при значении радиуса - 6370000.
Как такое могло получиться - ума не приложу ... Переправить, конечно, нетрудно, но - все-таки ...
А если еще учесть ваше :

А вот сфер всяких в разных программных продуктах используется куча

то получается, что шар земной вообще резиновый ...
М - да - а - а, без поллитры не разберешь ...

[stout]

то получается, что шар земной вообще резиновый ...
М - да - а - а, без поллитры не разберешь ...

Дан эллипсоид. Найти радиус эквивалентной сферы.
Ответ зависит от того, что понимать под эквивалентной сферой. Это сфера той же площади или того же объёма?
Или сфера с радиусом, равным среднему радиусу кривизны эллипсоида в точке?
Разные условия дают разный ответ. Насколько помню, на практике различают 5 разных радиусов. Для каждой задачи – свой радиус.
Сайт GeographicLib хорош ещё и тем, что там есть библиография, начиная от Ньютона и до наших дней. Практически всё в открытом доступе.

[Alexa111]

Камрады, прошу помощи.
На картах мелких масштабов ( карты стран, регионов и мира ) решаются задачи типа прямой и обратной геодезических, и ведутся разного рода геометрические построения. С целью упрощения расчетов все построения ведутся по методам сферической тригонометрии.
Интересует вопрос - насколько велики будут ошибки и искажения вследствие применения методов сферической тригонометрии для эллиптической поверхности ? имеются ли какие-либо численные выражения подобных ошибок ?
Для максимальной определенности максимально конкретизирую задачу : строится линия кратчайшего расстояния ( ортодрома ) между двумя точками земной поверхности, разнесенными друг от друга на расстояние, скажем, в четверть земной дуги. Насколько велико тут Читы будет отклонение этой линии, построенной на сфере, от этой же линии, построенной на эллипсоиде ?
Или другая задача - для заданной точки найти координаты точки, отстоящей от первой на заданных азимуте и дистанции ( дистанция, как и в первом случае - четверть земной окружности ). Насколько велики расхождения вторых точек в случае расчетов по указанным методам ?
Я знаю, что на эллипсоиде эта линия строится методом численного интегрирования дифференциальных
уравнений геодезической линии, и на основе этого можно получить численные оценки расхождений.
Но меня интересует - не имеются ли уже готовые ответы на этот вопрос ?

Тоже интересует данный вопрос. Как можно с вами связаться? Оставьте пожалуйста какой-нибудь контакт! Не могу послать личное сообщение..

[KotAlex]

Прошу прощения за долгое молчание. В последние недели не было возможности заглядывать сюда.

Дан эллипсоид. Найти радиус эквивалентной сферы.
Ответ зависит от того, что понимать под эквивалентной сферой. Это сфера той же площади или того же объёма?
Или сфера с радиусом, равным среднему радиусу кривизны эллипсоида в точке?
Разные условия дают разный ответ. Насколько помню, на практике различают 5 разных радиусов. Для каждой задачи – свой радиус.

В данной рассматриваемой задаче эквивалентная сфера - это такая сфера, на которой решение первой геодезической ( ставим вторую точку на заданных азимуте и дистанции от первой ) и проведение линии между двумя точками даст минимальное отклонение от решения этих же задач на эллипсоиде.

Тоже интересует данный вопрос. Как можно с вами связаться? Оставьте пожалуйста какой-нибудь контакт! Не могу послать личное сообщение..

Странно, очень странно ... Сколько я сам писал в личку другим, сколько мне писали - никогда никаких проблем не было ... Скорее всего, вы просто не ту кнопку там нажали. Просто поковыряйтесь в личке еще. В конце концов, вы же можете изложить свою проблему и здесь. Всегда рад помочь.
Ну, если совсем уж никак с личкой, посигнальте здесь, дам свой почтовый ящик.

Ну, теперь в заключение, - по обсуждению проблемы. Если очень кратко - проблему решил. Средние цифры такие - сравнение решений указанных задач дает их отклонения друг от друга, в среднем, - 1.0 ... 1.5 ... 2.0
километра на каждую тысячу километров дистанции. То есть, на 10 тысячах километрах - это будет до 20 километров отклонения. Величина-то, конечно, лишенная математической строгости, и, надо сказать, весьма условная, но меня вполне устраивает...
Фу - у - у, кажись, кончились мои мучения ...

[nickleb]

В R есть библиотека geosphere сферической тригонометрии:
https://cran.r-project.org/web/packages ... index.html
https://cran.r-project.org/web/packages ... sphere.pdf

[ericsson]

Если вас такой вопрос точности интересовал, то про это на Хабре статья была https://habrahabr.ru/post/239251/

[KotAlex]

В R есть библиотека geosphere сферической тригонометрии:
https://cran.r-project.org/web/packages ... index.html
https://cran.r-project.org/web/packages ... sphere.pdf

Спасибо за библиотеку. Вообще-то, у меня расчеты сферических треугольников реализованы в моей собственной библиотеке, но, как я бегло посмотрел описание по вашей ссылке, geosphere вдобавок к этому решает еще дополнительные задачи, которые могут представить интерес. Так что буду иметь в виду на будущее.