Как подключить карты в эллипсоидной проекции, если это не предусмотрено?

Или как подогнать тайлы Яндекс карт под проекцию OpenStreetMaps?

Вступление

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



Большинство онлайн-карт “считают”, что Земля имеет форму шара. Среди них, например, карты Google и OpenStreetMaps. Некоторые же, более дотошные, учитывают тот факт, что планета не является правильным шаром: как минимум, она сплющена у полюсов. Такая эллипсоидная проекция применяется, например, у карт Яндекса.



В результате клеточка с одним и тем же номером в разных проекциях будет показывать совершенно разные места. К примеру, вот тайл с номером 10427 по оси Х, 5119 по оси Y. Уровень масштаба 14. Слева – OSM, справа Яндекс. Вместо города — какой-то лес.



И хотя большинство картографических движков умеют автоматически подгонять тайлы к нужной проекции, иногда может потребоваться сделать это вручную. Но как? Наиболее простой способ – просто сдвинуть тайлы на некоторое количество пикселей. В результате мы увидим на карте нужную местность. Конечно, если всматриваться, то можно разглядеть некоторые искажения. Но думаю, все таки, что бытовых задач, подобной точности будет более чем достаточно. Так что пора заканчивать со вступлением и начинать делать конвертер.

Методика

Для работы нам потребуется формула преобразования. Насколько я понял, ее вытянули прямо из кода страницы Яндекс карт, в те далекие времена, когда такое еще было вполне реально сделать. Ссылку на первоисточник я сейчас не найду, но на хабре эту формулу уже публиковали. Я ее практически не трогал: просто переписал на Swift и дал непонятным однобуквенным переменным более «говорящие» имена. По крайней мере, тем из них, которые удалось опознать. (Спасибо Erelen'у за помощь)

Чтож, задача следующая. Нужно сделать конвертер, который принимает на вход номер тайла в стандартной проекции, а на выходе – номер тайла в эллипсоидной проекции и количество пикселей, на которое его необходимо сместить.

Итак. Для примера возьмет тайл с номером X 10427, Y 5119, Z 14.

Действовать будем в два шага. Для начала, нужно найти координаты (широту и долготу) этого тайла. Например, координаты его левого-верхнего угла.

func tileNumberToCoordinates(tileX: Int, tileY: Int, mapZoom: Int) -> (lat_deg: Double, lon_deg: Double) {
        
        let n : Double = pow(2.0, Double(mapZoom))
        let lon = (Double(tileX) / n) * 360.0 - 180.0
        let lat = atan( sinh (.pi - (Double(tileY) / n) * 2 * Double.pi)) * (180.0 / .pi)
        
        return (lat, lon)
    }

Получаем на выходе (55.7889 49.1088). Теперь подставим полученные значения в нашу формулу. Уровень зума все тот же: 14-й.

func getWGS84Position(latitude: Double, longitude: Double, zoom: Int) -> (x:Int, y:Int, offsetX:Int, offsetY:Int) {
        
        // Earth vertical and horisontal radiuses
        let radiusA = 6378137.0
        let radiusB = 6356752.0
        
        let latitudeInRadians = latitude * Double.pi / 180
        
        let yCompressionOfEllipsoid = sqrt( pow(radiusA, 2.0) - pow(radiusB, 2.0)) / radiusA
        
        // I really don't know what the name of this variable mean =(
        let m2 = log((1 + sin(latitudeInRadians)) / (1 - sin(latitudeInRadians))) / 2 - yCompressionOfEllipsoid * log((1 + yCompressionOfEllipsoid * sin(latitudeInRadians)) / (1 - yCompressionOfEllipsoid * sin(latitudeInRadians))) / 2
        
        // x count = y count
        let xTilesCountForThisZoom = Double(1 << zoom)
        
        //Tile numbers in WGS-84 proection
        let xTileNumber = floor((longitude + 180) / 360 * xTilesCountForThisZoom)
        let yTileNumber = floor(xTilesCountForThisZoom / 2 - m2 * xTilesCountForThisZoom / 2 / Double.pi)
        
        //Offset in pixels of the coordinate of the
        //left-top corner of the OSM tile
        //from the left-top corner of the WGS-84 tile
        let offsetX = floor(((longitude + 180) / 360 * xTilesCountForThisZoom - xTileNumber) * 256)
        let offsetY = floor(((xTilesCountForThisZoom / 2 - m2 * xTilesCountForThisZoom / 2 / Double.pi) - yTileNumber) * 256)
        
        return (Int(xTileNumber), Int(yTileNumber), Int(offsetX), Int(offsetY))
    }

Получаем (10427, 5133, 0, 117). Это значит, что нам нужен Яндекс тайл с номером X 10427, Y 5133, Z 14. И если его сместить его на 0 пикселей влево и на 117 пикселей вверх, то он займет нужное место.

И что с этим делать?

Если вы пишете свой навигатор и у вас есть возможность повлиять на отображение карты, то все сравнительно просто. Рассчитайте смещение для любого из тайлов на экране. И сместите элемент с картой на это количество пикселей любым удобным способом.

А вот если доступа к коду у вас нет, то придется вводить промежуточное звено между сервером карты и навигатором. Например, я сделал для этого простенький сервер. Он принимает на вход номер искомого тайла, вычисляет номер тайла в эллипсоидной проекции. Скачивает его и три соседние тайла. Склеивает эти четыре тайла в один большой, а затем вырезает из него нужный фрагмент и возвращает навигатору пользователя.



Результат и оригинал:



Разумеется, все эти операции требуют дополнительных затрат по времени. По этим ссылкам можно оценить, с какой скоростью сервер “фотошопит” карту в реальном времени:

https://anygis.ru/api/v1/Yandex_map/{x}/{y}/{z}

https://anygis.ru/api/v1/Yandex_sat_clean/{x}/{y}/{z}

Что ж, надеюсь, кому-нибудь пригодятся изложенные здесь сведения. Желаю удачи с вашими экспериментами.