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

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

Вступление


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

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

6kdcageqgsw_y7j6xbhpehkku-y.jpeg

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

zgsmfrsatg1gwtyofaanw-ymqfk.jpeg

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

7c-jcrrc1osrir4ttvlk-q2gsq8.jpeg

Методика


Для работы нам потребуется формула преобразования. Насколько я понял, ее вытянули прямо из кода страницы Яндекс карт, в те далекие времена, когда такое еще было вполне реально сделать. Ссылку на первоисточник я сейчас не найду, но на хабре эту формулу уже публиковали. Я практически не трогал: просто переписал на 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 пикселей вверх, то он займет нужное место.

kzvpfty0jfvlc_igwvzugk2qktw.jpeg

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


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

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

ns5c6kacnerrjgsi_gyf2rxf99u.jpeg

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

jie0ankna-4goxvdxvxsbe-qzhc.jpeg

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

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

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

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

© Habrahabr.ru