Ich versuche, eine Reihe von Bildkacheln zu georeferenzieren, die von einer Google Maps-Website heruntergeladen wurden. Nur basierend auf jedem Bild „s Google Maps-Rasterreferenz und Zoomstufe, wie kann ich die realen Koordinaten berechnen?

Die Seite unter http://www.maptiler.org/google-maps-coordinates-tile-bounds-projection zeigt die Koordinaten und die Google Maps-Rasterreferenz für Web-Mercator-Kacheln an.

Die mittlere Kachel lautet beispielsweise x=49, y=30, z=6 und hat die Eckkoordinaten 10644926.307106785, 626172.1357121654, 11271098.44281895, 1252344.27142432:

Geben Sie hier die Bildbeschreibung ein

Können die Eckkoordinaten ausschließlich anhand der Gitterreferenz und der Zoomstufe (49, 30, 6) berechnet werden?

Dies ähnelt Wie kann ich eine Web-Mercator-Kachel mit gdal korrekt georeferenzieren? Ich benötige jedoch eine vollständig programmatische Lösung, ohne online nachschlagen zu müssen. (Wenn ich online nach Werten suchen muss, muss dies programmgesteuert sein.)

Bearbeiten: Dieser Screenshot zeigt die Position einer Beispielkachel (Zoomstufe 11, x = 1593, y = 933). . Wie würde ich diese Kachel georeferenzieren?

Geben Sie hier die Bildbeschreibung ein

Kommentare

  • Kommt es als echte Kachel? Keine Umgebung?
  • @FelixIP Ja, alles was ich weiß ist die Gitterreferenz und die Zoomstufe jedes einzelnen JPG. Ich ‚ habe einen Screenshot hinzugefügt, um dies zu demonstrieren. PS Mir ist klar, dass dies etwas zwielichtig ist, aber manchmal muss es getan werden – und es ist ‚ eine interessante Übung!
  • Ich denke, Sie sollten in der Lage sein, sie zu bekommen Die Eckkoordinaten verwenden meines Wissens nur die Zoomstufen und die Gitterreferenz. Das Gitterreferenzsystem bleibt also konstant, daher sollte es eine Möglichkeit geben, dies zu tun. Da Zoomstufe 0 die Welt in 1 Kachel und Zoom 1 in 4 usw. hat, sollte es möglich sein, die Koordinaten zu erhalten.
  • @StephenLead versucht möglicherweise, die Welt in jedem Zoom in 4 Teile zu unterteilen ( Iteration) und erhalten die Koordinaten damit irgendwie. Nach meinem Verständnis muss es ‚ eine iterative Funktion sein.
  • dürfen Sie dies tun? Das heißt, erlauben die Google Maps-Begriffe dies?

Antwort

Diese Seite http://www.maptiler.org/google-maps-coordinates-tile-bounds-projection bietet in Python geschriebene Algorithmen zum Berechnen von Kachelgrenzen in EPSG: 900913-Koordinaten und in Latutude / Longitude unter Verwendung des WGS84-Datums.

Die Klasse GlobalMercator in diesem Skript http://www.maptiler.org/google-maps-coordinates-tile-bounds-projection / globalmaptiles .py enthält zwei Methoden: TileBounds() und TileLatLonBounds(), um dies zu tun.

Eine modifizierte Version des Skripts zum Anzeigen der Koordinaten der Kachelgrenzen:

 #!/usr/bin/env python ############################################################################### # $Id$ # # Project: GDAL2Tiles, Google Summer of Code 2007 & 2008 # Global Map Tiles Classes # Purpose: Convert a raster into TMS tiles, create KML SuperOverlay EPSG:4326, # generate a simple HTML viewers based on Google Maps and OpenLayers # Author: Klokan Petr Pridal, klokan at klokan dot cz # Web: http://www.klokan.cz/projects/gdal2tiles/ # ############################################################################### # Copyright (c) 2008 Klokan Petr Pridal. All rights reserved. # # Permission is hereby granted, free of charge, to any person obtaining a # copy of this software and associated documentation files (the "Software"), # to deal in the Software without restriction, including without limitation # the rights to use, copy, modify, merge, publish, distribute, sublicense, # and/or sell copies of the Software, and to permit persons to whom the # Software is furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included # in all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS # OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL # THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING # FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER # DEALINGS IN THE SOFTWARE. ############################################################################### """ tilebounds.py Adapted from: http://www.maptiler.org/google-maps-coordinates-tile-bounds-projection/globalmaptiles.py Global Map Tiles as defined in Tile Map Service (TMS) Profiles ============================================================== Functions necessary for generation of global tiles used on the web. It contains classes implementing coordinate conversions for: - GlobalMercator (based on EPSG:900913 = EPSG:3785) for Google Maps, Yahoo Maps, Microsoft Maps compatible tiles - GlobalGeodetic (based on EPSG:4326) for OpenLayers Base Map and Google Earth compatible tiles More info at: http://wiki.osgeo.org/wiki/Tile_Map_Service_Specification http://wiki.osgeo.org/wiki/WMS_Tiling_Client_Recommendation http://msdn.microsoft.com/en-us/library/bb259689.aspx http://code.google.com/apis/maps/documentation/overlays.html#Google_Maps_Coordinates Created by Klokan Petr Pridal on 2008-07-03. Google Summer of Code 2008, project GDAL2Tiles for OSGEO. In case you use this class in your product, translate it to another language or find it usefull for your project please let me know. My email: klokan at klokan dot cz. I would like to know where it was used. Class is available under the open-source GDAL license (www.gdal.org). """ import math class GlobalMercatorLight(object): def __init__(self, tileSize=256): "Initialize the TMS Global Mercator pyramid" self.tileSize = tileSize self.initialResolution = 2 * math.pi * 6378137 / self.tileSize # 156543.03392804062 for tileSize 256 pixels self.originShift = 2 * math.pi * 6378137 / 2.0 # 20037508.342789244 def MetersToLatLon(self, mx, my ): "Converts XY point from Spherical Mercator EPSG:900913 to lat/lon in WGS84 Datum" lon = (mx / self.originShift) * 180.0 lat = (my / self.originShift) * 180.0 lat = 180 / math.pi * (2 * math.atan( math.exp( lat * math.pi / 180.0)) - math.pi / 2.0) return lat, lon def PixelsToMeters(self, px, py, zoom): "Converts pixel coordinates in given zoom level of pyramid to EPSG:900913" res = self.Resolution( zoom ) mx = px * res - self.originShift my = py * res - self.originShift return mx, my def TileBounds(self, tx, ty, zoom): "Returns bounds of the given tile in EPSG:900913 coordinates" minx, miny = self.PixelsToMeters( tx*self.tileSize, ty*self.tileSize, zoom ) maxx, maxy = self.PixelsToMeters( (tx+1)*self.tileSize, (ty+1)*self.tileSize, zoom ) return ( minx, miny, maxx, maxy ) def TileLatLonBounds(self, tx, ty, zoom ): "Returns bounds of the given tile in latutude/longitude using WGS84 datum" bounds = self.TileBounds( tx, ty, zoom) minLat, minLon = self.MetersToLatLon(bounds[0], bounds[1]) maxLat, maxLon = self.MetersToLatLon(bounds[2], bounds[3]) return ( minLat, minLon, maxLat, maxLon ) def Resolution(self, zoom ): "Resolution (meters/pixel) for given zoom level (measured at Equator)" # return (2 * math.pi * 6378137) / (self.tileSize * 2**zoom) return self.initialResolution / (2**zoom) #--------------------- if __name__ == "__main__": import sys, os def Usage(s = ""): print "Usage: tilebounds.py tx ty zoomlevel" print sys.exit(1) profile = "mercator" zoomlevel = None tx, ty = None, None argv = sys.argv i = 1 while i < len(argv): arg = argv[i] if tx is None: tx = float(argv[i]) elif ty is None: ty = float(argv[i]) elif zoomlevel is None: zoomlevel = int(argv[i]) else: Usage("ERROR: Too many parameters") i = i + 1 if profile != "mercator": Usage("ERROR: Sorry, given profile is not implemented yet.") if zoomlevel == None or tx == None or ty == None: Usage("ERROR: Specify at least "lat", "lon" and "zoomlevel".") tz = zoomlevel mercator = GlobalMercatorLight() minx, miny, maxx, maxy = mercator.TileBounds( tx, ty, tz ) print "Bounds of the given tile in EPSG:900913 coordinates: " print " upper-left corner: " print (minx, miny) print " bottom-right corner: " print (maxx, maxy) print minLat, minLon, maxLat, maxLon = mercator.TileLatLonBounds( tx, ty, tz ) print "Bounds of the given tile in latitude/longitude using WGS84 datum: " print " upper-left corner: " print (minLat, minLon) print " bottom-right corner: " print (maxLat, maxLon) print  

Verwendung: tilebounds xTile yTile zoom.

Zum Beispiel die Ausgabe für die Kachel x=49, y=30, z=6 lautet wie folgt:

$./tilebounds.py 49 30 6 Bounds of the given tile in EPSG:900913 coordinates: upper-left corner: (10644926.307106785, -1252344.271424327) bottom-right corner: (11271098.44281895, -626172.1357121654) Bounds of the given tile in latitude/longitude using WGS84 datum: upper-left corner: (-11.178401873711772, 95.625) bottom-right corner: (-5.61598581915534, 101.25000000000001) 

Eine Kachel kann mit der URL http://mt.google.com/vt?x=xTile&y=yTile&z=zoom heruntergeladen werden Es ist verboten, direkt darauf zuzugreifen.

Die Software wurde ursprünglich von Klokan Petr Přidal geschrieben.

Ich hoffe, das könnte helfen!

Kommentare

  • Willkommen bei GIS Stack Exchange, und vielen Dank dafür. Meine Anforderung ist erfüllt, da ich ‚ nicht mehr an diesem Projekt arbeite, aber diese Informationen sollten sich hoffentlich in Zukunft für andere als nützlich erweisen

Antwort

Wirf dies einfach als Teil der Diskussion über Kacheln, zwischengespeicherte Kacheln, Georeferenzierung usw. ein.

Kachelschemata von was Ich weiß, dass Sie kein XY-Georeferenzierungssystem per se verwenden …

Wenn Sie beispielsweise die ESRI-Grundkarte in QGIS laden, wird eine Kachel-Stache-Ordnerstruktur wie die folgende erstellt:

Geben Sie hier die Bildbeschreibung ein.

Wenn Sie die Bilder der obersten Ebene anzeigen, haben Sie hat die Kacheln angefordert und sie werden in diesen verschiedenen Ordnern gespeichert.

Hier ist der Inhalt des f / d-Ordners:

Geben Sie hier die Bildbeschreibung ein.

Beachten Sie, dass die Datei (Weltdatei usw.) keine Georeferenzierungsinformationen enthält.

Wann Sie bringen diese Datei in MS Beim Malen können Sie beispielsweise das Bild bearbeiten und erneut speichern:

Geben Sie hier die Bildbeschreibung ein

…Wenn Sie also die Karte in QGIS erneut anzeigen – das die Kacheln rendert -, sehen Sie Ihre „bearbeitete“ Kachel:

Geben Sie hier die Bildbeschreibung ein.

Ich denke also, ich sage nur, wenn Sie einen Satz Kacheln von Google Earth hätten, könnten Sie die Kachel theoretisch neu erstellen strukturieren und einfach die Kacheln in diese Ordner legen, und sie würden in QGIS zeichnen …

Es ist die Software selbst, die weiß, wo die Kacheln basierend auf dem Caching-Schema und der Ordnerstruktur gezeichnet werden sollen … Dies geht wieder auf meine Tage zurück, als ich ArcGIS Server-Kurse unterrichtete.

Sie können sie dann mit Georeferenzierungsinformationen aus QGIS exportieren und an anderer Stelle verwenden.

Andererseits können Sie auch einfach Google Maps in QGIS laden und das Bild als TIF mit einer Weltdatei exportieren …

(Wiederum keine wirkliche Antwort, nur einige Diskussionen …)

Kommentare

  • Entschuldigung für die späte Antwort. Um ehrlich zu sein, kann ich mich ‚ nicht einmal daran erinnern, warum ich das gefragt habe, da das Projekt schon lange vorbei ist;) Ich erinnere mich jedoch, dass ich jedes Voodoo neu erstellen musste, das die Software wann macht Es weiß, wie Bilder basierend auf der Ordnerstruktur platziert werden.

Antwort

Dieses Problem ist kürzlich aufgetreten Der Code in wiki.openstreetmap.org

Hier ist der Java-Code:

 class BoundingBox { double north; double south; double east; double west; } BoundingBox tile2boundingBox(final int x, final int y, final int zoom) { BoundingBox bb = new BoundingBox(); bb.north = tile2lat(y, zoom); bb.south = tile2lat(y + 1, zoom); bb.west = tile2lon(x, zoom); bb.east = tile2lon(x + 1, zoom); return bb; } static double tile2lon(int x, int z) { return x / Math.pow(2.0, z) * 360.0 - 180; } static double tile2lat(int y, int z) { double n = Math.PI - (2.0 * Math.PI * y) / Math.pow(2.0, z); return Math.toDegrees(Math.atan(Math.sinh(n))); } 

Antwort

Den zugrunde liegenden mathematischen Ansatz finden Sie in meiner Antwort auf eine verwandte Frage hier: Berechnen des Begrenzungsrahmens aus bekannten Mittelkoordinaten und Zoom

Sie müssen in der Lage sein, von Kachelkoordinaten und Pixeln innerhalb einer Kachel in Pixelkoordinaten zu konvertieren, aber das ist es nicht zu schwer.

Schreibe einen Kommentar

Deine E-Mail-Adresse wird nicht veröffentlicht. Erforderliche Felder sind mit * markiert.