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

Или как подогнать тайлы Яндекс карт под проекцию 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}

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