Регулярно мы сталкиваемся с необходимостью преобразования из одной системы
координат в другую. Например, из простой и понятной долготы-широты в проекцию
меркатора.
Довольно часто эту задачу пытаются решать путем поиска формул пересчета, которые
в некоторых случаях довольно простые, а в некоторых довольно запутанные. Однако
для решения этой задачи есть простой и универсальный способ - библиотека
[PROJ.4](http://trac.osgeo.org/proj/).
Библиотека поддерживает множество разных типов проекций, вместе с ней
устанавливается обширная библиотека предопределенных проекций, но если
предустановленных проекций недостаточно, то возможно определение и
пользовательских проекций.
В качестве каталога предопределенных проекций можно использовать соответсвующий
раздел [spatialreference.org](http://www.spatialreference.org/ref/epsg/).
Библиотека поставляется с несколькими утилитами коммандной строки (пакет proj-
bin): `cs2cs`, `geod`. В качестве примера можно посмотреть список доступных
типов проекций:
$ cs2cs -lp
aea : Albers Equal Area
aeqd : Azimuthal Equidistant
...
wink2 : Winkel II
wintri : Winkel Tripel
Список доступных эллипсоидов:
$ cs2cs -le
MERIT a=6378137.0 rf=298.257 MERIT 1983
GS85 a=6378136.0 rf=298.257 Soviet Geodetic System 85
...
WGS84 a=6378137.0 rf=298.257223563 WGS 84
sphere a=6370997.0 b=6370997.0 Normal Sphere (r=6370997)
Рассмотрим примеры ее использования на python при помощи библиотеки
[pyproj](http://code.google.com/p/pyproj/).
В дистрибутивах, основанных на Debian эта библиотека расположенна в пакете
python-pyproj, но при необходимости она может быть так же установлена
в виртуальное окружение virtualenv через pypi - пакет pyproj.
Перепроецирование между системами координат
-------------------------------------------
Выполним преобразование из долготы-широты в проекцию меркатора на сфере,
которая наиболее часто используется в веб-картографии:
>>> from pyproj import Proj, transform
>>> lonlat = Proj(init="epsg:4326")
>>> sphmerc = Proj(init="epsg:3857")
>>> ll = (30, 59)
>>> sm = transform(lonlat, sphmerc, *ll)
>>> sm
(3372980.571036189, 8191201.661249711)
>>> transform(sphmerc, lonlat, *sm)
(29.999999999999996, 58.99999999999999)
Прямое, а затем обратное преобразование привело к исходному результату в
пределах погрешности.
Однако даже в веб-картографии не всегда используется одна проекция, например
[Яндекс-карты](http://maps.yandex.ru/) используют другую проекцию - меркатор на эллипсоиде. Зачем Яндекс
принял такое решение - остается только догадываться, но это делает исользование
слоев Яндекс-карт совместно с другими слоями довольно затруднительным.
Посмотрим чем отличается координаты в меркаторе на сфере и на эллипсоиде для
нашей точки:
>>> merc = Proj(proj="merc", ellps="WGS84")
>>> em = transform(lonlat, merc, *ll)
>>> sm[0] - em[0], sm[1] - em[1]
(0.0, 36659.2319375)
Помимо преобразования между спроецированными системами координат с помощью
библиотеки PROJ.4 можно выполнять преобразования между различными
географическими системами координат.
В качестве примера попробуем пересчитать координаты из WGS84 в систему координат
Пулково 1942 года. С параметрами перехода от WGS84 к Pulkovo-1942 есть некоторая
неопределенность, воспользуемся встроенным в PROJ.4 определением:
>>> pulkovo = Proj(init="epsg:4284")
>>> pp = transform(lonlat, pulkovo, *ll)
>>> pp
(30.002189268226005, 59.0000523360966)
**Задача**: Пепероецировать нашу тестовую точку в 36-ю зону проекции UTM.
Вычисления на геоиде
--------------------
Помимо работы с проекциями библиотека pyproj предоставляет возможность выполнять
прямое и обратное геодезическое преобразование. Грубо говоря, прямое
преобразованиие дает ответ на вопрос: где я окажусь, если будут вдигаться с
азимутом A из точки P расстояние L. Обратное, соответственно вычисляет
расстояние и азимут между двумя точками.
Исользуем обратное преобразование, для того чтобы оценить разницу между WGS84 и
Pulkovo 1942:
>>> from pyproj import Geod
>>> wgs84 = Geod(ellps="WGS84")
>>> wgs84.inv(*(ll + pp))
(87.34629057165381, -92.65183285869115, 125.96375775772009)
Первый элемент в кортеже: прямой азимут, второй - обратный, третий - расстояние
в метрах. Это как раз то расхождение в ~120 метров, с которым сталкиваются когда
пытаются использовать данные в СК Pulkovo совместно с данными WGS84.
**Задача**: Вычислить расстояние между Москвой и Санкт-Петербургом на геоиде WGS84 и
приближенной к нему сфере.