Формулы проекции Гаусса-Крюгера

[Vladislav Kalinnikov]

Посоветуйте хорошие формулы для проекции Гаусса-Крюгера (для прямого и обратного преобразования), которые позволяли бы вести пересчет координат не менее чем в 10 градусных зонах.

[stout]

Почти все современные алгоритмы для широкой зоны так или иначе основаны на оригинальной работе Крюгера.
Он же описан в книге Морозова Курс сфероидической геодезии.

Одна из лучших библиотек:
http://geographiclib.sourceforge.net/
Огромный список литературы, тестовые примеры, аналитически точный алгоритм проекции ТМ, хорошая реализация.

В PROJ.4 тоже добавили алгоритм Крюгера.
http://trac.osgeo.org/proj/browser/trun ... j_etmerc.c
Причём, в отличии от довольно нового стандарта NGA
( http://earth-info.nga.mil/GandG/publica ... UTMUPS.pdf ),
суммирование выполняется по алгоритму Кленшоу.

Одними из первых начали применять этот алгоритм шведы
http://www.lantmateriet.se/globalassets ... ection.pdf

Интересный алгоритм приведён в
C++ and Java Code for Geodesic and Meridian Arc Computations -- by Klaus Hehl
http://www.ngs.noaa.gov/gps-toolbox/Hehl1.htm
Саму статью можно скачать с sci-hub.io

Наивная реализация – EPSG Wide-Zone TM Investigation and Recommendation
http://www.hydrometronics.com/downloads ... M%20r2.pdf

Но лично мне нравится алгоритм в статье: Wide Zone Transverse Mercator Projection
http://www.dfo-mpo.gc.ca/Library/337182.pdf
Прост в реализации, достаточно быстр, в отличии от формул Крюгера нет надобности в вычислении коэффициентов как рядов по третьему сжатию.
Если посмотреть на PROJ.4, то где-то там есть сообщение об исправлении ошибки в вычислении одного из коэффициентов. (через несколько лет после первого релиза)

вот только в подпрограмме SUBROUTINE eGud( E,Psi, Phi, Dphi)
начальное значение для SinPhi
SinPhi = CDTANH( Psi)

считал бы как
sinφ = t + e^2×t×(1-t^2)×[1 + G1×e^2 + G2×e^4]
G1 = 1 - (5/3)×t^2
G2 = 1 - 4×t^2 + (16/5)×t^4

где t = tanh(ψ)

З.Ы. В статье по странной привычке геодезистов интегралом Валлиса называется то, что им не является.

[Vladislav Kalinnikov]

Спасибо Александр Юрьевич за ответ настоящего специалиста!

[stout]

Не за что.
Есть у меня ещё один любопытный алгоритм, как мне кажется.
Вот только статью уже лет пять всё не соберусь написать.
Если коротко, то формулы очень простые.

MyGK.png

Уже порядок n=2 обеспечивает более чем достаточную для практики точность.
На подмножестве тестовых данных для GeographicLib (при ограничении удаления от осевого меридиана в 75°)
получаются такие результаты

Northing_order_2_Lon_75.png

Easting__order_2_Lon_75.png

А при n=3

Northing_order_3_Lon_75.png

Easting__order_3_Lon_75.png

Причём, все "выбросы" есть следствие близости точки проектирования к полюсу.
Расстояния от полюса – сотые и тысячные доли угловой секунды.
То есть то, чего на практике быть не может.

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

Order = 3

Cphi =  1.0050525017882656148
C[0] = -7.0235486287468349323
C[0] = -7.023548628746834937     D[0] =  595.51444712599999998
C[1] = 0.066302967647480009766   D[1] = 3171.9701604184324712
C[2] = 0.45018388791371950894    D[2] =  975.51373108852711702
C[3] = 5.5020092713973698037     D[3] =  626.6421830947292918
VERSION : 5.1.0
Для продолжения нажмите любую клавишу . . .
n = 209950    Lat = 89.999999989961     Lon = 17.453254127432
Xn = 9997964.941951752   Ye = 0.0003361730433   x-Xn = 0.00106925   y-Ye = -0.000336173
n = 209986    Lat = 89.999999989211     Lon = 65.343754462988
Xn = 9997964.942518477   Ye = 0.0010947587778   x-Xn = 0.00050252   y-Ye = -0.00109476
n = 210001    Lat = 89.999999982901     Lon = 8.73145127768
Xn = 9997964.941134031   Ye = 0.0002898065862   x-Xn = 0.00107105   y-Ye = 0.000526105
n = 210072    Lat = 89.999999978859     Lon = 67.301195263
Xn = 9997964.942110159   Ye = 0.0021775575134   x-Xn = -0.000832504   y-Ye = -0.00127155
n = 210119    Lat = 89.999999989722     Lon = 23.814316771947
Xn = 9997964.941971168   Ye = 0.0004633432899   x-Xn = 0.00104983   y-Ye = -0.000463343
n = 210133    Lat = 89.99999998947     Lon = 21.13572398521
Xn = 9997964.941924419   Ye = 0.0004239202218   x-Xn = 0.00109658   y-Ye = -0.00042392
n = 210142    Lat = 89.999999989913     Lon = 3.269758812007
Xn = 9997964.941896625   Ye = 6.42355374e-005   x-Xn = 0.00112437   y-Ye = -6.42355e-005
n = 210165    Lat = 89.999999987618     Lon = 60.910821572258
Xn = 9997964.942348896   Ye = 0.0012080660827   x-Xn = 0.000672102   y-Ye = -0.00120807
n = 210322    Lat = 89.999999987322     Lon = 57.551979936139
Xn = 9997964.942261539   Ye = 0.0011945015199   x-Xn = 0.000759459   y-Ye = -0.0011945
n = 210379    Lat = 89.999999988312     Lon = 2.8546406241
Xn = 9997964.94171766   Ye = 6.49898199e-005   x-Xn = 0.00130334   y-Ye = -6.49898e-005
n = 210391    Lat = 89.999999981135     Lon = 8.148489986313
Xn = 9997964.940935998   Ye = 0.0002985398662   x-Xn = 0.001296   y-Ye = 0.000490457
n = 210478    Lat = 89.999999980993     Lon = 9.14770220425
Xn = 9997964.940925869   Ye = 0.0003373745642   x-Xn = 0.00126063   y-Ye = 0.000497124
n = 210506    Lat = 89.999999983887     Lon = 65.838301416234
Xn = 9997964.942284641   Ye = 0.001641401445   x-Xn = -0.00102521   y-Ye = -0.000700799
n = 210531    Lat = 89.999999989125     Lon = 9.646819943799
Xn = 9997964.941823981   Ye = 0.0002034665287   x-Xn = 0.00119702   y-Ye = -0.000203467
n = 210534    Lat = 89.999999986907     Lon = 50.351290215997
Xn = 9997964.942088239   Ye = 0.0011255624153   x-Xn = 0.000932759   y-Ye = -0.00112556
n = 210616    Lat = 89.999999980596     Lon = 66.643720386703
Xn = 9997964.942162117   Ye = 0.001988919096   x-Xn = -0.000892764   y-Ye = -0.00106703
n = 210649    Lat = 89.999999987874     Lon = 7.699729002035
Xn = 9997964.941679345   Ye = 0.0001813921156   x-Xn = 0.00134165   y-Ye = -0.000181392
n = 210715    Lat = 89.999999982664     Lon = 72.135639013664
Xn = 9997964.942427239   Ye = 0.0018422304213   x-Xn = -0.00108392   y-Ye = -0.00107008

Здесь n – порядковый номер точки в тестовом файле; Lat,Lon – широта и долгота в градусах, соответственно;
Xn,Ye – тестовые значения, а (x-Xn) и (y-Ye) соответствующие ошибки.