[9028] | 1 | #!/usr/bin/env python
|
---|
| 2 | ###############################################################################
|
---|
| 3 | # $Id$
|
---|
| 4 | #
|
---|
| 5 | # Project: GDAL2Tiles, Google Summer of Code 2007 & 2008
|
---|
| 6 | # Global Map Tiles Classes
|
---|
| 7 | # Purpose: Convert a raster into TMS tiles, create KML SuperOverlay EPSG:4326,
|
---|
| 8 | # generate a simple HTML viewers based on Google Maps and OpenLayers
|
---|
| 9 | # Author: Klokan Petr Pridal, klokan at klokan dot cz
|
---|
| 10 | # Web: http://www.klokan.cz/projects/gdal2tiles/
|
---|
| 11 | #
|
---|
| 12 | ###############################################################################
|
---|
| 13 | # Copyright (c) 2008 Klokan Petr Pridal. All rights reserved.
|
---|
| 14 | #
|
---|
| 15 | # Permission is hereby granted, free of charge, to any person obtaining a
|
---|
| 16 | # copy of this software and associated documentation files (the "Software"),
|
---|
| 17 | # to deal in the Software without restriction, including without limitation
|
---|
| 18 | # the rights to use, copy, modify, merge, publish, distribute, sublicense,
|
---|
| 19 | # and/or sell copies of the Software, and to permit persons to whom the
|
---|
| 20 | # Software is furnished to do so, subject to the following conditions:
|
---|
| 21 | #
|
---|
| 22 | # The above copyright notice and this permission notice shall be included
|
---|
| 23 | # in all copies or substantial portions of the Software.
|
---|
| 24 | #
|
---|
| 25 | # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
|
---|
| 26 | # OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
---|
| 27 | # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
|
---|
| 28 | # THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
---|
| 29 | # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
---|
| 30 | # FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
|
---|
| 31 | # DEALINGS IN THE SOFTWARE.
|
---|
| 32 | ###############################################################################
|
---|
| 33 |
|
---|
| 34 | """
|
---|
| 35 | globalmaptiles.py
|
---|
| 36 |
|
---|
| 37 | Global Map Tiles as defined in Tile Map Service (TMS) Profiles
|
---|
| 38 | ==============================================================
|
---|
| 39 |
|
---|
| 40 | Functions necessary for generation of global tiles used on the web.
|
---|
| 41 | It contains classes implementing coordinate conversions for:
|
---|
| 42 |
|
---|
| 43 | - GlobalMercator (based on EPSG:900913 = EPSG:3785)
|
---|
| 44 | for Google Maps, Yahoo Maps, Microsoft Maps compatible tiles
|
---|
| 45 | - GlobalGeodetic (based on EPSG:4326)
|
---|
| 46 | for OpenLayers Base Map and Google Earth compatible tiles
|
---|
| 47 |
|
---|
| 48 | More info at:
|
---|
| 49 |
|
---|
| 50 | http://wiki.osgeo.org/wiki/Tile_Map_Service_Specification
|
---|
| 51 | http://wiki.osgeo.org/wiki/WMS_Tiling_Client_Recommendation
|
---|
| 52 | http://msdn.microsoft.com/en-us/library/bb259689.aspx
|
---|
| 53 | http://code.google.com/apis/maps/documentation/overlays.html#Google_Maps_Coordinates
|
---|
| 54 |
|
---|
| 55 | Created by Klokan Petr Pridal on 2008-07-03.
|
---|
| 56 | Google Summer of Code 2008, project GDAL2Tiles for OSGEO.
|
---|
| 57 |
|
---|
| 58 | In case you use this class in your product, translate it to another language
|
---|
| 59 | or find it usefull for your project please let me know.
|
---|
| 60 | My email: klokan at klokan dot cz.
|
---|
| 61 | I would like to know where it was used.
|
---|
| 62 |
|
---|
| 63 | Class is available under the open-source GDAL license (www.gdal.org).
|
---|
| 64 | """
|
---|
| 65 |
|
---|
| 66 | import math
|
---|
| 67 |
|
---|
| 68 | class GlobalMercator(object):
|
---|
| 69 | """
|
---|
| 70 | TMS Global Mercator Profile
|
---|
| 71 | ---------------------------
|
---|
| 72 |
|
---|
| 73 | Functions necessary for generation of tiles in Spherical Mercator projection,
|
---|
| 74 | EPSG:900913 (EPSG:gOOglE, Google Maps Global Mercator), EPSG:3785, OSGEO:41001.
|
---|
| 75 |
|
---|
| 76 | Such tiles are compatible with Google Maps, Microsoft Virtual Earth, Yahoo Maps,
|
---|
| 77 | UK Ordnance Survey OpenSpace API, ...
|
---|
| 78 | and you can overlay them on top of base maps of those web mapping applications.
|
---|
| 79 |
|
---|
| 80 | Pixel and tile coordinates are in TMS notation (origin [0,0] in bottom-left).
|
---|
| 81 |
|
---|
| 82 | What coordinate conversions do we need for TMS Global Mercator tiles::
|
---|
| 83 |
|
---|
| 84 | LatLon <-> Meters <-> Pixels <-> Tile
|
---|
| 85 |
|
---|
| 86 | WGS84 coordinates Spherical Mercator Pixels in pyramid Tiles in pyramid
|
---|
| 87 | lat/lon XY in metres XY pixels Z zoom XYZ from TMS
|
---|
| 88 | EPSG:4326 EPSG:900913
|
---|
| 89 | .----. --------- -- TMS
|
---|
| 90 | / \ <-> | | <-> /----/ <-> Google
|
---|
| 91 | \ / | | /--------/ QuadTree
|
---|
| 92 | ----- --------- /------------/
|
---|
| 93 | KML, public WebMapService Web Clients TileMapService
|
---|
| 94 |
|
---|
| 95 | What is the coordinate extent of Earth in EPSG:900913?
|
---|
| 96 |
|
---|
| 97 | [-20037508.342789244, -20037508.342789244, 20037508.342789244, 20037508.342789244]
|
---|
| 98 | Constant 20037508.342789244 comes from the circumference of the Earth in meters,
|
---|
| 99 | which is 40 thousand kilometers, the coordinate origin is in the middle of extent.
|
---|
| 100 | In fact you can calculate the constant as: 2 * math.pi * 6378137 / 2.0
|
---|
| 101 | $ echo 180 85 | gdaltransform -s_srs EPSG:4326 -t_srs EPSG:900913
|
---|
| 102 | Polar areas with abs(latitude) bigger then 85.05112878 are clipped off.
|
---|
| 103 |
|
---|
| 104 | What are zoom level constants (pixels/meter) for pyramid with EPSG:900913?
|
---|
| 105 |
|
---|
| 106 | whole region is on top of pyramid (zoom=0) covered by 256x256 pixels tile,
|
---|
| 107 | every lower zoom level resolution is always divided by two
|
---|
| 108 | initialResolution = 20037508.342789244 * 2 / 256 = 156543.03392804062
|
---|
| 109 |
|
---|
| 110 | What is the difference between TMS and Google Maps/QuadTree tile name convention?
|
---|
| 111 |
|
---|
| 112 | The tile raster itself is the same (equal extent, projection, pixel size),
|
---|
| 113 | there is just different identification of the same raster tile.
|
---|
| 114 | Tiles in TMS are counted from [0,0] in the bottom-left corner, id is XYZ.
|
---|
| 115 | Google placed the origin [0,0] to the top-left corner, reference is XYZ.
|
---|
| 116 | Microsoft is referencing tiles by a QuadTree name, defined on the website:
|
---|
| 117 | http://msdn2.microsoft.com/en-us/library/bb259689.aspx
|
---|
| 118 |
|
---|
| 119 | The lat/lon coordinates are using WGS84 datum, yeh?
|
---|
| 120 |
|
---|
| 121 | Yes, all lat/lon we are mentioning should use WGS84 Geodetic Datum.
|
---|
| 122 | Well, the web clients like Google Maps are projecting those coordinates by
|
---|
| 123 | Spherical Mercator, so in fact lat/lon coordinates on sphere are treated as if
|
---|
| 124 | the were on the WGS84 ellipsoid.
|
---|
| 125 |
|
---|
| 126 | From MSDN documentation:
|
---|
| 127 | To simplify the calculations, we use the spherical form of projection, not
|
---|
| 128 | the ellipsoidal form. Since the projection is used only for map display,
|
---|
| 129 | and not for displaying numeric coordinates, we don't need the extra precision
|
---|
| 130 | of an ellipsoidal projection. The spherical projection causes approximately
|
---|
| 131 | 0.33 percent scale distortion in the Y direction, which is not visually noticable.
|
---|
| 132 |
|
---|
| 133 | How do I create a raster in EPSG:900913 and convert coordinates with PROJ.4?
|
---|
| 134 |
|
---|
| 135 | You can use standard GIS tools like gdalwarp, cs2cs or gdaltransform.
|
---|
| 136 | All of the tools supports -t_srs 'epsg:900913'.
|
---|
| 137 |
|
---|
| 138 | For other GIS programs check the exact definition of the projection:
|
---|
| 139 | More info at http://spatialreference.org/ref/user/google-projection/
|
---|
| 140 | The same projection is degined as EPSG:3785. WKT definition is in the official
|
---|
| 141 | EPSG database.
|
---|
| 142 |
|
---|
| 143 | Proj4 Text:
|
---|
| 144 | +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0
|
---|
| 145 | +k=1.0 +units=m +nadgrids=@null +no_defs
|
---|
| 146 |
|
---|
| 147 | Human readable WKT format of EPGS:900913:
|
---|
| 148 | PROJCS["Google Maps Global Mercator",
|
---|
| 149 | GEOGCS["WGS 84",
|
---|
| 150 | DATUM["WGS_1984",
|
---|
| 151 | SPHEROID["WGS 84",6378137,298.2572235630016,
|
---|
| 152 | AUTHORITY["EPSG","7030"]],
|
---|
| 153 | AUTHORITY["EPSG","6326"]],
|
---|
| 154 | PRIMEM["Greenwich",0],
|
---|
| 155 | UNIT["degree",0.0174532925199433],
|
---|
| 156 | AUTHORITY["EPSG","4326"]],
|
---|
| 157 | PROJECTION["Mercator_1SP"],
|
---|
| 158 | PARAMETER["central_meridian",0],
|
---|
| 159 | PARAMETER["scale_factor",1],
|
---|
| 160 | PARAMETER["false_easting",0],
|
---|
| 161 | PARAMETER["false_northing",0],
|
---|
| 162 | UNIT["metre",1,
|
---|
| 163 | AUTHORITY["EPSG","9001"]]]
|
---|
| 164 | """
|
---|
| 165 |
|
---|
| 166 | def __init__(self, tileSize=256):
|
---|
| 167 | "Initialize the TMS Global Mercator pyramid"
|
---|
| 168 | self.tileSize = tileSize
|
---|
| 169 | self.initialResolution = 2 * math.pi * 6378137 / self.tileSize
|
---|
| 170 | # 156543.03392804062 for tileSize 256 pixels
|
---|
| 171 | self.originShift = 2 * math.pi * 6378137 / 2.0
|
---|
| 172 | # 20037508.342789244
|
---|
| 173 |
|
---|
| 174 | def LatLonToMeters(self, lat, lon ):
|
---|
| 175 | "Converts given lat/lon in WGS84 Datum to XY in Spherical Mercator EPSG:900913"
|
---|
| 176 |
|
---|
| 177 | mx = lon * self.originShift / 180.0
|
---|
| 178 | my = math.log( math.tan((90 + lat) * math.pi / 360.0 )) / (math.pi / 180.0)
|
---|
| 179 |
|
---|
| 180 | my = my * self.originShift / 180.0
|
---|
| 181 | return mx, my
|
---|
| 182 |
|
---|
| 183 | def MetersToLatLon(self, mx, my ):
|
---|
| 184 | "Converts XY point from Spherical Mercator EPSG:900913 to lat/lon in WGS84 Datum"
|
---|
| 185 |
|
---|
| 186 | lon = (mx / self.originShift) * 180.0
|
---|
| 187 | lat = (my / self.originShift) * 180.0
|
---|
| 188 |
|
---|
| 189 | lat = 180 / math.pi * (2 * math.atan( math.exp( lat * math.pi / 180.0)) - math.pi / 2.0)
|
---|
| 190 | return lat, lon
|
---|
| 191 |
|
---|
| 192 | def PixelsToMeters(self, px, py, zoom):
|
---|
| 193 | "Converts pixel coordinates in given zoom level of pyramid to EPSG:900913"
|
---|
| 194 |
|
---|
| 195 | res = self.Resolution( zoom )
|
---|
| 196 | mx = px * res - self.originShift
|
---|
| 197 | my = py * res - self.originShift
|
---|
| 198 | return mx, my
|
---|
| 199 |
|
---|
| 200 | def MetersToPixels(self, mx, my, zoom):
|
---|
| 201 | "Converts EPSG:900913 to pyramid pixel coordinates in given zoom level"
|
---|
| 202 |
|
---|
| 203 | res = self.Resolution( zoom )
|
---|
| 204 | px = (mx + self.originShift) / res
|
---|
| 205 | py = (my + self.originShift) / res
|
---|
| 206 | return px, py
|
---|
| 207 |
|
---|
| 208 | def PixelsToTile(self, px, py):
|
---|
| 209 | "Returns a tile covering region in given pixel coordinates"
|
---|
| 210 |
|
---|
| 211 | tx = int( math.ceil( px / float(self.tileSize) ) - 1 )
|
---|
| 212 | ty = int( math.ceil( py / float(self.tileSize) ) - 1 )
|
---|
| 213 | return tx, ty
|
---|
| 214 |
|
---|
| 215 | def PixelsToRaster(self, px, py, zoom):
|
---|
| 216 | "Move the origin of pixel coordinates to top-left corner"
|
---|
| 217 |
|
---|
| 218 | mapSize = self.tileSize << zoom
|
---|
| 219 | return px, mapSize - py
|
---|
| 220 |
|
---|
| 221 | def MetersToTile(self, mx, my, zoom):
|
---|
| 222 | "Returns tile for given mercator coordinates"
|
---|
| 223 |
|
---|
| 224 | px, py = self.MetersToPixels( mx, my, zoom)
|
---|
| 225 | return self.PixelsToTile( px, py)
|
---|
| 226 |
|
---|
| 227 | def TileBounds(self, tx, ty, zoom):
|
---|
| 228 | "Returns bounds of the given tile in EPSG:900913 coordinates"
|
---|
| 229 |
|
---|
| 230 | minx, miny = self.PixelsToMeters( tx*self.tileSize, ty*self.tileSize, zoom )
|
---|
| 231 | maxx, maxy = self.PixelsToMeters( (tx+1)*self.tileSize, (ty+1)*self.tileSize, zoom )
|
---|
| 232 | return ( minx, miny, maxx, maxy )
|
---|
| 233 |
|
---|
| 234 | def TileLatLonBounds(self, tx, ty, zoom ):
|
---|
| 235 | "Returns bounds of the given tile in latutude/longitude using WGS84 datum"
|
---|
| 236 |
|
---|
| 237 | bounds = self.TileBounds( tx, ty, zoom)
|
---|
| 238 | minLat, minLon = self.MetersToLatLon(bounds[0], bounds[1])
|
---|
| 239 | maxLat, maxLon = self.MetersToLatLon(bounds[2], bounds[3])
|
---|
| 240 |
|
---|
| 241 | return ( minLat, minLon, maxLat, maxLon )
|
---|
| 242 |
|
---|
| 243 | def Resolution(self, zoom ):
|
---|
| 244 | "Resolution (meters/pixel) for given zoom level (measured at Equator)"
|
---|
| 245 |
|
---|
| 246 | # return (2 * math.pi * 6378137) / (self.tileSize * 2**zoom)
|
---|
| 247 | return self.initialResolution / (2**zoom)
|
---|
| 248 |
|
---|
| 249 | def ZoomForPixelSize(self, pixelSize ):
|
---|
| 250 | "Maximal scaledown zoom of the pyramid closest to the pixelSize."
|
---|
| 251 |
|
---|
| 252 | for i in range(30):
|
---|
| 253 | if pixelSize > self.Resolution(i):
|
---|
| 254 | return i-1 if i!=0 else 0 # We don't want to scale up
|
---|
| 255 |
|
---|
| 256 | def GoogleTile(self, tx, ty, zoom):
|
---|
| 257 | "Converts TMS tile coordinates to Google Tile coordinates"
|
---|
| 258 |
|
---|
| 259 | # coordinate origin is moved from bottom-left to top-left corner of the extent
|
---|
| 260 | return tx, (2**zoom - 1) - ty
|
---|
| 261 |
|
---|
| 262 | def QuadTree(self, tx, ty, zoom ):
|
---|
| 263 | "Converts TMS tile coordinates to Microsoft QuadTree"
|
---|
| 264 |
|
---|
| 265 | quadKey = ""
|
---|
| 266 | ty = (2**zoom - 1) - ty
|
---|
| 267 | for i in range(zoom, 0, -1):
|
---|
| 268 | digit = 0
|
---|
| 269 | mask = 1 << (i-1)
|
---|
| 270 | if (tx & mask) != 0:
|
---|
| 271 | digit += 1
|
---|
| 272 | if (ty & mask) != 0:
|
---|
| 273 | digit += 2
|
---|
| 274 | quadKey += str(digit)
|
---|
| 275 |
|
---|
| 276 | return quadKey
|
---|
| 277 |
|
---|
| 278 | #---------------------
|
---|
| 279 |
|
---|
| 280 | class GlobalGeodetic(object):
|
---|
| 281 | """
|
---|
| 282 | TMS Global Geodetic Profile
|
---|
| 283 | ---------------------------
|
---|
| 284 |
|
---|
| 285 | Functions necessary for generation of global tiles in Plate Carre projection,
|
---|
| 286 | EPSG:4326, "unprojected profile".
|
---|
| 287 |
|
---|
| 288 | Such tiles are compatible with Google Earth (as any other EPSG:4326 rasters)
|
---|
| 289 | and you can overlay the tiles on top of OpenLayers base map.
|
---|
| 290 |
|
---|
| 291 | Pixel and tile coordinates are in TMS notation (origin [0,0] in bottom-left).
|
---|
| 292 |
|
---|
| 293 | What coordinate conversions do we need for TMS Global Geodetic tiles?
|
---|
| 294 |
|
---|
| 295 | Global Geodetic tiles are using geodetic coordinates (latitude,longitude)
|
---|
| 296 | directly as planar coordinates XY (it is also called Unprojected or Plate
|
---|
| 297 | Carre). We need only scaling to pixel pyramid and cutting to tiles.
|
---|
| 298 | Pyramid has on top level two tiles, so it is not square but rectangle.
|
---|
| 299 | Area [-180,-90,180,90] is scaled to 512x256 pixels.
|
---|
| 300 | TMS has coordinate origin (for pixels and tiles) in bottom-left corner.
|
---|
| 301 | Rasters are in EPSG:4326 and therefore are compatible with Google Earth.
|
---|
| 302 |
|
---|
| 303 | LatLon <-> Pixels <-> Tiles
|
---|
| 304 |
|
---|
| 305 | WGS84 coordinates Pixels in pyramid Tiles in pyramid
|
---|
| 306 | lat/lon XY pixels Z zoom XYZ from TMS
|
---|
| 307 | EPSG:4326
|
---|
| 308 | .----. ----
|
---|
| 309 | / \ <-> /--------/ <-> TMS
|
---|
| 310 | \ / /--------------/
|
---|
| 311 | ----- /--------------------/
|
---|
| 312 | WMS, KML Web Clients, Google Earth TileMapService
|
---|
| 313 | """
|
---|
| 314 |
|
---|
| 315 | def __init__(self, tileSize = 256):
|
---|
| 316 | self.tileSize = tileSize
|
---|
| 317 |
|
---|
| 318 | def LatLonToPixels(self, lat, lon, zoom):
|
---|
| 319 | "Converts lat/lon to pixel coordinates in given zoom of the EPSG:4326 pyramid"
|
---|
| 320 |
|
---|
| 321 | res = 180 / 256.0 / 2**zoom
|
---|
| 322 | px = (180 + lat) / res
|
---|
| 323 | py = (90 + lon) / res
|
---|
| 324 | return px, py
|
---|
| 325 |
|
---|
| 326 | def PixelsToTile(self, px, py):
|
---|
| 327 | "Returns coordinates of the tile covering region in pixel coordinates"
|
---|
| 328 |
|
---|
| 329 | tx = int( math.ceil( px / float(self.tileSize) ) - 1 )
|
---|
| 330 | ty = int( math.ceil( py / float(self.tileSize) ) - 1 )
|
---|
| 331 | return tx, ty
|
---|
| 332 |
|
---|
| 333 | def Resolution(self, zoom ):
|
---|
| 334 | "Resolution (arc/pixel) for given zoom level (measured at Equator)"
|
---|
| 335 |
|
---|
| 336 | return 180 / 256.0 / 2**zoom
|
---|
| 337 | #return 180 / float( 1 << (8+zoom) )
|
---|
| 338 |
|
---|
| 339 | def TileBounds(tx, ty, zoom):
|
---|
| 340 | "Returns bounds of the given tile"
|
---|
| 341 | res = 180 / 256.0 / 2**zoom
|
---|
| 342 | return (
|
---|
| 343 | tx*256*res - 180,
|
---|
| 344 | ty*256*res - 90,
|
---|
| 345 | (tx+1)*256*res - 180,
|
---|
| 346 | (ty+1)*256*res - 90
|
---|
| 347 | )
|
---|
| 348 |
|
---|
| 349 | if __name__ == "__main__":
|
---|
| 350 | import sys, os
|
---|
| 351 |
|
---|
| 352 | def Usage(s = ""):
|
---|
| 353 | print "Usage: globalmaptiles.py [-profile 'mercator'|'geodetic'] zoomlevel lat lon [latmax lonmax]"
|
---|
| 354 | print
|
---|
| 355 | if s:
|
---|
| 356 | print s
|
---|
| 357 | print
|
---|
| 358 | print "This utility prints for given WGS84 lat/lon coordinates (or bounding box) the list of tiles"
|
---|
| 359 | print "covering specified area. Tiles are in the given 'profile' (default is Google Maps 'mercator')"
|
---|
| 360 | print "and in the given pyramid 'zoomlevel'."
|
---|
| 361 | print "For each tile several information is printed including bonding box in EPSG:900913 and WGS84."
|
---|
| 362 | sys.exit(1)
|
---|
| 363 |
|
---|
| 364 | profile = 'mercator'
|
---|
| 365 | zoomlevel = None
|
---|
| 366 | lat, lon, latmax, lonmax = None, None, None, None
|
---|
| 367 | boundingbox = False
|
---|
| 368 |
|
---|
| 369 | argv = sys.argv
|
---|
| 370 | i = 1
|
---|
| 371 | while i < len(argv):
|
---|
| 372 | arg = argv[i]
|
---|
| 373 |
|
---|
| 374 | if arg == '-profile':
|
---|
| 375 | i = i + 1
|
---|
| 376 | profile = argv[i]
|
---|
| 377 |
|
---|
| 378 | if zoomlevel is None:
|
---|
| 379 | zoomlevel = int(argv[i])
|
---|
| 380 | elif lat is None:
|
---|
| 381 | lat = float(argv[i])
|
---|
| 382 | elif lon is None:
|
---|
| 383 | lon = float(argv[i])
|
---|
| 384 | elif latmax is None:
|
---|
| 385 | latmax = float(argv[i])
|
---|
| 386 | elif lonmax is None:
|
---|
| 387 | lonmax = float(argv[i])
|
---|
| 388 | else:
|
---|
| 389 | Usage("ERROR: Too many parameters")
|
---|
| 390 |
|
---|
| 391 | i = i + 1
|
---|
| 392 |
|
---|
| 393 | if profile != 'mercator':
|
---|
| 394 | Usage("ERROR: Sorry, given profile is not implemented yet.")
|
---|
| 395 |
|
---|
| 396 | if zoomlevel == None or lat == None or lon == None:
|
---|
| 397 | Usage("ERROR: Specify at least 'zoomlevel', 'lat' and 'lon'.")
|
---|
| 398 | if latmax is not None and lonmax is None:
|
---|
| 399 | Usage("ERROR: Both 'latmax' and 'lonmax' must be given.")
|
---|
| 400 |
|
---|
| 401 | if latmax != None and lonmax != None:
|
---|
| 402 | if latmax < lat:
|
---|
| 403 | Usage("ERROR: 'latmax' must be bigger then 'lat'")
|
---|
| 404 | if lonmax < lon:
|
---|
| 405 | Usage("ERROR: 'lonmax' must be bigger then 'lon'")
|
---|
| 406 | boundingbox = (lon, lat, lonmax, latmax)
|
---|
| 407 |
|
---|
| 408 | tz = zoomlevel
|
---|
| 409 | mercator = GlobalMercator()
|
---|
| 410 |
|
---|
| 411 | mx, my = mercator.LatLonToMeters( lat, lon )
|
---|
| 412 | print "Spherical Mercator (ESPG:900913) coordinates for lat/lon: "
|
---|
| 413 | print (mx, my)
|
---|
| 414 | tminx, tminy = mercator.MetersToTile( mx, my, tz )
|
---|
| 415 |
|
---|
| 416 | if boundingbox:
|
---|
| 417 | mx, my = mercator.LatLonToMeters( latmax, lonmax )
|
---|
| 418 | print "Spherical Mercator (ESPG:900913) cooridnate for maxlat/maxlon: "
|
---|
| 419 | print (mx, my)
|
---|
| 420 | tmaxx, tmaxy = mercator.MetersToTile( mx, my, tz )
|
---|
| 421 | else:
|
---|
| 422 | tmaxx, tmaxy = tminx, tminy
|
---|
| 423 |
|
---|
| 424 | for ty in range(tminy, tmaxy+1):
|
---|
| 425 | for tx in range(tminx, tmaxx+1):
|
---|
| 426 | tilefilename = "%s/%s/%s" % (tz, tx, ty)
|
---|
| 427 | print tilefilename, "( TileMapService: z / x / y )"
|
---|
| 428 |
|
---|
| 429 | gx, gy = mercator.GoogleTile(tx, ty, tz)
|
---|
| 430 | print "\tGoogle:", gx, gy
|
---|
| 431 | quadkey = mercator.QuadTree(tx, ty, tz)
|
---|
| 432 | print "\tQuadkey:", quadkey, '(',int(quadkey, 4),')'
|
---|
| 433 | bounds = mercator.TileBounds( tx, ty, tz)
|
---|
| 434 | print
|
---|
| 435 | print "\tEPSG:900913 Extent: ", bounds
|
---|
| 436 | wgsbounds = mercator.TileLatLonBounds( tx, ty, tz)
|
---|
| 437 | print "\tWGS84 Extent:", wgsbounds
|
---|
| 438 | print "\tgdalwarp -ts 256 256 -te %s %s %s %s %s %s_%s_%s.tif" % (
|
---|
| 439 | bounds[0], bounds[1], bounds[2], bounds[3], "<your-raster-file-in-epsg900913.ext>", tz, tx, ty)
|
---|
| 440 | print
|
---|