Как рисовать прямоугольные полигоны на карте без искажений в проекции web mercator

Отрисовка прямоугольных полигонов на интерактивной карте оказалась не такой тривиальной задачей, как кажется на первый взгляд. Мы строили сетку из ячеек и для каждой центральной точки должны были получить прямоугольник фиксированного положения внутри глобальной сетки. Проблема проявилась сразу: чем севернее находился центр, тем визуально меньше становился прямоугольник; на юге же клетки начинали заметно перекрываться. Источник — масштабные искажения проекции Псевдо‑Меркатора (Web Mercator, EPSG:3857), где масштаб по широте растет как 1/cos φ, а картографический y нелинеен по широте.

Чтобы получить предсказуемые полигоны, мы полностью перешли к логике «сначала сетка в проекции, затем обратные преобразования в географические координаты». Земля разбивается на равномерную сетку по осям X и Y в пространстве Web Mercator, а затем каждая грань ячейки переводится обратно в долготу/широту. Ключевая идея: индексы row/col должны соответствовать равномерным шагам по осям проекции (в метрах или в нормированных координатах), а не в географических градусах, поскольку именно в меркаторском пространстве шаги равны визуально.

Определения:
- row и col — индексы ряда и колонки на глобальной сетке, покрывающей весь диапазон Web Mercator.
- resolution — степень разбиения по каждой оси; по сути, число логических делений 2^resolution на ось, что дает (2^resolution × 2^resolution) ячеек. Чем выше resolution, тем более детальна сетка.
- Проекция — WGS 84 / Pseudo‑Mercator, сферический вариант, достаточно точный и быстрый для отрисовки.

Почему возникало искажение. Если формировать прямоугольники «по градусам» вокруг центральной точки, одинаковые приращения по долготе дают разную физическую ширину в метрах на разных широтах (сжатие к полюсам). В Web Mercator долгота трансформируется линейно в x, а широта — нелинейно в y: y = R ln tan(π/4 + φ/2) для сферической модели. Значит, равные приращения φ не соответствуют равным шагам y. Решение — строить шаги в пространстве y, а не φ.

Мы используем две функции преобразования индексов в географические координаты. Сначала определяем нормированные доли по осям для текущего полигона, затем применяем обратные формулы проекции.

row_to_lat(row, resolution, shift):
- Вычисляем нормированную вертикальную долю v = (row + shift)/N, где N = 2^resolution — количество рядов.
- Преобразуем v в меркаторский y в нормализованной форме, затем применяем обратную функцию к широте: φ = atan(sinh(π · (1 − 2v))) в радианах, далее переводим в градусы. Это обеспечивает корректное распределение точек по широте с учетом нелинейности.

col_to_lon(col, resolution, shift):
- Горизонтальная доля u = (col + shift)/N.
- Долгота λ = 360° · u − 180°. Чтобы держать значения в диапазоне [−180°, 180°], при необходимости сдвигаем λ на ±360°.

Проверка пределов: верхняя граница самого верхнего ряда дает примерно +85.05113°, нижняя граница самого нижнего — около −85.05113°. Это естественное ограничение Web Mercator, который не охватывает полюса.

Далее функция сборки четырех углов прямоугольника — get_polygon(row, col, is_space_between):
- Берем верхнюю грань как row_to_lat(row, resolution, shift=1), нижнюю — row_to_lat(row, resolution, shift=0).
- Аналогично по долготе: левый край — col_to_lon(col, resolution, shift=0), правый — col_to_lon(col, resolution, shift=1).
- Если включен is_space_between, делаем небольшой отступ внутрь полигона. Проще всего уменьшить диапазоны по долготе и широте на фиксированную долю шага сетки (например, по 5–10% с каждой стороны), вычислив промежуточные lat/lon через дополнительные shift, расположенные ближе к центру.

Таким образом мы получаем четыре вершины: (lon_left, lat_top), (lon_right, lat_top), (lon_right, lat_bottom), (lon_left, lat_bottom). Эти значения применимы напрямую к полилиниям/полигонам на карте. Поскольку границы получены через обратную проекцию, визуальный размер ячейки на экране остаётся согласованным вдоль долгот, а вертикальные интервалы учитывают нелинейность по широте.

Альтернативный рабочий путь — перейти в метры проекции EPSG:3857 и оперировать реальными размерами:
- Конвертируете центральную точку (lat, lon) в (x, y) в метрах.
- Строите прямоугольник фиксированного размера в метрах вокруг (x, y).
- Конвертируете четыре угла назад в (lat, lon).
Этот вариант полезен, когда нужно задавать размер полигонов в физических единицах (например, 1×1 км) и гарантировать постоянство площади/ширины на карте.

Нюансы и грабли:
- Пересечение линии смены даты. Для ячеек вблизи ±180° может визуально возникать разрыв. Лучше хранить долготные границы как два сегмента, если полигоны «перескакивают» через −180°/180°. Либо нормализовать координаты во время рендеринга тайла.
- Стыковка соседних ячеек. Если отступ is_space_between выключен, плавающая точность может приводить к микрощелям или перекрытиям на экстремальном масштабе. Решение — единый источник вычисления углов для соседних полигонов, либо аккуратное округление координат.
- Производительность. Частые вызовы тригонометрии (atan, sinh) для миллионов ячеек дорогие. Эффективна векторизация (NumPy), кэширование строк lat по row и столбцов lon по col, предвычисление границ для каждого row/col один раз на слой.
- Ограничения широт. Не забывайте, что все значения должны строго лежать в пределах ±85.05113°. Любую «вылезающую» широту нужно обрезать до этих границ.
- Высокая resolution. При 2^15 и выше координатные шаги становятся очень маленькими. Используйте двойную точность и избегайте избыточных преобразований туда‑сюда, чтобы не накопить погрешность.

Пример логики без привязки к конкретной библиотеке:
1) N = 2^resolution.
2) lat_top = rad2deg(atan(sinh(π · (1 − 2 · (row + 1)/N))))
3) lat_bottom = rad2deg(atan(sinh(π · (1 − 2 · (row + 0)/N))))
4) lon_left = 360 · (col + 0)/N − 180
5) lon_right = 360 · (col + 1)/N − 180
6) Если нужен зазор g в долях ячейки, заменяем 0 на g и 1 на 1 − g в формулах 2–5.

Почему это устраняет исходное искажение. Мы больше не пытаемся рисовать «равноградусные» прямоугольники — вместо этого мы используем равные интервалы вдоль осей собственно проекции. Визуальный шаг по долготе и «меркаторской высоте» становится постоянным для заданного уровня resolution. Тем самым полигоны перестают сжиматься к северу и перестают расползаться к югу, а их взаимные границы совпадают по всей карте.

Если же цель — одинаковая физическая ширина/высота в метрах для всех ячеек независимо от широты, используйте метры проекции. Для заданного размера dx × dy в метрах:
- y_top = y_center + dy/2, y_bottom = y_center − dy/2
- x_left = x_center − dx/2, x_right = x_center + dx/2
- Далее обратная конверсия в (lat, lon).
Такой подход инвариантен к широте в смысле реальных размеров, что важно для задач плотности, расчёта площадей и корректной заливки.

Практические советы для интеграции:
- Рендер по тайлам. Предвычисляйте полигоны на уровне тайла/зум‑уровня и храните в кэше. При смене масштаба перестраивайте только требуемые ряды/колонки.
- Порядок вершин. Всегда используйте согласованный порядок (например, по часовой стрелке), чтобы не получить «вывернутые» полигоны в графическом движке и чтобы корректно работал заливочный алгоритм.
- Обход антимередиана. Для визуализации многоугольников, пересекающих 180°, разбивайте геометрию на две с корректной нормализацией долгот. Это заметно снижает артефакты в окрестностях TMS/XYZ тайлов.
- Параметр зазора. Даже минимальный зазор 1–2% убирает мерцание общих границ при сглаживании и устраняет «лесенки» на плотных масштабах.
- Точность цвета и статистики. Если ячейки кодируют статистику, храните ключевые метаданные вместе с индексами row/col. Это позволит быстро обновлять заливку без пересчётов геометрии.

Коротко о математике Web Mercator (сфера):
- x = R · λ, где λ в радианах (линейно по долготе).
- y = R · ln tan(π/4 + φ/2), φ — широта в радианах.
- Обратные: λ = x/R, φ = 2 · arctan(exp(y/R)) − π/2, или эквивалентно φ = atan(sinh(y/R)).
- Ограничение φ ∈ [−85.05113°, +85.05113°], поскольку tan(π/4 + φ/2) стремится к бесконечности у полюсов.

Итог. Чтобы рисовать прямоугольные полигоны без искажений, соответствующих ожиданиям интерфейса, опирайтесь на равномерную сетку в пространстве проекции и переводите её границы в широту/долготу через обратные формулы. Для задач с физическим размером используйте метры проекции, для задач «логической» сетки — доли индексов с нормализацией. Дополнительные улучшения — кэширование, контроль границ и обработка даталинии — позволяют добиться стабильной и гладкой отрисовки на всех масштабах.

1
1
Прокрутить вверх