Jag försöker georeferens en uppsättning bildkakel nedladdade från en Google Maps-webbplats. Baseras endast på varje bild ”s Google Maps nätreferens och zoomnivå, hur kan jag beräkna de verkliga koordinaterna?

Sidan vid http://www.maptiler.org/google-maps-coordinates-tile-bounds-projection visar koordinaterna och Google Maps rutnätreferens för webbsäljare.

Till exempel är mittplattan x=49, y=30, z=6 hörnkoordinater 10644926.307106785, 626172.1357121654, 11271098.44281895, 1252344.27142432:

ange bildbeskrivning här

Är det möjligt att beräkna hörnkoordinaterna enbart baserat på nätreferensen och zoomnivån (49, 30, 6)?

Detta liknar Hur georefererar jag en webcommercatorplatta korrekt med gdal? men jag behöver en helt programmatisk lösning utan att behöva leta upp något online. (Eller om jag behöver leta upp värden online måste det vara programmatiskt.).

Redigera: den här skärmdumpen visar placeringen av en provplatta (zoomnivå 11, x = 1593, y = 933) . Hur skulle jag georferens om den brickan?

ange bildbeskrivning här

Kommentarer

  • Kommer det som en riktig sida? Ingen omgivning?
  • @ FelixIP ja allt jag vet är nätreferens och zoomnivå för varje enskild JPG. Jag ’ har lagt till en skärmdump för att visa detta. PS Jag inser att detta är lite tvivelaktigt, men ibland måste det göras – plus att det ’ är en intressant övning!
  • Jag tycker att du borde kunna få hörnet kordinerar med bara zoomnivåer och rutnätsreferens, så vitt jag vet förblir nätreferenssystemet konstant, så det borde finnas ett sätt att göra detta. Eftersom zoomnivå 0 har världen i 1 kakel och zoom 1 i 4 och så vidare bör det vara möjligt att få koordinaterna.
  • @StephenLead kanske försök att dela upp världen i fyra delar i varje zoom ( iteration) och få koordinaterna på något sätt. Det ’ måste vara en iterativ funktion men enligt min förståelse.
  • får du göra detta? det vill säga, tillåter Google Maps termer detta?

Svar

Denna sida http://www.maptiler.org/google-maps-coordinates-tile-bounds-projection tillhandahåller algoritmer, skrivna i python, för att beräkna brickgränser i EPSG: 900913-koordinater och i latutude / longitud med WGS84-referens.

Klassen GlobalMercator i detta skript http://www.maptiler.org/google-maps-coordinates-tile-bounds-projection / globalmaptiles .py innehåller två metoder, TileBounds() och TileLatLonBounds() för att göra detta.

En modifierad version av skriptet för att visa kordkoordinater för kakel:

 #!/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  

Användning: tilebounds xTile yTile zoom.

Till exempel utgången för brickan x=49, y=30, z=6 är som följer:

$./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) 

En sida kan laddas ned med url http://mt.google.com/vt?x=xTile&y=yTile&z=zoom, men detta är förbjudet att komma åt det direkt.

Programvaran skrevs ursprungligen av Klokan Petr Přidal.

Jag hoppas att det här kan hjälpa!

Kommentarer

  • Välkommen till GIS Stack Exchange, och tack för detta. Mitt krav har passerat eftersom jag ’ inte längre arbetar med det här projektet, men denna information ska förhoppningsvis vara användbar för någon annan i framtiden

Svar

Kasta bara in det här som en del av diskussionen om brickor, cachade brickor, georeferenser osv.

Tegelplaner från vad Jag vet att jag inte använder ett XY-georeferenssystem i sig …

När du till exempel läser in ESRI-baskartan i QGIS, skapar den en struktur för mappstapelmappen så här:

ange bildbeskrivning här

När du tittar på toppnivåbilder har du begärde brickorna och de lagras i de olika mapparna.

Här är innehållet i mappen f / d:

ange bildbeskrivning här

Observera ingen georeferensinformation i filen (världsfil etc.)

När du tar med den här filen till MS måla, till exempel kan du manipulera bilden och spara den igen:

ange bildbeskrivning här

…så när du visar kartan igen i QGIS – som återger brickorna – kan du se din ”redigerade” sida:

ange bildbeskrivning här

Så jag antar att allt jag säger är att om du hade en uppsättning brickor från Google Earth, kan du i teorin återskapa brickan strukturera och lägg helt enkelt brickorna i dessa mappar, och de skulle rita i QGIS …

Det är själva programvaran som vet var plattorna ska räknas baserat på cachingschemat och mappstrukturen … detta går tillbaka till mina dagar då jag undervisade ArcGIS Server-kurser.

Du kan sedan exportera dem från QGIS med georeferensinformation och använda dem någon annanstans.

Men då kan du också bara ladda google maps i QGIS och exportera bilden som en TIF med en världsfil …

(Återigen, inte riktigt ett svar, bara lite diskussion …)

Kommentarer

  • ledsen för det sena svaret. För att vara ärlig kan jag ’ inte ens komma ihåg varför jag frågade det här, eftersom projektet är över:) Jag minns dock att jag behövde återskapa vilken voodoo som programvaran gör när den vet hur man placerar bilder baserat på mappstrukturen

Svar

Påträffade detta problem nyligen kan du använda koden i wiki.openstreetmap.org

Här är Java-koden:

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

Svar

För den underliggande matematiska metoden, se mitt svar på en relaterad fråga här: beräkning av avgränsningsruta från känd centrumkoordinat och zoom

Du måste kunna konvertera från kakelkoordinater och pixlar inom en kakel till pixelkoordinater, men det inte för hårt.

Lämna ett svar

Din e-postadress kommer inte publiceras. Obligatoriska fält är märkta *