Проекция географической точки на изображение карты

Задача о выгрузке информации о странах получила свое продолжение. Собранные данные о геолокации стран (координаты столицы, крупных городов) нужно применить, а именно, сделать соответствующие отметки на карте.

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

📚 С решением задачи нам поможет картографическая проекция Меркатора, а точнее ее вариант WGS 84 / Pseudo-Mercator, который используется в современных картографических приложениях и сервисах. WGS - всемирная геодезическая система, в настоящее время являющаяся стандартном в картографии, геодезии и спутниковой навигации. Согласно этой системе, Земля - эллипсоид, а экватор и нулевой меридиан - опорные точки для широты и долготы соответственно.

❗️❗️ Обязательным условием является работа с изображением карты, которое соответствует WGS 84.

Подходящее изображение карты можно взять, например, здесь.

🚀 Перейдем к решению.

Поскольку точка отсчета - это пересечение экватора и нулевого меридиана, а нам удобнее было бы работать с левым верхним углом изображения, как началом координат, то для смещения по оси X используем ложный сдвиг в восточном направлении, а по оси Y ложный сдвиг в северном направлении. Ложный сдвиг - это линейное значение, применяемое к началу координат, которое позволяет гарантировать, что значения по конкретной оси будут положительными, что нам и нужно при работе с изображением.

Сама формула расчета давно известна, я просто приведу вариант ее реализации на языке python применительно к задаче, а также чтение и отрисовку тестовых данных с использованием библиотеки opencv.

✅ Тестовые данные

data = [
    {"lat": 55.75, "lon": 37.6, "label": "Moscow"},
    {"lat": 55.0415, "lon": 82.9346, "label": "Novosibirsk"},
    {"lat": 56.8519, "lon": 60.6122, "label": "Yekaterinburg"},
    {"lat": 51.3026, "lon": 0.0739, "label": "London"},
    {"lat": 48.48, "lon": 2.20, "label": "Paris"},
    {"lat": 40.24, "lon": -3.41, "label": "Madrid"},
]

✅ Функция построения проекции

# lat - широта
# lon - долгота
# width - ширина изображения карты
# height - высота изображения карты
def latlon_to_offset(lat, lon, width, height):
    # ложный сдвиг в восточном направлении
    FALSE_EASTING = 180

    # долгота меридиана в точке отсчета
    # 0, т.к. у нас точкой отсчета является левый край карты
    lon0 = 0

    # радиус Земли в проекции
    # ширина изображения здесь является длиной окружности
    radius = width / (2 * math.pi)

    # переводим широту и долготу в радианы
    radian_lat = math.radians(lat)
    radian_lon = math.radians(lon - lon0 + FALSE_EASTING)

    # определяем положение по X
    offset_x = radius * radian_lon

    # определяем положение по Y относительно экватора
    offset_y_from_equator = radius * math.log(math.tan(math.pi / 4 + radian_lat / 2))

    # определяем положение по Y относительно верхней границы изображения
    offset_y = height / 2 - offset_y_from_equator

    return int(offset_x), int(offset_y)

✅ Читаем изображение, применяем функцию на тестовых данных, отрисовываем результат и сохраняем

src = cv2.imread("1084px-Web_maps_Mercator_projection_SW.jpg", cv2.IMREAD_COLOR)
h, w, _ = src.shape

for item in data:
    x, y = latlon_to_offset(item["lat"], item["lon"], w, h)
    src = cv2.circle(src, (x, y), 5, (0, 0, 255), 2)

cv2.imwrite("result.jpg", src)

✅ Проверяем результат