Wie verbinde ich Karten in einer Ellipsoidprojektion, wenn dies nicht vorgesehen ist?

Oder wie passt man Yandex-Kartenkacheln an die OpenStreetMaps-Projektion an?

Eintrag

Jedes Mal, wenn Sie eine Art Online-Karte öffnen, wird diese nicht vollständig heruntergeladen. Um das Laden der Karte zu beschleunigen, ist diese in kleine Teile (Kacheln) unterteilt, so dass nur der gewünschte Bereich heruntergeladen werden kann. Das Problem ist, dass Sie auf verschiedene Arten in diese Quadrate schneiden können.



Die meisten Online-Karten „denken“, dass die Erde die Form einer Kugel hat. Darunter zum Beispiel Google Maps und OpenStreetMaps. Einige, akribischer, berücksichtigen die Tatsache, dass der Planet nicht der richtige Ball ist: Zumindest ist er an den Polen abgeflacht. Eine solche Ellipsoidprojektion wird beispielsweise in Yandex-Karten verwendet.



Infolgedessen zeigt eine Zelle mit derselben Nummer in verschiedenen Projektionen völlig unterschiedliche Stellen. Beispiel: Hier ist eine Kachel mit der Nummer 10427 auf der X-Achse und 5119 auf der Y-Achse. Maßstabsebene 14. Links - OSM, rechts Yandex. Anstelle einer Stadt - eine Art Wald.



Obwohl die meisten Kartenmodule Kacheln automatisch an die gewünschte Projektion anpassen können, müssen Sie dies manchmal manuell tun. Aber wie? Am einfachsten ist es, die Kacheln einfach um eine bestimmte Anzahl von Pixeln zu verschieben. Als Ergebnis sehen wir den gewünschten Bereich auf der Karte. Wenn Sie genau hinschauen, können Sie natürlich einige Verzerrungen feststellen. Trotzdem denke ich, dass alltägliche Aufgaben von ähnlicher Genauigkeit mehr als genug sein werden. Also ist es Zeit, mit der Einführung fertig zu werden und mit der Herstellung des Konverters zu beginnen.

Methodik

Für die Arbeit benötigen wir eine Umrechnungsformel. Soweit ich weiß, wurde es in jenen Tagen, als es noch möglich war, direkt aus dem Seitencode von Yandex Maps gezogen. Ich werde den Link zur Quelle momentan nicht finden, aber diese Formel wurde bereits auf dem Hub veröffentlicht. Ich habe es praktisch nicht angerührt: Ich habe es einfach auf Swift umgeschrieben und unverständlichen Variablen mit einem Buchstaben mehr "sprechende" Namen gegeben. Zumindest diejenigen, die sich identifizieren konnten. (Danke an Erelen für die Hilfe)

Nun, die Aufgabe ist wie folgt. Wir müssen einen Konverter erstellen, der die Kachelnummer in der Standardprojektion als Eingabe und die Kachelnummer in der Ellipsoidprojektion sowie die Anzahl der Pixel, um die sie in die Ausgabe verschoben werden muss, verwendet.

Also Nehmen Sie zum Beispiel ein Plättchen mit der Nummer X 10427, Y 5119, Z 14.

Wir werden in zwei Schritten handeln. Zunächst müssen Sie die Koordinaten (Breiten- und Längengrade) dieser Kachel ermitteln. Zum Beispiel die Koordinaten der linken oberen Ecke.

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) } 

Wir erhalten die Ausgabe (55.7889 49.1088). Jetzt ersetzen wir die erhaltenen Werte in unserer Formel. Die Zoomstufe ist immer noch dieselbe: 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)) } 

Wir erhalten (10427, 5133, 0, 117). Dies bedeutet, dass wir ein Yandex- Plättchen mit der Nummer X 10427, Y 5133, Z 14 benötigen. Wenn Sie es um 0 Pixel nach links und 117 Pixel nach oben verschieben, wird es die richtige Stelle einnehmen.

Und was tun?

Wenn Sie Ihren Navigator schreiben und die Anzeige der Karte beeinflussen können, ist alles relativ einfach. Berechnen Sie den Versatz für alle Kacheln auf dem Bildschirm. Verschieben Sie das Element mit der Karte auf eine beliebige Weise um diese Pixelanzahl.

Wenn Sie jedoch keinen Zugriff auf den Code haben, müssen Sie eine Zwischenverbindung zwischen dem Kartenserver und dem Navigator eingeben. Zum Beispiel habe ich einen einfachen Server dafür gemacht. Es übernimmt die Eingabe der Nummer der gewünschten Kachel, berechnet die Nummer der Kachel in der Ellipsoidprojektion. Außerdem werden drei benachbarte Kacheln heruntergeladen. Klebt diese vier Kacheln zu einer großen, schneidet dann das gewünschte Fragment aus und gibt es an den Benutzer-Navigator zurück.



Ergebnis und Original:



Natürlich erfordern alle diese Operationen zusätzliche Zeitkosten. Mit diesen Links können Sie die Geschwindigkeit, mit der der Server die Karte "Photoshop" in Echtzeit:

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

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

Nun, ich hoffe, dass die hier präsentierten Informationen für jemanden nützlich sind. Viel Glück bei deinen Experimenten.