Как подключить карты в эллипсоидной проекции, если это не предусмотрено?
Или как подогнать тайлы Яндекс карт под проекцию OpenStreetMaps?
Вступление
Каждый раз, когда открываете какую-нибудь онлайн-карту, вы не скачиваете ее целиком. Для ускорения загрузки карты разделена на небольшие кусочки (тайлы), чтобы можно было скачать только нужную область. Проблем в том, что разрезать на эти квадратики можно несколькими способами. 
Большинство онлайн-карт «считают», что Земля имеет форму шара. Среди них, например, карты Google и OpenStreetMaps. Некоторые же, более дотошные, учитывают тот факт, что планета не является правильным шаром: как минимум, она сплющена у полюсов. Такая эллипсоидная проекция применяется, например, у карт Яндекса.
В результате клеточка с одним и тем же номером в разных проекциях будет показывать совершенно разные места. К примеру, вот тайл с номером 10427 по оси Х, 5119 по оси Y. Уровень масштаба 14. Слева — OSM, справа Яндекс.
И хотя большинство картографических движков умеют автоматически подгонять тайлы к нужной проекции, иногда может потребоваться сделать это вручную. Но как? Наиболее простой способ — просто сдвинуть тайлы на некоторое количество пикселей. В результате мы увидим на карте нужную местность. Конечно, если всматриваться, то можно разглядеть некоторые искажения. Но думаю, все таки, я думая, что бытовых задач, подобной точности будет более чем достаточно. Так что пора заканчивать со вступлением и начинать делать конвертер.
Методика
Для работы нам потребуется формула преобразования. Насколько я понял, ее вытянули прямо из кода страницы Яндекс карт, в те далекие времена, когда такое еще было вполне реально сделать. Ссылку на первоисточник я сейчас не найду, но на хабре эту формулу уже публиковали. Я практически не трогал: просто переписал на Swift и дал однобуквенным переменным более «говорящие» имена. По крайней мере, тем из них, которые удалось опознать.
Чтож, задача следующая. Нужно сделать делать конвертер, который принимает на вход номер тайла в стандартной проекции, а на выходе — номер тайла в эллипсоидной проекции и количество пикселей, на которое его необходимо сместить.
Итак. Для примера возьмет тайл с номером 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}
Что ж, надеюсь, кому-нибудь пригодятся изложенные здесь сведения. Желаю удачи с вашими экспериментами.
