Проекция Меркатора: формула и код

Урок реализует проекцию Меркатора на чистом Python: из географических координат получаем плоские метры, которые ложатся на квадратную карту.

Прямая проекция Меркатора задаётся формулами $x = R\,\lambda$ и $y = R\,\ln\!\left(\tan\!\left(\frac{\pi}{4} + \frac{\phi}{2}\right)\right)$, где $\phi, \lambda$ — широта и долгота в радианах.

В прошлом уроке мы поняли, почему Меркатор искажает площади. Теперь реализуем саму проекцию. Идея проста: долгота переходит в горизонтальную координату линейно, а широта — нелинейно, через логарифм тангенса, и именно этот логарифм растягивает высокие широты.

Разбор формулы

Горизонтальная координата: $x = R\,\lambda$. Долгота просто масштабируется радиусом — меридианы на карте Меркатора равноотстоящие вертикальные линии.

Вертикальная координата: $y = R\,\ln\!\left(\tan\!\left(\frac{\pi}{4} + \frac{\phi}{2}\right)\right)$. Когда $\phi \to 90°$, аргумент тангенса стремится к $90°$, тангенс — к бесконечности, логарифм — тоже. Поэтому $y$ неограниченно растёт у полюсов, и карту приходится обрезать примерно на $\pm 85°$ — отсюда «квадратная» карта веба.

Реализация

import math

def mercator(lat_deg, lon_deg):
    R = 6378137.0  # радиус WGS84 в метрах
    lat = math.radians(lat_deg)
    lon = math.radians(lon_deg)
    x = R * lon
    y = R * math.log(math.tan(math.pi / 4 + lat / 2))
    return x, y

# Москва
x, y = mercator(55.7558, 37.6173)
print(f"x = {x:.1f} м")
print(f"y = {y:.1f} м")

# Точка на экваторе на той же долготе
x0, y0 = mercator(0.0, 37.6173)
print(f"y на экваторе = {abs(y0):.1f} м")

Вывод:

x = 4187538.7 м
y = 7509955.1 м
y на экваторе = 0.0 м

На экваторе $y = 0$ — это «нулевая линия» карты, как и должно быть: $\ln(\tan 45°) = \ln 1 = 0$. Для Москвы $y$ превышает $x$ — широта 56° уже заметно «вытянута» логарифмом.

Как работает под капотом

Существует красивое тождество: $\ln\!\left(\tan\!\left(\frac{\pi}{4} + \frac{\phi}{2}\right)\right) = \operatorname{arctanh}(\sin\phi)$. Оба варианта дают одно и то же, и в библиотеках встречается любой. Обратная проекция (из метров обратно в широту) использует функцию Гудермана и нужна, когда по клику на карте надо понять, какая это широта-долгота. Web Mercator (EPSG:3857) — это ровно эта формула с фиксированным $R = 6378137$, плюс упрощение «считаем Землю шаром, а не эллипсоидом» ради скорости тайлов.

Частые ошибки

Забыть перевести градусы в радианы. math.tan и math.log работают с радианами; подача градусов даёт мусор.
Подать широту $\pm 90°$. Логарифм уйдёт в бесконечность или выдаст ошибку — полюса в Меркаторе недостижимы, обрезайте до $\pm 85°$.
Перепутать порядок аргументов. Здесь функция принимает (широта, долгота); в GeoJSON порядок обратный.

Зачем уметь проецировать руками

Может показаться, что писать проекцию самому незачем — есть pyproj. Но понимание формулы окупается в практических ситуациях. Когда точки на вашей веб-карте «уезжают» от настоящего положения, причина почти всегда в путанице систем координат, и, зная, что Web Mercator берёт градусы и выдаёт метры в диапазоне примерно $\pm 20$ миллионов, вы сразу отличите «координаты в градусах подали как метры» от «перепутан порядок широты и долготы». Когда нужно вычислить, какой пиксель тайла соответствует точке, вы применяете ровно эту формулу плюс перевод метров в номера тайлов. А когда отлаживаете чужой геокод, способность в уме прикинуть, что Москва в Web Mercator это примерно $(4{,}19$ млн, $7{,}51$ млн$)$ метров, экономит часы.

Полезно осознавать и предел применимости. Меркатор хорош для навигации и тайлов, но катастрофичен для любой задачи про площади: считать по нему «сколько леса вырубили» — значит завысить вырубку в высоких широтах в разы. Поэтому реальный конвейер почти всегда двухсистемный: данные хранят и отдают в градусах WGS84, на экран показывают в Web Mercator, а измеряют — перепроецировав в локальную метровую зону. Меркатор в этой схеме отвечает строго за «показать», и попытка измерять прямо на нём — классическая ошибка новичка, которую теперь вы распознаете.

Итог

Меркатор: $x = R\lambda$ линейно, $y = R\ln(\tan(\frac{\pi}{4}+\frac{\phi}{2}))$ нелинейно.
На экваторе $y = 0$; у полюсов $y \to \infty$, поэтому карту обрезают на $\pm 85°$.
Перед тригонометрией градусы переводят в радианы.
Web Mercator (EPSG:3857) — эта формула с шаровым радиусом WGS84.