[23] | 1 | #!/usr/bin/env python |
---|
| 2 | #****************************************************************************** |
---|
| 3 | # $Id: gdal2tiles.py 18194 2009-12-06 20:07:45Z rouault $ |
---|
| 4 | # |
---|
| 5 | # Project: Google Summer of Code 2007, 2008 (http://code.google.com/soc/) |
---|
| 6 | # Support: BRGM (http://www.brgm.fr) |
---|
| 7 | # Purpose: Convert a raster into TMS (Tile Map Service) tiles in a directory. |
---|
| 8 | # - generate Google Earth metadata (KML SuperOverlay) |
---|
| 9 | # - generate simple HTML viewer based on Google Maps and OpenLayers |
---|
| 10 | # - support of global tiles (Spherical Mercator) for compatibility |
---|
| 11 | # with interactive web maps a la Google Maps |
---|
| 12 | # Author: Klokan Petr Pridal, klokan at klokan dot cz |
---|
| 13 | # Web: http://www.klokan.cz/projects/gdal2tiles/ |
---|
| 14 | # GUI: http://www.maptiler.org/ |
---|
| 15 | # |
---|
| 16 | ############################################################################### |
---|
| 17 | # Copyright (c) 2008, Klokan Petr Pridal |
---|
| 18 | # |
---|
| 19 | # Permission is hereby granted, free of charge, to any person obtaining a |
---|
| 20 | # copy of this software and associated documentation files (the "Software"), |
---|
| 21 | # to deal in the Software without restriction, including without limitation |
---|
| 22 | # the rights to use, copy, modify, merge, publish, distribute, sublicense, |
---|
| 23 | # and/or sell copies of the Software, and to permit persons to whom the |
---|
| 24 | # Software is furnished to do so, subject to the following conditions: |
---|
| 25 | # |
---|
| 26 | # The above copyright notice and this permission notice shall be included |
---|
| 27 | # in all copies or substantial portions of the Software. |
---|
| 28 | # |
---|
| 29 | # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS |
---|
| 30 | # OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, |
---|
| 31 | # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL |
---|
| 32 | # THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER |
---|
| 33 | # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING |
---|
| 34 | # FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER |
---|
| 35 | # DEALINGS IN THE SOFTWARE. |
---|
| 36 | #****************************************************************************** |
---|
| 37 | |
---|
| 38 | from osgeo import gdal |
---|
| 39 | from osgeo import osr |
---|
| 40 | |
---|
| 41 | import sys |
---|
| 42 | import os |
---|
| 43 | import math |
---|
| 44 | |
---|
| 45 | try: |
---|
| 46 | from PIL import Image |
---|
| 47 | import numpy |
---|
| 48 | import osgeo.gdal_array as gdalarray |
---|
| 49 | except: |
---|
| 50 | # 'antialias' resampling is not available |
---|
| 51 | pass |
---|
| 52 | |
---|
| 53 | __version__ = "$Id: gdal2tiles.py 18194 2009-12-06 20:07:45Z rouault $" |
---|
| 54 | |
---|
| 55 | resampling_list = ('average','near','bilinear','cubic','cubicspline','lanczos','antialias') |
---|
| 56 | profile_list = ('mercator','geodetic','raster') #,'zoomify') |
---|
| 57 | webviewer_list = ('all','google','openlayers','none') |
---|
| 58 | |
---|
| 59 | # ============================================================================= |
---|
| 60 | # ============================================================================= |
---|
| 61 | # ============================================================================= |
---|
| 62 | |
---|
| 63 | __doc__globalmaptiles = """ |
---|
| 64 | globalmaptiles.py |
---|
| 65 | |
---|
| 66 | Global Map Tiles as defined in Tile Map Service (TMS) Profiles |
---|
| 67 | ============================================================== |
---|
| 68 | |
---|
| 69 | Functions necessary for generation of global tiles used on the web. |
---|
| 70 | It contains classes implementing coordinate conversions for: |
---|
| 71 | |
---|
| 72 | - GlobalMercator (based on EPSG:900913 = EPSG:3785) |
---|
| 73 | for Google Maps, Yahoo Maps, Microsoft Maps compatible tiles |
---|
| 74 | - GlobalGeodetic (based on EPSG:4326) |
---|
| 75 | for OpenLayers Base Map and Google Earth compatible tiles |
---|
| 76 | |
---|
| 77 | More info at: |
---|
| 78 | |
---|
| 79 | http://wiki.osgeo.org/wiki/Tile_Map_Service_Specification |
---|
| 80 | http://wiki.osgeo.org/wiki/WMS_Tiling_Client_Recommendation |
---|
| 81 | http://msdn.microsoft.com/en-us/library/bb259689.aspx |
---|
| 82 | http://code.google.com/apis/maps/documentation/overlays.html#Google_Maps_Coordinates |
---|
| 83 | |
---|
| 84 | Created by Klokan Petr Pridal on 2008-07-03. |
---|
| 85 | Google Summer of Code 2008, project GDAL2Tiles for OSGEO. |
---|
| 86 | |
---|
| 87 | In case you use this class in your product, translate it to another language |
---|
| 88 | or find it usefull for your project please let me know. |
---|
| 89 | My email: klokan at klokan dot cz. |
---|
| 90 | I would like to know where it was used. |
---|
| 91 | |
---|
| 92 | Class is available under the open-source GDAL license (www.gdal.org). |
---|
| 93 | """ |
---|
| 94 | |
---|
| 95 | import math |
---|
| 96 | |
---|
| 97 | MAXZOOMLEVEL = 32 |
---|
| 98 | |
---|
| 99 | class GlobalMercator(object): |
---|
| 100 | """ |
---|
| 101 | TMS Global Mercator Profile |
---|
| 102 | --------------------------- |
---|
| 103 | |
---|
| 104 | Functions necessary for generation of tiles in Spherical Mercator projection, |
---|
| 105 | EPSG:900913 (EPSG:gOOglE, Google Maps Global Mercator), EPSG:3785, OSGEO:41001. |
---|
| 106 | |
---|
| 107 | Such tiles are compatible with Google Maps, Microsoft Virtual Earth, Yahoo Maps, |
---|
| 108 | UK Ordnance Survey OpenSpace API, ... |
---|
| 109 | and you can overlay them on top of base maps of those web mapping applications. |
---|
| 110 | |
---|
| 111 | Pixel and tile coordinates are in TMS notation (origin [0,0] in bottom-left). |
---|
| 112 | |
---|
| 113 | What coordinate conversions do we need for TMS Global Mercator tiles:: |
---|
| 114 | |
---|
| 115 | LatLon <-> Meters <-> Pixels <-> Tile |
---|
| 116 | |
---|
| 117 | WGS84 coordinates Spherical Mercator Pixels in pyramid Tiles in pyramid |
---|
| 118 | lat/lon XY in metres XY pixels Z zoom XYZ from TMS |
---|
| 119 | EPSG:4326 EPSG:900913 |
---|
| 120 | .----. --------- -- TMS |
---|
| 121 | / \ <-> | | <-> /----/ <-> Google |
---|
| 122 | \ / | | /--------/ QuadTree |
---|
| 123 | ----- --------- /------------/ |
---|
| 124 | KML, public WebMapService Web Clients TileMapService |
---|
| 125 | |
---|
| 126 | What is the coordinate extent of Earth in EPSG:900913? |
---|
| 127 | |
---|
| 128 | [-20037508.342789244, -20037508.342789244, 20037508.342789244, 20037508.342789244] |
---|
| 129 | Constant 20037508.342789244 comes from the circumference of the Earth in meters, |
---|
| 130 | which is 40 thousand kilometers, the coordinate origin is in the middle of extent. |
---|
| 131 | In fact you can calculate the constant as: 2 * math.pi * 6378137 / 2.0 |
---|
| 132 | $ echo 180 85 | gdaltransform -s_srs EPSG:4326 -t_srs EPSG:900913 |
---|
| 133 | Polar areas with abs(latitude) bigger then 85.05112878 are clipped off. |
---|
| 134 | |
---|
| 135 | What are zoom level constants (pixels/meter) for pyramid with EPSG:900913? |
---|
| 136 | |
---|
| 137 | whole region is on top of pyramid (zoom=0) covered by 256x256 pixels tile, |
---|
| 138 | every lower zoom level resolution is always divided by two |
---|
| 139 | initialResolution = 20037508.342789244 * 2 / 256 = 156543.03392804062 |
---|
| 140 | |
---|
| 141 | What is the difference between TMS and Google Maps/QuadTree tile name convention? |
---|
| 142 | |
---|
| 143 | The tile raster itself is the same (equal extent, projection, pixel size), |
---|
| 144 | there is just different identification of the same raster tile. |
---|
| 145 | Tiles in TMS are counted from [0,0] in the bottom-left corner, id is XYZ. |
---|
| 146 | Google placed the origin [0,0] to the top-left corner, reference is XYZ. |
---|
| 147 | Microsoft is referencing tiles by a QuadTree name, defined on the website: |
---|
| 148 | http://msdn2.microsoft.com/en-us/library/bb259689.aspx |
---|
| 149 | |
---|
| 150 | The lat/lon coordinates are using WGS84 datum, yeh? |
---|
| 151 | |
---|
| 152 | Yes, all lat/lon we are mentioning should use WGS84 Geodetic Datum. |
---|
| 153 | Well, the web clients like Google Maps are projecting those coordinates by |
---|
| 154 | Spherical Mercator, so in fact lat/lon coordinates on sphere are treated as if |
---|
| 155 | the were on the WGS84 ellipsoid. |
---|
| 156 | |
---|
| 157 | From MSDN documentation: |
---|
| 158 | To simplify the calculations, we use the spherical form of projection, not |
---|
| 159 | the ellipsoidal form. Since the projection is used only for map display, |
---|
| 160 | and not for displaying numeric coordinates, we don't need the extra precision |
---|
| 161 | of an ellipsoidal projection. The spherical projection causes approximately |
---|
| 162 | 0.33 percent scale distortion in the Y direction, which is not visually noticable. |
---|
| 163 | |
---|
| 164 | How do I create a raster in EPSG:900913 and convert coordinates with PROJ.4? |
---|
| 165 | |
---|
| 166 | You can use standard GIS tools like gdalwarp, cs2cs or gdaltransform. |
---|
| 167 | All of the tools supports -t_srs 'epsg:900913'. |
---|
| 168 | |
---|
| 169 | For other GIS programs check the exact definition of the projection: |
---|
| 170 | More info at http://spatialreference.org/ref/user/google-projection/ |
---|
| 171 | The same projection is degined as EPSG:3785. WKT definition is in the official |
---|
| 172 | EPSG database. |
---|
| 173 | |
---|
| 174 | Proj4 Text: |
---|
| 175 | +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 |
---|
| 176 | +k=1.0 +units=m +nadgrids=@null +no_defs |
---|
| 177 | |
---|
| 178 | Human readable WKT format of EPGS:900913: |
---|
| 179 | PROJCS["Google Maps Global Mercator", |
---|
| 180 | GEOGCS["WGS 84", |
---|
| 181 | DATUM["WGS_1984", |
---|
| 182 | SPHEROID["WGS 84",6378137,298.257223563, |
---|
| 183 | AUTHORITY["EPSG","7030"]], |
---|
| 184 | AUTHORITY["EPSG","6326"]], |
---|
| 185 | PRIMEM["Greenwich",0], |
---|
| 186 | UNIT["degree",0.0174532925199433], |
---|
| 187 | AUTHORITY["EPSG","4326"]], |
---|
| 188 | PROJECTION["Mercator_1SP"], |
---|
| 189 | PARAMETER["central_meridian",0], |
---|
| 190 | PARAMETER["scale_factor",1], |
---|
| 191 | PARAMETER["false_easting",0], |
---|
| 192 | PARAMETER["false_northing",0], |
---|
| 193 | UNIT["metre",1, |
---|
| 194 | AUTHORITY["EPSG","9001"]]] |
---|
| 195 | """ |
---|
| 196 | |
---|
| 197 | def __init__(self, tileSize=256): |
---|
| 198 | "Initialize the TMS Global Mercator pyramid" |
---|
| 199 | self.tileSize = tileSize |
---|
| 200 | self.initialResolution = 2 * math.pi * 6378137 / self.tileSize |
---|
| 201 | # 156543.03392804062 for tileSize 256 pixels |
---|
| 202 | self.originShift = 2 * math.pi * 6378137 / 2.0 |
---|
| 203 | # 20037508.342789244 |
---|
| 204 | |
---|
| 205 | def LatLonToMeters(self, lat, lon ): |
---|
| 206 | "Converts given lat/lon in WGS84 Datum to XY in Spherical Mercator EPSG:900913" |
---|
| 207 | |
---|
| 208 | mx = lon * self.originShift / 180.0 |
---|
| 209 | my = math.log( math.tan((90 + lat) * math.pi / 360.0 )) / (math.pi / 180.0) |
---|
| 210 | |
---|
| 211 | my = my * self.originShift / 180.0 |
---|
| 212 | return mx, my |
---|
| 213 | |
---|
| 214 | def MetersToLatLon(self, mx, my ): |
---|
| 215 | "Converts XY point from Spherical Mercator EPSG:900913 to lat/lon in WGS84 Datum" |
---|
| 216 | |
---|
| 217 | lon = (mx / self.originShift) * 180.0 |
---|
| 218 | lat = (my / self.originShift) * 180.0 |
---|
| 219 | |
---|
| 220 | lat = 180 / math.pi * (2 * math.atan( math.exp( lat * math.pi / 180.0)) - math.pi / 2.0) |
---|
| 221 | return lat, lon |
---|
| 222 | |
---|
| 223 | def PixelsToMeters(self, px, py, zoom): |
---|
| 224 | "Converts pixel coordinates in given zoom level of pyramid to EPSG:900913" |
---|
| 225 | |
---|
| 226 | res = self.Resolution( zoom ) |
---|
| 227 | mx = px * res - self.originShift |
---|
| 228 | my = py * res - self.originShift |
---|
| 229 | return mx, my |
---|
| 230 | |
---|
| 231 | def MetersToPixels(self, mx, my, zoom): |
---|
| 232 | "Converts EPSG:900913 to pyramid pixel coordinates in given zoom level" |
---|
| 233 | |
---|
| 234 | res = self.Resolution( zoom ) |
---|
| 235 | px = (mx + self.originShift) / res |
---|
| 236 | py = (my + self.originShift) / res |
---|
| 237 | return px, py |
---|
| 238 | |
---|
| 239 | def PixelsToTile(self, px, py): |
---|
| 240 | "Returns a tile covering region in given pixel coordinates" |
---|
| 241 | |
---|
| 242 | tx = int( math.ceil( px / float(self.tileSize) ) - 1 ) |
---|
| 243 | ty = int( math.ceil( py / float(self.tileSize) ) - 1 ) |
---|
| 244 | return tx, ty |
---|
| 245 | |
---|
| 246 | def PixelsToRaster(self, px, py, zoom): |
---|
| 247 | "Move the origin of pixel coordinates to top-left corner" |
---|
| 248 | |
---|
| 249 | mapSize = self.tileSize << zoom |
---|
| 250 | return px, mapSize - py |
---|
| 251 | |
---|
| 252 | def MetersToTile(self, mx, my, zoom): |
---|
| 253 | "Returns tile for given mercator coordinates" |
---|
| 254 | |
---|
| 255 | px, py = self.MetersToPixels( mx, my, zoom) |
---|
| 256 | return self.PixelsToTile( px, py) |
---|
| 257 | |
---|
| 258 | def TileBounds(self, tx, ty, zoom): |
---|
| 259 | "Returns bounds of the given tile in EPSG:900913 coordinates" |
---|
| 260 | |
---|
| 261 | minx, miny = self.PixelsToMeters( tx*self.tileSize, ty*self.tileSize, zoom ) |
---|
| 262 | maxx, maxy = self.PixelsToMeters( (tx+1)*self.tileSize, (ty+1)*self.tileSize, zoom ) |
---|
| 263 | return ( minx, miny, maxx, maxy ) |
---|
| 264 | |
---|
| 265 | def TileLatLonBounds(self, tx, ty, zoom ): |
---|
| 266 | "Returns bounds of the given tile in latutude/longitude using WGS84 datum" |
---|
| 267 | |
---|
| 268 | bounds = self.TileBounds( tx, ty, zoom) |
---|
| 269 | minLat, minLon = self.MetersToLatLon(bounds[0], bounds[1]) |
---|
| 270 | maxLat, maxLon = self.MetersToLatLon(bounds[2], bounds[3]) |
---|
| 271 | |
---|
| 272 | return ( minLat, minLon, maxLat, maxLon ) |
---|
| 273 | |
---|
| 274 | def Resolution(self, zoom ): |
---|
| 275 | "Resolution (meters/pixel) for given zoom level (measured at Equator)" |
---|
| 276 | |
---|
| 277 | # return (2 * math.pi * 6378137) / (self.tileSize * 2**zoom) |
---|
| 278 | return self.initialResolution / (2**zoom) |
---|
| 279 | |
---|
| 280 | def ZoomForPixelSize(self, pixelSize ): |
---|
| 281 | "Maximal scaledown zoom of the pyramid closest to the pixelSize." |
---|
| 282 | |
---|
| 283 | for i in range(MAXZOOMLEVEL): |
---|
| 284 | if pixelSize > self.Resolution(i): |
---|
| 285 | if i!=0: |
---|
| 286 | return i-1 |
---|
| 287 | else: |
---|
| 288 | return 0 # We don't want to scale up |
---|
| 289 | |
---|
| 290 | def GoogleTile(self, tx, ty, zoom): |
---|
| 291 | "Converts TMS tile coordinates to Google Tile coordinates" |
---|
| 292 | |
---|
| 293 | # coordinate origin is moved from bottom-left to top-left corner of the extent |
---|
| 294 | return tx, (2**zoom - 1) - ty |
---|
| 295 | |
---|
| 296 | def QuadTree(self, tx, ty, zoom ): |
---|
| 297 | "Converts TMS tile coordinates to Microsoft QuadTree" |
---|
| 298 | |
---|
| 299 | quadKey = "" |
---|
| 300 | ty = (2**zoom - 1) - ty |
---|
| 301 | for i in range(zoom, 0, -1): |
---|
| 302 | digit = 0 |
---|
| 303 | mask = 1 << (i-1) |
---|
| 304 | if (tx & mask) != 0: |
---|
| 305 | digit += 1 |
---|
| 306 | if (ty & mask) != 0: |
---|
| 307 | digit += 2 |
---|
| 308 | quadKey += str(digit) |
---|
| 309 | |
---|
| 310 | return quadKey |
---|
| 311 | |
---|
| 312 | #--------------------- |
---|
| 313 | |
---|
| 314 | class GlobalGeodetic(object): |
---|
| 315 | """ |
---|
| 316 | TMS Global Geodetic Profile |
---|
| 317 | --------------------------- |
---|
| 318 | |
---|
| 319 | Functions necessary for generation of global tiles in Plate Carre projection, |
---|
| 320 | EPSG:4326, "unprojected profile". |
---|
| 321 | |
---|
| 322 | Such tiles are compatible with Google Earth (as any other EPSG:4326 rasters) |
---|
| 323 | and you can overlay the tiles on top of OpenLayers base map. |
---|
| 324 | |
---|
| 325 | Pixel and tile coordinates are in TMS notation (origin [0,0] in bottom-left). |
---|
| 326 | |
---|
| 327 | What coordinate conversions do we need for TMS Global Geodetic tiles? |
---|
| 328 | |
---|
| 329 | Global Geodetic tiles are using geodetic coordinates (latitude,longitude) |
---|
| 330 | directly as planar coordinates XY (it is also called Unprojected or Plate |
---|
| 331 | Carre). We need only scaling to pixel pyramid and cutting to tiles. |
---|
| 332 | Pyramid has on top level two tiles, so it is not square but rectangle. |
---|
| 333 | Area [-180,-90,180,90] is scaled to 512x256 pixels. |
---|
| 334 | TMS has coordinate origin (for pixels and tiles) in bottom-left corner. |
---|
| 335 | Rasters are in EPSG:4326 and therefore are compatible with Google Earth. |
---|
| 336 | |
---|
| 337 | LatLon <-> Pixels <-> Tiles |
---|
| 338 | |
---|
| 339 | WGS84 coordinates Pixels in pyramid Tiles in pyramid |
---|
| 340 | lat/lon XY pixels Z zoom XYZ from TMS |
---|
| 341 | EPSG:4326 |
---|
| 342 | .----. ---- |
---|
| 343 | / \ <-> /--------/ <-> TMS |
---|
| 344 | \ / /--------------/ |
---|
| 345 | ----- /--------------------/ |
---|
| 346 | WMS, KML Web Clients, Google Earth TileMapService |
---|
| 347 | """ |
---|
| 348 | |
---|
| 349 | def __init__(self, tileSize = 256): |
---|
| 350 | self.tileSize = tileSize |
---|
| 351 | |
---|
| 352 | def LatLonToPixels(self, lat, lon, zoom): |
---|
| 353 | "Converts lat/lon to pixel coordinates in given zoom of the EPSG:4326 pyramid" |
---|
| 354 | |
---|
| 355 | res = 180.0 / self.tileSize / 2**zoom |
---|
| 356 | px = (180 + lat) / res |
---|
| 357 | py = (90 + lon) / res |
---|
| 358 | return px, py |
---|
| 359 | |
---|
| 360 | def PixelsToTile(self, px, py): |
---|
| 361 | "Returns coordinates of the tile covering region in pixel coordinates" |
---|
| 362 | |
---|
| 363 | tx = int( math.ceil( px / float(self.tileSize) ) - 1 ) |
---|
| 364 | ty = int( math.ceil( py / float(self.tileSize) ) - 1 ) |
---|
| 365 | return tx, ty |
---|
| 366 | |
---|
| 367 | def LatLonToTile(self, lat, lon, zoom): |
---|
| 368 | "Returns the tile for zoom which covers given lat/lon coordinates" |
---|
| 369 | |
---|
| 370 | px, py = self.LatLonToPixels( lat, lon, zoom) |
---|
| 371 | return self.PixelsToTile(px,py) |
---|
| 372 | |
---|
| 373 | def Resolution(self, zoom ): |
---|
| 374 | "Resolution (arc/pixel) for given zoom level (measured at Equator)" |
---|
| 375 | |
---|
| 376 | return 180.0 / self.tileSize / 2**zoom |
---|
| 377 | #return 180 / float( 1 << (8+zoom) ) |
---|
| 378 | |
---|
| 379 | def ZoomForPixelSize(self, pixelSize ): |
---|
| 380 | "Maximal scaledown zoom of the pyramid closest to the pixelSize." |
---|
| 381 | |
---|
| 382 | for i in range(MAXZOOMLEVEL): |
---|
| 383 | if pixelSize > self.Resolution(i): |
---|
| 384 | if i!=0: |
---|
| 385 | return i-1 |
---|
| 386 | else: |
---|
| 387 | return 0 # We don't want to scale up |
---|
| 388 | |
---|
| 389 | def TileBounds(self, tx, ty, zoom): |
---|
| 390 | "Returns bounds of the given tile" |
---|
| 391 | res = 180.0 / self.tileSize / 2**zoom |
---|
| 392 | return ( |
---|
| 393 | tx*self.tileSize*res - 180, |
---|
| 394 | ty*self.tileSize*res - 90, |
---|
| 395 | (tx+1)*self.tileSize*res - 180, |
---|
| 396 | (ty+1)*self.tileSize*res - 90 |
---|
| 397 | ) |
---|
| 398 | |
---|
| 399 | def TileLatLonBounds(self, tx, ty, zoom): |
---|
| 400 | "Returns bounds of the given tile in the SWNE form" |
---|
| 401 | b = self.TileBounds(tx, ty, zoom) |
---|
| 402 | return (b[1],b[0],b[3],b[2]) |
---|
| 403 | |
---|
| 404 | #--------------------- |
---|
| 405 | # TODO: Finish Zoomify implemtentation!!! |
---|
| 406 | class Zoomify(object): |
---|
| 407 | """ |
---|
| 408 | Tiles compatible with the Zoomify viewer |
---|
| 409 | ---------------------------------------- |
---|
| 410 | """ |
---|
| 411 | |
---|
| 412 | def __init__(self, width, height, tilesize = 256, tileformat='jpg'): |
---|
| 413 | """Initialization of the Zoomify tile tree""" |
---|
| 414 | |
---|
| 415 | self.tilesize = tilesize |
---|
| 416 | self.tileformat = tileformat |
---|
| 417 | imagesize = (width, height) |
---|
| 418 | tiles = ( math.ceil( width / tilesize ), math.ceil( height / tilesize ) ) |
---|
| 419 | |
---|
| 420 | # Size (in tiles) for each tier of pyramid. |
---|
| 421 | self.tierSizeInTiles = [] |
---|
| 422 | self.tierSizeInTiles.push( tiles ) |
---|
| 423 | |
---|
| 424 | # Image size in pixels for each pyramid tierself |
---|
| 425 | self.tierImageSize = [] |
---|
| 426 | self.tierImageSize.append( imagesize ); |
---|
| 427 | |
---|
| 428 | while (imagesize[0] > tilesize or imageSize[1] > tilesize ): |
---|
| 429 | imagesize = (math.floor( imagesize[0] / 2 ), math.floor( imagesize[1] / 2) ) |
---|
| 430 | tiles = ( math.ceil( imagesize[0] / tilesize ), math.ceil( imagesize[1] / tilesize ) ) |
---|
| 431 | self.tierSizeInTiles.append( tiles ) |
---|
| 432 | self.tierImageSize.append( imagesize ) |
---|
| 433 | |
---|
| 434 | self.tierSizeInTiles.reverse() |
---|
| 435 | self.tierImageSize.reverse() |
---|
| 436 | |
---|
| 437 | # Depth of the Zoomify pyramid, number of tiers (zoom levels) |
---|
| 438 | self.numberOfTiers = len(self.tierSizeInTiles) |
---|
| 439 | |
---|
| 440 | # Number of tiles up to the given tier of pyramid. |
---|
| 441 | self.tileCountUpToTier = [] |
---|
| 442 | self.tileCountUpToTier[0] = 0 |
---|
| 443 | for i in range(1, self.numberOfTiers+1): |
---|
| 444 | self.tileCountUpToTier.append( |
---|
| 445 | self.tierSizeInTiles[i-1][0] * self.tierSizeInTiles[i-1][1] + self.tileCountUpToTier[i-1] |
---|
| 446 | ) |
---|
| 447 | |
---|
| 448 | def tilefilename(self, x, y, z): |
---|
| 449 | """Returns filename for tile with given coordinates""" |
---|
| 450 | |
---|
| 451 | tileIndex = x + y * self.tierSizeInTiles[z][0] + self.tileCountUpToTier[z] |
---|
| 452 | return os.path.join("TileGroup%.0f" % math.floor( tileIndex / 256 ), |
---|
| 453 | "%s-%s-%s.%s" % ( z, x, y, self.tileformat)) |
---|
| 454 | |
---|
| 455 | # ============================================================================= |
---|
| 456 | # ============================================================================= |
---|
| 457 | # ============================================================================= |
---|
| 458 | |
---|
| 459 | class GDAL2Tiles(object): |
---|
| 460 | |
---|
| 461 | # ------------------------------------------------------------------------- |
---|
| 462 | def process(self): |
---|
| 463 | """The main processing function, runs all the main steps of processing""" |
---|
| 464 | |
---|
| 465 | # Opening and preprocessing of the input file |
---|
| 466 | self.open_input() |
---|
| 467 | |
---|
| 468 | # Generation of main metadata files and HTML viewers |
---|
| 469 | self.generate_metadata() |
---|
| 470 | |
---|
| 471 | # Generation of the lowest tiles |
---|
| 472 | self.generate_base_tiles() |
---|
| 473 | |
---|
| 474 | # Generation of the overview tiles (higher in the pyramid) |
---|
| 475 | self.generate_overview_tiles() |
---|
| 476 | |
---|
| 477 | # ------------------------------------------------------------------------- |
---|
| 478 | def error(self, msg, details = "" ): |
---|
| 479 | """Print an error message and stop the processing""" |
---|
| 480 | |
---|
| 481 | if details: |
---|
| 482 | self.parser.error(msg + "\n\n" + details) |
---|
| 483 | else: |
---|
| 484 | self.parser.error(msg) |
---|
| 485 | |
---|
| 486 | # ------------------------------------------------------------------------- |
---|
| 487 | def progressbar(self, complete = 0.0): |
---|
| 488 | """Print progressbar for float value 0..1""" |
---|
| 489 | |
---|
| 490 | gdal.TermProgress_nocb(complete) |
---|
| 491 | |
---|
| 492 | # ------------------------------------------------------------------------- |
---|
| 493 | def stop(self): |
---|
| 494 | """Stop the rendering immediately""" |
---|
| 495 | self.stopped = True |
---|
| 496 | |
---|
| 497 | # ------------------------------------------------------------------------- |
---|
| 498 | def __init__(self, arguments ): |
---|
| 499 | """Constructor function - initialization""" |
---|
| 500 | |
---|
| 501 | self.stopped = False |
---|
| 502 | self.input = None |
---|
| 503 | self.output = None |
---|
| 504 | |
---|
| 505 | # Tile format |
---|
| 506 | self.tilesize = 256 |
---|
| 507 | self.tiledriver = 'PNG' |
---|
| 508 | self.tileext = 'png' |
---|
| 509 | |
---|
| 510 | # Should we read bigger window of the input raster and scale it down? |
---|
| 511 | # Note: Modified leter by open_input() |
---|
| 512 | # Not for 'near' resampling |
---|
| 513 | # Not for Wavelet based drivers (JPEG2000, ECW, MrSID) |
---|
| 514 | # Not for 'raster' profile |
---|
| 515 | self.scaledquery = True |
---|
| 516 | # How big should be query window be for scaling down |
---|
| 517 | # Later on reset according the chosen resampling algorightm |
---|
| 518 | self.querysize = 4 * self.tilesize |
---|
| 519 | |
---|
| 520 | # Should we use Read on the input file for generating overview tiles? |
---|
| 521 | # Note: Modified later by open_input() |
---|
| 522 | # Otherwise the overview tiles are generated from existing underlying tiles |
---|
| 523 | self.overviewquery = False |
---|
| 524 | |
---|
| 525 | # RUN THE ARGUMENT PARSER: |
---|
| 526 | |
---|
| 527 | self.optparse_init() |
---|
| 528 | self.options, self.args = self.parser.parse_args(args=arguments) |
---|
| 529 | if not self.args: |
---|
| 530 | self.error("No input file specified") |
---|
| 531 | |
---|
| 532 | # POSTPROCESSING OF PARSED ARGUMENTS: |
---|
| 533 | |
---|
| 534 | # Workaround for old versions of GDAL |
---|
| 535 | try: |
---|
| 536 | if (self.options.verbose and self.options.resampling == 'near') or gdal.TermProgress_nocb: |
---|
| 537 | pass |
---|
| 538 | except: |
---|
| 539 | self.error("This version of GDAL is not supported. Please upgrade to 1.6+.") |
---|
| 540 | #,"You can try run crippled version of gdal2tiles with parameters: -v -r 'near'") |
---|
| 541 | |
---|
| 542 | # Is output directory the last argument? |
---|
| 543 | |
---|
| 544 | # Test output directory, if it doesn't exist |
---|
| 545 | if os.path.isdir(self.args[-1]) or ( len(self.args) > 1 and not os.path.exists(self.args[-1])): |
---|
| 546 | self.output = self.args[-1] |
---|
| 547 | self.args = self.args[:-1] |
---|
| 548 | |
---|
| 549 | # More files on the input not directly supported yet |
---|
| 550 | |
---|
| 551 | if (len(self.args) > 1): |
---|
| 552 | self.error("Processing of several input files is not supported.", |
---|
| 553 | """Please first use a tool like gdal_vrtmerge.py or gdal_merge.py on the files: |
---|
| 554 | gdal_vrtmerge.py -o merged.vrt %s""" % " ".join(self.args)) |
---|
| 555 | # TODO: Call functions from gdal_vrtmerge.py directly |
---|
| 556 | |
---|
| 557 | self.input = self.args[0] |
---|
| 558 | |
---|
| 559 | # Default values for not given options |
---|
| 560 | |
---|
| 561 | if not self.output: |
---|
| 562 | # Directory with input filename without extension in actual directory |
---|
| 563 | self.output = os.path.splitext(os.path.basename( self.input ))[0] |
---|
| 564 | |
---|
| 565 | if not self.options.title: |
---|
| 566 | self.options.title = os.path.basename( self.input ) |
---|
| 567 | |
---|
| 568 | if self.options.url and not self.options.url.endswith('/'): |
---|
| 569 | self.options.url += '/' |
---|
| 570 | if self.options.url: |
---|
| 571 | self.options.url += os.path.basename( self.output ) + '/' |
---|
| 572 | |
---|
| 573 | # Supported options |
---|
| 574 | |
---|
| 575 | self.resampling = None |
---|
| 576 | |
---|
| 577 | if self.options.resampling == 'average': |
---|
| 578 | try: |
---|
| 579 | if gdal.RegenerateOverview: |
---|
| 580 | pass |
---|
| 581 | except: |
---|
| 582 | self.error("'average' resampling algorithm is not available.", "Please use -r 'near' argument or upgrade to newer version of GDAL.") |
---|
| 583 | |
---|
| 584 | elif self.options.resampling == 'antialias': |
---|
| 585 | try: |
---|
| 586 | if numpy: |
---|
| 587 | pass |
---|
| 588 | except: |
---|
| 589 | self.error("'antialias' resampling algorithm is not available.", "Install PIL (Python Imaging Library) and numpy.") |
---|
| 590 | |
---|
| 591 | elif self.options.resampling == 'near': |
---|
| 592 | self.resampling = gdal.GRA_NearestNeighbour |
---|
| 593 | self.querysize = self.tilesize |
---|
| 594 | |
---|
| 595 | elif self.options.resampling == 'bilinear': |
---|
| 596 | self.resampling = gdal.GRA_Bilinear |
---|
| 597 | self.querysize = self.tilesize * 2 |
---|
| 598 | |
---|
| 599 | elif self.options.resampling == 'cubic': |
---|
| 600 | self.resampling = gdal.GRA_Cubic |
---|
| 601 | |
---|
| 602 | elif self.options.resampling == 'cubicspline': |
---|
| 603 | self.resampling = gdal.GRA_CubicSpline |
---|
| 604 | |
---|
| 605 | elif self.options.resampling == 'lanczos': |
---|
| 606 | self.resampling = gdal.GRA_Lanczos |
---|
| 607 | |
---|
| 608 | # User specified zoom levels |
---|
| 609 | self.tminz = None |
---|
| 610 | self.tmaxz = None |
---|
| 611 | if self.options.zoom: |
---|
| 612 | minmax = self.options.zoom.split('-',1) |
---|
| 613 | minmax.extend(['']) |
---|
| 614 | min, max = minmax[:2] |
---|
| 615 | self.tminz = int(min) |
---|
| 616 | if max: |
---|
| 617 | self.tmaxz = int(max) |
---|
| 618 | else: |
---|
| 619 | self.tmaxz = int(min) |
---|
| 620 | |
---|
| 621 | # KML generation |
---|
| 622 | self.kml = self.options.kml |
---|
| 623 | |
---|
| 624 | # Output the results |
---|
| 625 | |
---|
| 626 | if self.options.verbose: |
---|
| 627 | print("Options:", self.options) |
---|
| 628 | print("Input:", self.input) |
---|
| 629 | print("Output:", self.output) |
---|
| 630 | print("Cache: %s MB" % (gdal.GetCacheMax() / 1024 / 1024)) |
---|
| 631 | print('') |
---|
| 632 | |
---|
| 633 | # ------------------------------------------------------------------------- |
---|
| 634 | def optparse_init(self): |
---|
| 635 | """Prepare the option parser for input (argv)""" |
---|
| 636 | |
---|
| 637 | from optparse import OptionParser, OptionGroup |
---|
| 638 | usage = "Usage: %prog [options] input_file(s) [output]" |
---|
| 639 | p = OptionParser(usage, version="%prog "+ __version__) |
---|
| 640 | p.add_option("-p", "--profile", dest='profile', type='choice', choices=profile_list, |
---|
| 641 | help="Tile cutting profile (%s) - default 'mercator' (Google Maps compatible)" % ",".join(profile_list)) |
---|
| 642 | p.add_option("-r", "--resampling", dest="resampling", type='choice', choices=resampling_list, |
---|
| 643 | help="Resampling method (%s) - default 'average'" % ",".join(resampling_list)) |
---|
| 644 | p.add_option('-s', '--s_srs', dest="s_srs", metavar="SRS", |
---|
| 645 | help="The spatial reference system used for the source input data") |
---|
| 646 | p.add_option('-z', '--zoom', dest="zoom", |
---|
| 647 | help="Zoom levels to render (format:'2-5' or '10').") |
---|
| 648 | p.add_option('-e', '--resume', dest="resume", action="store_true", |
---|
| 649 | help="Resume mode. Generate only missing files.") |
---|
| 650 | p.add_option('-a', '--srcnodata', dest="srcnodata", metavar="NODATA", |
---|
| 651 | help="NODATA transparency value to assign to the input data") |
---|
| 652 | p.add_option("-v", "--verbose", |
---|
| 653 | action="store_true", dest="verbose", |
---|
| 654 | help="Print status messages to stdout") |
---|
| 655 | |
---|
| 656 | # KML options |
---|
| 657 | g = OptionGroup(p, "KML (Google Earth) options", "Options for generated Google Earth SuperOverlay metadata") |
---|
| 658 | g.add_option("-k", "--force-kml", dest='kml', action="store_true", |
---|
| 659 | help="Generate KML for Google Earth - default for 'geodetic' profile and 'raster' in EPSG:4326. For a dataset with different projection use with caution!") |
---|
| 660 | g.add_option("-n", "--no-kml", dest='kml', action="store_false", |
---|
| 661 | help="Avoid automatic generation of KML files for EPSG:4326") |
---|
| 662 | g.add_option("-u", "--url", dest='url', |
---|
| 663 | help="URL address where the generated tiles are going to be published") |
---|
| 664 | p.add_option_group(g) |
---|
| 665 | |
---|
| 666 | # HTML options |
---|
| 667 | g = OptionGroup(p, "Web viewer options", "Options for generated HTML viewers a la Google Maps") |
---|
| 668 | g.add_option("-w", "--webviewer", dest='webviewer', type='choice', choices=webviewer_list, |
---|
| 669 | help="Web viewer to generate (%s) - default 'all'" % ",".join(webviewer_list)) |
---|
| 670 | g.add_option("-t", "--title", dest='title', |
---|
| 671 | help="Title of the map") |
---|
| 672 | g.add_option("-c", "--copyright", dest='copyright', |
---|
| 673 | help="Copyright for the map") |
---|
| 674 | g.add_option("-g", "--googlekey", dest='googlekey', |
---|
| 675 | help="Google Maps API key from http://code.google.com/apis/maps/signup.html") |
---|
| 676 | g.add_option("-y", "--yahookey", dest='yahookey', |
---|
| 677 | help="Yahoo Application ID from http://developer.yahoo.com/wsregapp/") |
---|
| 678 | p.add_option_group(g) |
---|
| 679 | |
---|
| 680 | # TODO: MapFile + TileIndexes per zoom level for efficient MapServer WMS |
---|
| 681 | #g = OptionGroup(p, "WMS MapServer metadata", "Options for generated mapfile and tileindexes for MapServer") |
---|
| 682 | #g.add_option("-i", "--tileindex", dest='wms', action="store_true" |
---|
| 683 | # help="Generate tileindex and mapfile for MapServer (WMS)") |
---|
| 684 | # p.add_option_group(g) |
---|
| 685 | |
---|
| 686 | p.set_defaults(verbose=False, profile="mercator", kml=False, url='', |
---|
| 687 | webviewer='all', copyright='', resampling='average', resume=False, |
---|
| 688 | googlekey='INSERT_YOUR_KEY_HERE', yahookey='INSERT_YOUR_YAHOO_APP_ID_HERE') |
---|
| 689 | |
---|
| 690 | self.parser = p |
---|
| 691 | |
---|
| 692 | # ------------------------------------------------------------------------- |
---|
| 693 | def open_input(self): |
---|
| 694 | """Initialization of the input raster, reprojection if necessary""" |
---|
| 695 | |
---|
| 696 | gdal.AllRegister() |
---|
| 697 | |
---|
| 698 | # Initialize necessary GDAL drivers |
---|
| 699 | |
---|
| 700 | self.out_drv = gdal.GetDriverByName( self.tiledriver ) |
---|
| 701 | self.mem_drv = gdal.GetDriverByName( 'MEM' ) |
---|
| 702 | |
---|
| 703 | if not self.out_drv: |
---|
| 704 | raise Exception("The '%s' driver was not found, is it available in this GDAL build?", self.tiledriver) |
---|
| 705 | if not self.mem_drv: |
---|
| 706 | raise Exception("The 'MEM' driver was not found, is it available in this GDAL build?") |
---|
| 707 | |
---|
| 708 | # Open the input file |
---|
| 709 | |
---|
| 710 | if self.input: |
---|
| 711 | self.in_ds = gdal.Open(self.input, gdal.GA_ReadOnly) |
---|
| 712 | else: |
---|
| 713 | raise Exception("No input file was specified") |
---|
| 714 | |
---|
| 715 | if self.options.verbose: |
---|
| 716 | print("Input file:", "( %sP x %sL - %s bands)" % (self.in_ds.RasterXSize, self.in_ds.RasterYSize, self.in_ds.RasterCount)) |
---|
| 717 | |
---|
| 718 | if not self.in_ds: |
---|
| 719 | # Note: GDAL prints the ERROR message too |
---|
| 720 | self.error("It is not possible to open the input file '%s'." % self.input ) |
---|
| 721 | |
---|
| 722 | # Read metadata from the input file |
---|
| 723 | if self.in_ds.RasterCount == 0: |
---|
| 724 | self.error( "Input file '%s' has no raster band" % self.input ) |
---|
| 725 | |
---|
| 726 | if self.in_ds.GetRasterBand(1).GetRasterColorTable(): |
---|
| 727 | # TODO: Process directly paletted dataset by generating VRT in memory |
---|
| 728 | self.error( "Please convert this file to RGB/RGBA and run gdal2tiles on the result.", |
---|
| 729 | """From paletted file you can create RGBA file (temp.vrt) by: |
---|
| 730 | gdal_translate -of vrt -expand rgba %s temp.vrt |
---|
| 731 | then run: |
---|
| 732 | gdal2tiles temp.vrt""" % self.input ) |
---|
| 733 | |
---|
| 734 | # Get NODATA value |
---|
| 735 | self.in_nodata = [] |
---|
| 736 | for i in range(1, self.in_ds.RasterCount+1): |
---|
| 737 | if self.in_ds.GetRasterBand(i).GetNoDataValue() != None: |
---|
| 738 | self.in_nodata.append( self.in_ds.GetRasterBand(i).GetNoDataValue() ) |
---|
| 739 | if self.options.srcnodata: |
---|
| 740 | nds = list(map( float, self.options.srcnodata.split(','))) |
---|
| 741 | if len(nds) < self.in_ds.RasterCount: |
---|
| 742 | self.in_nodata = (nds * self.in_ds.RasterCount)[:self.in_ds.RasterCount] |
---|
| 743 | else: |
---|
| 744 | self.in_nodata = nds |
---|
| 745 | |
---|
| 746 | if self.options.verbose: |
---|
| 747 | print("NODATA: %s" % self.in_nodata) |
---|
| 748 | |
---|
| 749 | # |
---|
| 750 | # Here we should have RGBA input dataset opened in self.in_ds |
---|
| 751 | # |
---|
| 752 | |
---|
| 753 | if self.options.verbose: |
---|
| 754 | print("Preprocessed file:", "( %sP x %sL - %s bands)" % (self.in_ds.RasterXSize, self.in_ds.RasterYSize, self.in_ds.RasterCount)) |
---|
| 755 | |
---|
| 756 | # Spatial Reference System of the input raster |
---|
| 757 | |
---|
| 758 | |
---|
| 759 | self.in_srs = None |
---|
| 760 | |
---|
| 761 | if self.options.s_srs: |
---|
| 762 | self.in_srs = osr.SpatialReference() |
---|
| 763 | self.in_srs.SetFromUserInput(self.options.s_srs) |
---|
| 764 | self.in_srs_wkt = self.in_srs.ExportToWkt() |
---|
| 765 | else: |
---|
| 766 | self.in_srs_wkt = self.in_ds.GetProjection() |
---|
| 767 | if not self.in_srs_wkt and self.in_ds.GetGCPCount() != 0: |
---|
| 768 | self.in_srs_wkt = self.in_ds.GetGCPProjection() |
---|
| 769 | if self.in_srs_wkt: |
---|
| 770 | self.in_srs = osr.SpatialReference() |
---|
| 771 | self.in_srs.ImportFromWkt(self.in_srs_wkt) |
---|
| 772 | #elif self.options.profile != 'raster': |
---|
| 773 | # self.error("There is no spatial reference system info included in the input file.","You should run gdal2tiles with --s_srs EPSG:XXXX or similar.") |
---|
| 774 | |
---|
| 775 | # Spatial Reference System of tiles |
---|
| 776 | |
---|
| 777 | self.out_srs = osr.SpatialReference() |
---|
| 778 | |
---|
| 779 | if self.options.profile == 'mercator': |
---|
| 780 | self.out_srs.ImportFromEPSG(900913) |
---|
| 781 | elif self.options.profile == 'geodetic': |
---|
| 782 | self.out_srs.ImportFromEPSG(4326) |
---|
| 783 | else: |
---|
| 784 | self.out_srs = self.in_srs |
---|
| 785 | |
---|
| 786 | # Are the reference systems the same? Reproject if necessary. |
---|
| 787 | |
---|
| 788 | self.out_ds = None |
---|
| 789 | |
---|
| 790 | if self.options.profile in ('mercator', 'geodetic'): |
---|
| 791 | |
---|
| 792 | if (self.in_ds.GetGeoTransform() == (0.0, 1.0, 0.0, 0.0, 0.0, 1.0)) and (self.in_ds.GetGCPCount() == 0): |
---|
| 793 | self.error("There is no georeference - neither affine transformation (worldfile) nor GCPs. You can generate only 'raster' profile tiles.", |
---|
| 794 | "Either gdal2tiles with parameter -p 'raster' or use another GIS software for georeference e.g. gdal_transform -gcp / -a_ullr / -a_srs") |
---|
| 795 | |
---|
| 796 | if self.in_srs: |
---|
| 797 | |
---|
| 798 | if (self.in_srs.ExportToProj4() != self.out_srs.ExportToProj4()) or (self.in_ds.GetGCPCount() != 0): |
---|
| 799 | |
---|
| 800 | # Generation of VRT dataset in tile projection, default 'nearest neighbour' warping |
---|
| 801 | self.out_ds = gdal.AutoCreateWarpedVRT( self.in_ds, self.in_srs_wkt, self.out_srs.ExportToWkt() ) |
---|
| 802 | |
---|
| 803 | # TODO: HIGH PRIORITY: Correction of AutoCreateWarpedVRT according the max zoomlevel for correct direct warping!!! |
---|
| 804 | |
---|
| 805 | if self.options.verbose: |
---|
| 806 | print("Warping of the raster by AutoCreateWarpedVRT (result saved into 'tiles.vrt')") |
---|
| 807 | self.out_ds.GetDriver().CreateCopy("tiles.vrt", self.out_ds) |
---|
| 808 | |
---|
| 809 | # Note: self.in_srs and self.in_srs_wkt contain still the non-warped reference system!!! |
---|
| 810 | |
---|
| 811 | # Correction of AutoCreateWarpedVRT for NODATA values |
---|
| 812 | if self.in_nodata != []: |
---|
| 813 | import tempfile |
---|
| 814 | tempfilename = tempfile.mktemp('-gdal2tiles.vrt') |
---|
| 815 | self.out_ds.GetDriver().CreateCopy(tempfilename, self.out_ds) |
---|
| 816 | # open as a text file |
---|
| 817 | s = open(tempfilename).read() |
---|
| 818 | # Add the warping options |
---|
| 819 | s = s.replace("""<GDALWarpOptions>""","""<GDALWarpOptions> |
---|
| 820 | <Option name="INIT_DEST">NO_DATA</Option> |
---|
| 821 | <Option name="UNIFIED_SRC_NODATA">YES</Option>""") |
---|
| 822 | # replace BandMapping tag for NODATA bands.... |
---|
| 823 | for i in range(len(self.in_nodata)): |
---|
| 824 | s = s.replace("""<BandMapping src="%i" dst="%i"/>""" % ((i+1),(i+1)),"""<BandMapping src="%i" dst="%i"> |
---|
| 825 | <SrcNoDataReal>%i</SrcNoDataReal> |
---|
| 826 | <SrcNoDataImag>0</SrcNoDataImag> |
---|
| 827 | <DstNoDataReal>%i</DstNoDataReal> |
---|
| 828 | <DstNoDataImag>0</DstNoDataImag> |
---|
| 829 | </BandMapping>""" % ((i+1), (i+1), self.in_nodata[i], self.in_nodata[i])) # Or rewrite to white by: , 255 )) |
---|
| 830 | # save the corrected VRT |
---|
| 831 | open(tempfilename,"w").write(s) |
---|
| 832 | # open by GDAL as self.out_ds |
---|
| 833 | self.out_ds = gdal.Open(tempfilename) #, gdal.GA_ReadOnly) |
---|
| 834 | # delete the temporary file |
---|
| 835 | os.unlink(tempfilename) |
---|
| 836 | |
---|
| 837 | # set NODATA_VALUE metadata |
---|
| 838 | self.out_ds.SetMetadataItem('NODATA_VALUES','%i %i %i' % (self.in_nodata[0],self.in_nodata[1],self.in_nodata[2])) |
---|
| 839 | |
---|
| 840 | if self.options.verbose: |
---|
| 841 | print("Modified warping result saved into 'tiles1.vrt'") |
---|
| 842 | open("tiles1.vrt","w").write(s) |
---|
| 843 | |
---|
| 844 | # ----------------------------------- |
---|
| 845 | # Correction of AutoCreateWarpedVRT for Mono (1 band) and RGB (3 bands) files without NODATA: |
---|
| 846 | # equivalent of gdalwarp -dstalpha |
---|
| 847 | if self.in_nodata == [] and self.out_ds.RasterCount in [1,3]: |
---|
| 848 | import tempfile |
---|
| 849 | tempfilename = tempfile.mktemp('-gdal2tiles.vrt') |
---|
| 850 | self.out_ds.GetDriver().CreateCopy(tempfilename, self.out_ds) |
---|
| 851 | # open as a text file |
---|
| 852 | s = open(tempfilename).read() |
---|
| 853 | # Add the warping options |
---|
| 854 | s = s.replace("""<BlockXSize>""","""<VRTRasterBand dataType="Byte" band="%i" subClass="VRTWarpedRasterBand"> |
---|
| 855 | <ColorInterp>Alpha</ColorInterp> |
---|
| 856 | </VRTRasterBand> |
---|
| 857 | <BlockXSize>""" % (self.out_ds.RasterCount + 1)) |
---|
| 858 | s = s.replace("""</GDALWarpOptions>""", """<DstAlphaBand>%i</DstAlphaBand> |
---|
| 859 | </GDALWarpOptions>""" % (self.out_ds.RasterCount + 1)) |
---|
| 860 | s = s.replace("""</WorkingDataType>""", """</WorkingDataType> |
---|
| 861 | <Option name="INIT_DEST">0</Option>""") |
---|
| 862 | # save the corrected VRT |
---|
| 863 | open(tempfilename,"w").write(s) |
---|
| 864 | # open by GDAL as self.out_ds |
---|
| 865 | self.out_ds = gdal.Open(tempfilename) #, gdal.GA_ReadOnly) |
---|
| 866 | # delete the temporary file |
---|
| 867 | os.unlink(tempfilename) |
---|
| 868 | |
---|
| 869 | if self.options.verbose: |
---|
| 870 | print("Modified -dstalpha warping result saved into 'tiles1.vrt'") |
---|
| 871 | open("tiles1.vrt","w").write(s) |
---|
| 872 | s = ''' |
---|
| 873 | ''' |
---|
| 874 | |
---|
| 875 | else: |
---|
| 876 | self.error("Input file has unknown SRS.", "Use --s_srs ESPG:xyz (or similar) to provide source reference system." ) |
---|
| 877 | |
---|
| 878 | if self.out_ds and self.options.verbose: |
---|
| 879 | print("Projected file:", "tiles.vrt", "( %sP x %sL - %s bands)" % (self.out_ds.RasterXSize, self.out_ds.RasterYSize, self.out_ds.RasterCount)) |
---|
| 880 | |
---|
| 881 | if not self.out_ds: |
---|
| 882 | self.out_ds = self.in_ds |
---|
| 883 | |
---|
| 884 | # |
---|
| 885 | # Here we should have a raster (out_ds) in the correct Spatial Reference system |
---|
| 886 | # |
---|
| 887 | |
---|
| 888 | # Get alpha band (either directly or from NODATA value) |
---|
| 889 | self.alphaband = self.out_ds.GetRasterBand(1).GetMaskBand() |
---|
| 890 | if (self.alphaband.GetMaskFlags() & gdal.GMF_ALPHA) or self.out_ds.RasterCount==4 or self.out_ds.RasterCount==2: |
---|
| 891 | # TODO: Better test for alpha band in the dataset |
---|
| 892 | self.dataBandsCount = self.out_ds.RasterCount - 1 |
---|
| 893 | else: |
---|
| 894 | self.dataBandsCount = self.out_ds.RasterCount |
---|
| 895 | |
---|
| 896 | # KML test |
---|
| 897 | self.isepsg4326 = False |
---|
| 898 | srs4326 = osr.SpatialReference() |
---|
| 899 | srs4326.ImportFromEPSG(4326) |
---|
| 900 | if self.out_srs and srs4326.ExportToProj4() == self.out_srs.ExportToProj4(): |
---|
| 901 | self.kml = True |
---|
| 902 | self.isepsg4326 = True |
---|
| 903 | if self.options.verbose: |
---|
| 904 | print("KML autotest OK!") |
---|
| 905 | |
---|
| 906 | # Read the georeference |
---|
| 907 | |
---|
| 908 | self.out_gt = self.out_ds.GetGeoTransform() |
---|
| 909 | |
---|
| 910 | #originX, originY = self.out_gt[0], self.out_gt[3] |
---|
| 911 | #pixelSize = self.out_gt[1] # = self.out_gt[5] |
---|
| 912 | |
---|
| 913 | # Test the size of the pixel |
---|
| 914 | |
---|
| 915 | # MAPTILER - COMMENTED |
---|
| 916 | #if self.out_gt[1] != (-1 * self.out_gt[5]) and self.options.profile != 'raster': |
---|
| 917 | # TODO: Process corectly coordinates with are have swichted Y axis (display in OpenLayers too) |
---|
| 918 | #self.error("Size of the pixel in the output differ for X and Y axes.") |
---|
| 919 | |
---|
| 920 | # Report error in case rotation/skew is in geotransform (possible only in 'raster' profile) |
---|
| 921 | if (self.out_gt[2], self.out_gt[4]) != (0,0): |
---|
| 922 | self.error("Georeference of the raster contains rotation or skew. Such raster is not supported. Please use gdalwarp first.") |
---|
| 923 | # TODO: Do the warping in this case automaticaly |
---|
| 924 | |
---|
| 925 | # |
---|
| 926 | # Here we expect: pixel is square, no rotation on the raster |
---|
| 927 | # |
---|
| 928 | |
---|
| 929 | # Output Bounds - coordinates in the output SRS |
---|
| 930 | self.ominx = self.out_gt[0] |
---|
| 931 | self.omaxx = self.out_gt[0]+self.out_ds.RasterXSize*self.out_gt[1] |
---|
| 932 | self.omaxy = self.out_gt[3] |
---|
| 933 | self.ominy = self.out_gt[3]-self.out_ds.RasterYSize*self.out_gt[1] |
---|
| 934 | # Note: maybe round(x, 14) to avoid the gdal_translate behaviour, when 0 becomes -1e-15 |
---|
| 935 | |
---|
| 936 | if self.options.verbose: |
---|
| 937 | print("Bounds (output srs):", round(self.ominx, 13), self.ominy, self.omaxx, self.omaxy) |
---|
| 938 | |
---|
| 939 | # |
---|
| 940 | # Calculating ranges for tiles in different zoom levels |
---|
| 941 | # |
---|
| 942 | |
---|
| 943 | if self.options.profile == 'mercator': |
---|
| 944 | |
---|
| 945 | self.mercator = GlobalMercator() # from globalmaptiles.py |
---|
| 946 | |
---|
| 947 | # Function which generates SWNE in LatLong for given tile |
---|
| 948 | self.tileswne = self.mercator.TileLatLonBounds |
---|
| 949 | |
---|
| 950 | # Generate table with min max tile coordinates for all zoomlevels |
---|
| 951 | self.tminmax = list(range(0,32)) |
---|
| 952 | for tz in range(0, 32): |
---|
| 953 | tminx, tminy = self.mercator.MetersToTile( self.ominx, self.ominy, tz ) |
---|
| 954 | tmaxx, tmaxy = self.mercator.MetersToTile( self.omaxx, self.omaxy, tz ) |
---|
| 955 | # crop tiles extending world limits (+-180,+-90) |
---|
| 956 | tminx, tminy = max(0, tminx), max(0, tminy) |
---|
| 957 | tmaxx, tmaxy = min(2**tz-1, tmaxx), min(2**tz-1, tmaxy) |
---|
| 958 | self.tminmax[tz] = (tminx, tminy, tmaxx, tmaxy) |
---|
| 959 | |
---|
| 960 | # TODO: Maps crossing 180E (Alaska?) |
---|
| 961 | |
---|
| 962 | # Get the minimal zoom level (map covers area equivalent to one tile) |
---|
| 963 | if self.tminz == None: |
---|
| 964 | self.tminz = self.mercator.ZoomForPixelSize( self.out_gt[1] * max( self.out_ds.RasterXSize, self.out_ds.RasterYSize) / float(self.tilesize) ) |
---|
| 965 | |
---|
| 966 | # Get the maximal zoom level (closest possible zoom level up on the resolution of raster) |
---|
| 967 | if self.tmaxz == None: |
---|
| 968 | self.tmaxz = self.mercator.ZoomForPixelSize( self.out_gt[1] ) |
---|
| 969 | |
---|
| 970 | if self.options.verbose: |
---|
| 971 | print("Bounds (latlong):", self.mercator.MetersToLatLon( self.ominx, self.ominy), self.mercator.MetersToLatLon( self.omaxx, self.omaxy)) |
---|
| 972 | print('MinZoomLevel:', self.tminz) |
---|
| 973 | print("MaxZoomLevel:", self.tmaxz, "(", self.mercator.Resolution( self.tmaxz ),")") |
---|
| 974 | |
---|
| 975 | if self.options.profile == 'geodetic': |
---|
| 976 | |
---|
| 977 | self.geodetic = GlobalGeodetic() # from globalmaptiles.py |
---|
| 978 | |
---|
| 979 | # Function which generates SWNE in LatLong for given tile |
---|
| 980 | self.tileswne = self.geodetic.TileLatLonBounds |
---|
| 981 | |
---|
| 982 | # Generate table with min max tile coordinates for all zoomlevels |
---|
| 983 | self.tminmax = list(range(0,32)) |
---|
| 984 | for tz in range(0, 32): |
---|
| 985 | tminx, tminy = self.geodetic.LatLonToTile( self.ominx, self.ominy, tz ) |
---|
| 986 | tmaxx, tmaxy = self.geodetic.LatLonToTile( self.omaxx, self.omaxy, tz ) |
---|
| 987 | # crop tiles extending world limits (+-180,+-90) |
---|
| 988 | tminx, tminy = max(0, tminx), max(0, tminy) |
---|
| 989 | tmaxx, tmaxy = min(2**(tz+1)-1, tmaxx), min(2**tz-1, tmaxy) |
---|
| 990 | self.tminmax[tz] = (tminx, tminy, tmaxx, tmaxy) |
---|
| 991 | |
---|
| 992 | # TODO: Maps crossing 180E (Alaska?) |
---|
| 993 | |
---|
| 994 | # Get the maximal zoom level (closest possible zoom level up on the resolution of raster) |
---|
| 995 | if self.tminz == None: |
---|
| 996 | self.tminz = self.geodetic.ZoomForPixelSize( self.out_gt[1] * max( self.out_ds.RasterXSize, self.out_ds.RasterYSize) / float(self.tilesize) ) |
---|
| 997 | |
---|
| 998 | # Get the maximal zoom level (closest possible zoom level up on the resolution of raster) |
---|
| 999 | if self.tmaxz == None: |
---|
| 1000 | self.tmaxz = self.geodetic.ZoomForPixelSize( self.out_gt[1] ) |
---|
| 1001 | |
---|
| 1002 | if self.options.verbose: |
---|
| 1003 | print("Bounds (latlong):", self.ominx, self.ominy, self.omaxx, self.omaxy) |
---|
| 1004 | |
---|
| 1005 | if self.options.profile == 'raster': |
---|
| 1006 | |
---|
| 1007 | log2 = lambda x: math.log10(x) / math.log10(2) # log2 (base 2 logarithm) |
---|
| 1008 | |
---|
| 1009 | self.nativezoom = int(max( math.ceil(log2(self.out_ds.RasterXSize/float(self.tilesize))), |
---|
| 1010 | math.ceil(log2(self.out_ds.RasterYSize/float(self.tilesize))))) |
---|
| 1011 | |
---|
| 1012 | if self.options.verbose: |
---|
| 1013 | print("Native zoom of the raster:", self.nativezoom) |
---|
| 1014 | |
---|
| 1015 | # Get the minimal zoom level (whole raster in one tile) |
---|
| 1016 | if self.tminz == None: |
---|
| 1017 | self.tminz = 0 |
---|
| 1018 | |
---|
| 1019 | # Get the maximal zoom level (native resolution of the raster) |
---|
| 1020 | if self.tmaxz == None: |
---|
| 1021 | self.tmaxz = self.nativezoom |
---|
| 1022 | |
---|
| 1023 | # Generate table with min max tile coordinates for all zoomlevels |
---|
| 1024 | self.tminmax = list(range(0, self.tmaxz+1)) |
---|
| 1025 | self.tsize = list(range(0, self.tmaxz+1)) |
---|
| 1026 | for tz in range(0, self.tmaxz+1): |
---|
| 1027 | tsize = 2.0**(self.nativezoom-tz)*self.tilesize |
---|
| 1028 | tminx, tminy = 0, 0 |
---|
| 1029 | tmaxx = int(math.ceil( self.out_ds.RasterXSize / tsize )) - 1 |
---|
| 1030 | tmaxy = int(math.ceil( self.out_ds.RasterYSize / tsize )) - 1 |
---|
| 1031 | self.tsize[tz] = math.ceil(tsize) |
---|
| 1032 | self.tminmax[tz] = (tminx, tminy, tmaxx, tmaxy) |
---|
| 1033 | |
---|
| 1034 | # Function which generates SWNE in LatLong for given tile |
---|
| 1035 | if self.kml and self.in_srs_wkt: |
---|
| 1036 | self.ct = osr.CoordinateTransformation(self.in_srs, srs4326) |
---|
| 1037 | |
---|
| 1038 | def rastertileswne(x,y,z): |
---|
| 1039 | pixelsizex = (2**(self.tmaxz-z) * self.out_gt[1]) # X-pixel size in level |
---|
| 1040 | pixelsizey = (2**(self.tmaxz-z) * self.out_gt[1]) # Y-pixel size in level (usually -1*pixelsizex) |
---|
| 1041 | west = self.out_gt[0] + x*self.tilesize*pixelsizex |
---|
| 1042 | east = west + self.tilesize*pixelsizex |
---|
| 1043 | south = self.ominy + y*self.tilesize*pixelsizex |
---|
| 1044 | north = south + self.tilesize*pixelsizex |
---|
| 1045 | if not self.isepsg4326: |
---|
| 1046 | # Transformation to EPSG:4326 (WGS84 datum) |
---|
| 1047 | west, south = self.ct.TransformPoint(west, south)[:2] |
---|
| 1048 | east, north = self.ct.TransformPoint(east, north)[:2] |
---|
| 1049 | return south, west, north, east |
---|
| 1050 | |
---|
| 1051 | self.tileswne = rastertileswne |
---|
| 1052 | else: |
---|
| 1053 | self.tileswne = lambda x, y, z: (0,0,0,0) |
---|
| 1054 | |
---|
| 1055 | # ------------------------------------------------------------------------- |
---|
| 1056 | def generate_metadata(self): |
---|
| 1057 | """Generation of main metadata files and HTML viewers (metadata related to particular tiles are generated during the tile processing).""" |
---|
| 1058 | |
---|
| 1059 | if not os.path.exists(self.output): |
---|
| 1060 | os.makedirs(self.output) |
---|
| 1061 | |
---|
| 1062 | if self.options.profile == 'mercator': |
---|
| 1063 | |
---|
| 1064 | south, west = self.mercator.MetersToLatLon( self.ominx, self.ominy) |
---|
| 1065 | north, east = self.mercator.MetersToLatLon( self.omaxx, self.omaxy) |
---|
| 1066 | south, west = max(-85.05112878, south), max(-180.0, west) |
---|
| 1067 | north, east = min(85.05112878, north), min(180.0, east) |
---|
| 1068 | self.swne = (south, west, north, east) |
---|
| 1069 | |
---|
| 1070 | # Generate googlemaps.html |
---|
| 1071 | if self.options.webviewer in ('all','google') and self.options.profile == 'mercator': |
---|
| 1072 | if not self.options.resume or not os.path.exists(os.path.join(self.output, 'googlemaps.html')): |
---|
| 1073 | f = open(os.path.join(self.output, 'googlemaps.html'), 'w') |
---|
| 1074 | f.write( self.generate_googlemaps() ) |
---|
| 1075 | f.close() |
---|
| 1076 | |
---|
| 1077 | # Generate openlayers.html |
---|
| 1078 | if self.options.webviewer in ('all','openlayers'): |
---|
| 1079 | if not self.options.resume or not os.path.exists(os.path.join(self.output, 'openlayers.html')): |
---|
| 1080 | f = open(os.path.join(self.output, 'openlayers.html'), 'w') |
---|
| 1081 | f.write( self.generate_openlayers() ) |
---|
| 1082 | f.close() |
---|
| 1083 | |
---|
| 1084 | elif self.options.profile == 'geodetic': |
---|
| 1085 | |
---|
| 1086 | west, south = self.ominx, self.ominy |
---|
| 1087 | east, north = self.omaxx, self.omaxy |
---|
| 1088 | south, west = max(-90.0, south), max(-180.0, west) |
---|
| 1089 | north, east = min(90.0, north), min(180.0, east) |
---|
| 1090 | self.swne = (south, west, north, east) |
---|
| 1091 | |
---|
| 1092 | # Generate openlayers.html |
---|
| 1093 | if self.options.webviewer in ('all','openlayers'): |
---|
| 1094 | if not self.options.resume or not os.path.exists(os.path.join(self.output, 'openlayers.html')): |
---|
| 1095 | f = open(os.path.join(self.output, 'openlayers.html'), 'w') |
---|
| 1096 | f.write( self.generate_openlayers() ) |
---|
| 1097 | f.close() |
---|
| 1098 | |
---|
| 1099 | elif self.options.profile == 'raster': |
---|
| 1100 | |
---|
| 1101 | west, south = self.ominx, self.ominy |
---|
| 1102 | east, north = self.omaxx, self.omaxy |
---|
| 1103 | |
---|
| 1104 | self.swne = (south, west, north, east) |
---|
| 1105 | |
---|
| 1106 | # Generate openlayers.html |
---|
| 1107 | if self.options.webviewer in ('all','openlayers'): |
---|
| 1108 | if not self.options.resume or not os.path.exists(os.path.join(self.output, 'openlayers.html')): |
---|
| 1109 | f = open(os.path.join(self.output, 'openlayers.html'), 'w') |
---|
| 1110 | f.write( self.generate_openlayers() ) |
---|
| 1111 | f.close() |
---|
| 1112 | |
---|
| 1113 | |
---|
| 1114 | # Generate tilemapresource.xml. |
---|
| 1115 | if not self.options.resume or not os.path.exists(os.path.join(self.output, 'tilemapresource.xml')): |
---|
| 1116 | f = open(os.path.join(self.output, 'tilemapresource.xml'), 'w') |
---|
| 1117 | f.write( self.generate_tilemapresource()) |
---|
| 1118 | f.close() |
---|
| 1119 | |
---|
| 1120 | if self.kml: |
---|
| 1121 | # TODO: Maybe problem for not automatically generated tminz |
---|
| 1122 | # The root KML should contain links to all tiles in the tminz level |
---|
| 1123 | children = [] |
---|
| 1124 | xmin, ymin, xmax, ymax = self.tminmax[self.tminz] |
---|
| 1125 | for x in range(xmin, xmax+1): |
---|
| 1126 | for y in range(ymin, ymax+1): |
---|
| 1127 | children.append( [ x, y, self.tminz ] ) |
---|
| 1128 | # Generate Root KML |
---|
| 1129 | if self.kml: |
---|
| 1130 | if not self.options.resume or not os.path.exists(os.path.join(self.output, 'doc.kml')): |
---|
| 1131 | f = open(os.path.join(self.output, 'doc.kml'), 'w') |
---|
| 1132 | f.write( self.generate_kml( None, None, None, children) ) |
---|
| 1133 | f.close() |
---|
| 1134 | |
---|
| 1135 | # ------------------------------------------------------------------------- |
---|
| 1136 | def generate_base_tiles(self): |
---|
| 1137 | """Generation of the base tiles (the lowest in the pyramid) directly from the input raster""" |
---|
| 1138 | |
---|
| 1139 | print("Generating Base Tiles:") |
---|
| 1140 | |
---|
| 1141 | if self.options.verbose: |
---|
| 1142 | #mx, my = self.out_gt[0], self.out_gt[3] # OriginX, OriginY |
---|
| 1143 | #px, py = self.mercator.MetersToPixels( mx, my, self.tmaxz) |
---|
| 1144 | #print "Pixel coordinates:", px, py, (mx, my) |
---|
| 1145 | print('') |
---|
| 1146 | print("Tiles generated from the max zoom level:") |
---|
| 1147 | print("----------------------------------------") |
---|
| 1148 | print('') |
---|
| 1149 | |
---|
| 1150 | |
---|
| 1151 | # Set the bounds |
---|
| 1152 | tminx, tminy, tmaxx, tmaxy = self.tminmax[self.tmaxz] |
---|
| 1153 | |
---|
| 1154 | # Just the center tile |
---|
| 1155 | #tminx = tminx+ (tmaxx - tminx)/2 |
---|
| 1156 | #tminy = tminy+ (tmaxy - tminy)/2 |
---|
| 1157 | #tmaxx = tminx |
---|
| 1158 | #tmaxy = tminy |
---|
| 1159 | |
---|
| 1160 | ds = self.out_ds |
---|
| 1161 | tilebands = self.dataBandsCount + 1 |
---|
| 1162 | querysize = self.querysize |
---|
| 1163 | |
---|
| 1164 | if self.options.verbose: |
---|
| 1165 | print("dataBandsCount: ", self.dataBandsCount) |
---|
| 1166 | print("tilebands: ", tilebands) |
---|
| 1167 | |
---|
| 1168 | #print tminx, tminy, tmaxx, tmaxy |
---|
| 1169 | tcount = (1+abs(tmaxx-tminx)) * (1+abs(tmaxy-tminy)) |
---|
| 1170 | #print tcount |
---|
| 1171 | ti = 0 |
---|
| 1172 | |
---|
| 1173 | tz = self.tmaxz |
---|
| 1174 | for ty in range(tmaxy, tminy-1, -1): #range(tminy, tmaxy+1): |
---|
| 1175 | for tx in range(tminx, tmaxx+1): |
---|
| 1176 | |
---|
| 1177 | if self.stopped: |
---|
| 1178 | break |
---|
| 1179 | ti += 1 |
---|
| 1180 | tilefilename = os.path.join(self.output, str(tz), str(tx), "%s.%s" % (ty, self.tileext)) |
---|
| 1181 | if self.options.verbose: |
---|
| 1182 | print(ti,'/',tcount, tilefilename) #, "( TileMapService: z / x / y )" |
---|
| 1183 | |
---|
| 1184 | if self.options.resume and os.path.exists(tilefilename): |
---|
| 1185 | if self.options.verbose: |
---|
| 1186 | print("Tile generation skiped because of --resume") |
---|
| 1187 | else: |
---|
| 1188 | self.progressbar( ti / float(tcount) ) |
---|
| 1189 | continue |
---|
| 1190 | |
---|
| 1191 | # Create directories for the tile |
---|
| 1192 | if not os.path.exists(os.path.dirname(tilefilename)): |
---|
| 1193 | os.makedirs(os.path.dirname(tilefilename)) |
---|
| 1194 | |
---|
| 1195 | if self.options.profile == 'mercator': |
---|
| 1196 | # Tile bounds in EPSG:900913 |
---|
| 1197 | b = self.mercator.TileBounds(tx, ty, tz) |
---|
| 1198 | elif self.options.profile == 'geodetic': |
---|
| 1199 | b = self.geodetic.TileBounds(tx, ty, tz) |
---|
| 1200 | |
---|
| 1201 | #print "\tgdalwarp -ts 256 256 -te %s %s %s %s %s %s_%s_%s.tif" % ( b[0], b[1], b[2], b[3], "tiles.vrt", tz, tx, ty) |
---|
| 1202 | |
---|
| 1203 | # Don't scale up by nearest neighbour, better change the querysize |
---|
| 1204 | # to the native resolution (and return smaller query tile) for scaling |
---|
| 1205 | |
---|
| 1206 | if self.options.profile in ('mercator','geodetic'): |
---|
| 1207 | rb, wb = self.geo_query( ds, b[0], b[3], b[2], b[1]) |
---|
| 1208 | nativesize = wb[0]+wb[2] # Pixel size in the raster covering query geo extent |
---|
| 1209 | if self.options.verbose: |
---|
| 1210 | print("\tNative Extent (querysize",nativesize,"): ", rb, wb) |
---|
| 1211 | |
---|
| 1212 | # Tile bounds in raster coordinates for ReadRaster query |
---|
| 1213 | rb, wb = self.geo_query( ds, b[0], b[3], b[2], b[1], querysize=querysize) |
---|
| 1214 | |
---|
| 1215 | rx, ry, rxsize, rysize = rb |
---|
| 1216 | wx, wy, wxsize, wysize = wb |
---|
| 1217 | |
---|
| 1218 | else: # 'raster' profile: |
---|
| 1219 | |
---|
| 1220 | tsize = int(self.tsize[tz]) # tilesize in raster coordinates for actual zoom |
---|
| 1221 | xsize = self.out_ds.RasterXSize # size of the raster in pixels |
---|
| 1222 | ysize = self.out_ds.RasterYSize |
---|
| 1223 | if tz >= self.nativezoom: |
---|
| 1224 | querysize = self.tilesize # int(2**(self.nativezoom-tz) * self.tilesize) |
---|
| 1225 | |
---|
| 1226 | rx = (tx) * tsize |
---|
| 1227 | rxsize = 0 |
---|
| 1228 | if tx == tmaxx: |
---|
| 1229 | rxsize = xsize % tsize |
---|
| 1230 | if rxsize == 0: |
---|
| 1231 | rxsize = tsize |
---|
| 1232 | |
---|
| 1233 | rysize = 0 |
---|
| 1234 | if ty == tmaxy: |
---|
| 1235 | rysize = ysize % tsize |
---|
| 1236 | if rysize == 0: |
---|
| 1237 | rysize = tsize |
---|
| 1238 | ry = ysize - (ty * tsize) - rysize |
---|
| 1239 | |
---|
| 1240 | wx, wy = 0, 0 |
---|
| 1241 | wxsize, wysize = int(rxsize/float(tsize) * self.tilesize), int(rysize/float(tsize) * self.tilesize) |
---|
| 1242 | if wysize != self.tilesize: |
---|
| 1243 | wy = self.tilesize - wysize |
---|
| 1244 | |
---|
| 1245 | if self.options.verbose: |
---|
| 1246 | print("\tReadRaster Extent: ", (rx, ry, rxsize, rysize), (wx, wy, wxsize, wysize)) |
---|
| 1247 | |
---|
| 1248 | # Query is in 'nearest neighbour' but can be bigger in then the tilesize |
---|
| 1249 | # We scale down the query to the tilesize by supplied algorithm. |
---|
| 1250 | |
---|
| 1251 | # Tile dataset in memory |
---|
| 1252 | dstile = self.mem_drv.Create('', self.tilesize, self.tilesize, tilebands) |
---|
| 1253 | data = ds.ReadRaster(rx, ry, rxsize, rysize, wxsize, wysize, band_list=list(range(1,self.dataBandsCount+1))) |
---|
| 1254 | alpha = self.alphaband.ReadRaster(rx, ry, rxsize, rysize, wxsize, wysize) |
---|
| 1255 | |
---|
| 1256 | if self.tilesize == querysize: |
---|
| 1257 | # Use the ReadRaster result directly in tiles ('nearest neighbour' query) |
---|
| 1258 | dstile.WriteRaster(wx, wy, wxsize, wysize, data, band_list=list(range(1,self.dataBandsCount+1))) |
---|
| 1259 | dstile.WriteRaster(wx, wy, wxsize, wysize, alpha, band_list=[tilebands]) |
---|
| 1260 | |
---|
| 1261 | # Note: For source drivers based on WaveLet compression (JPEG2000, ECW, MrSID) |
---|
| 1262 | # the ReadRaster function returns high-quality raster (not ugly nearest neighbour) |
---|
| 1263 | # TODO: Use directly 'near' for WaveLet files |
---|
| 1264 | else: |
---|
| 1265 | # Big ReadRaster query in memory scaled to the tilesize - all but 'near' algo |
---|
| 1266 | dsquery = self.mem_drv.Create('', querysize, querysize, tilebands) |
---|
| 1267 | # TODO: fill the null value in case a tile without alpha is produced (now only png tiles are supported) |
---|
| 1268 | #for i in range(1, tilebands+1): |
---|
| 1269 | # dsquery.GetRasterBand(1).Fill(tilenodata) |
---|
| 1270 | dsquery.WriteRaster(wx, wy, wxsize, wysize, data, band_list=list(range(1,self.dataBandsCount+1))) |
---|
| 1271 | dsquery.WriteRaster(wx, wy, wxsize, wysize, alpha, band_list=[tilebands]) |
---|
| 1272 | |
---|
| 1273 | self.scale_query_to_tile(dsquery, dstile, tilefilename) |
---|
| 1274 | del dsquery |
---|
| 1275 | |
---|
| 1276 | del data |
---|
| 1277 | |
---|
| 1278 | if self.options.resampling != 'antialias': |
---|
| 1279 | # Write a copy of tile to png/jpg |
---|
| 1280 | self.out_drv.CreateCopy(tilefilename, dstile, strict=0) |
---|
| 1281 | |
---|
| 1282 | del dstile |
---|
| 1283 | |
---|
| 1284 | # Create a KML file for this tile. |
---|
| 1285 | if self.kml: |
---|
| 1286 | kmlfilename = os.path.join(self.output, str(tz), str(tx), '%d.kml' % ty) |
---|
| 1287 | if not self.options.resume or not os.path.exists(kmlfilename): |
---|
| 1288 | f = open( kmlfilename, 'w') |
---|
| 1289 | f.write( self.generate_kml( tx, ty, tz )) |
---|
| 1290 | f.close() |
---|
| 1291 | |
---|
| 1292 | if not self.options.verbose: |
---|
| 1293 | self.progressbar( ti / float(tcount) ) |
---|
| 1294 | |
---|
| 1295 | # ------------------------------------------------------------------------- |
---|
| 1296 | def generate_overview_tiles(self): |
---|
| 1297 | """Generation of the overview tiles (higher in the pyramid) based on existing tiles""" |
---|
| 1298 | |
---|
| 1299 | print("Generating Overview Tiles:") |
---|
| 1300 | |
---|
| 1301 | tilebands = self.dataBandsCount + 1 |
---|
| 1302 | |
---|
| 1303 | # Usage of existing tiles: from 4 underlying tiles generate one as overview. |
---|
| 1304 | |
---|
| 1305 | tcount = 0 |
---|
| 1306 | for tz in range(self.tmaxz-1, self.tminz-1, -1): |
---|
| 1307 | tminx, tminy, tmaxx, tmaxy = self.tminmax[tz] |
---|
| 1308 | tcount += (1+abs(tmaxx-tminx)) * (1+abs(tmaxy-tminy)) |
---|
| 1309 | |
---|
| 1310 | ti = 0 |
---|
| 1311 | |
---|
| 1312 | # querysize = tilesize * 2 |
---|
| 1313 | |
---|
| 1314 | for tz in range(self.tmaxz-1, self.tminz-1, -1): |
---|
| 1315 | tminx, tminy, tmaxx, tmaxy = self.tminmax[tz] |
---|
| 1316 | for ty in range(tmaxy, tminy-1, -1): #range(tminy, tmaxy+1): |
---|
| 1317 | for tx in range(tminx, tmaxx+1): |
---|
| 1318 | |
---|
| 1319 | if self.stopped: |
---|
| 1320 | break |
---|
| 1321 | |
---|
| 1322 | ti += 1 |
---|
| 1323 | tilefilename = os.path.join( self.output, str(tz), str(tx), "%s.%s" % (ty, self.tileext) ) |
---|
| 1324 | |
---|
| 1325 | if self.options.verbose: |
---|
| 1326 | print(ti,'/',tcount, tilefilename) #, "( TileMapService: z / x / y )" |
---|
| 1327 | |
---|
| 1328 | if self.options.resume and os.path.exists(tilefilename): |
---|
| 1329 | if self.options.verbose: |
---|
| 1330 | print("Tile generation skiped because of --resume") |
---|
| 1331 | else: |
---|
| 1332 | self.progressbar( ti / float(tcount) ) |
---|
| 1333 | continue |
---|
| 1334 | |
---|
| 1335 | # Create directories for the tile |
---|
| 1336 | if not os.path.exists(os.path.dirname(tilefilename)): |
---|
| 1337 | os.makedirs(os.path.dirname(tilefilename)) |
---|
| 1338 | |
---|
| 1339 | dsquery = self.mem_drv.Create('', 2*self.tilesize, 2*self.tilesize, tilebands) |
---|
| 1340 | # TODO: fill the null value |
---|
| 1341 | #for i in range(1, tilebands+1): |
---|
| 1342 | # dsquery.GetRasterBand(1).Fill(tilenodata) |
---|
| 1343 | dstile = self.mem_drv.Create('', self.tilesize, self.tilesize, tilebands) |
---|
| 1344 | |
---|
| 1345 | # TODO: Implement more clever walking on the tiles with cache functionality |
---|
| 1346 | # probably walk should start with reading of four tiles from top left corner |
---|
| 1347 | # Hilbert curve... |
---|
| 1348 | |
---|
| 1349 | children = [] |
---|
| 1350 | # Read the tiles and write them to query window |
---|
| 1351 | for y in range(2*ty,2*ty+2): |
---|
| 1352 | for x in range(2*tx,2*tx+2): |
---|
| 1353 | minx, miny, maxx, maxy = self.tminmax[tz+1] |
---|
| 1354 | if x >= minx and x <= maxx and y >= miny and y <= maxy: |
---|
| 1355 | dsquerytile = gdal.Open( os.path.join( self.output, str(tz+1), str(x), "%s.%s" % (y, self.tileext)), gdal.GA_ReadOnly) |
---|
| 1356 | if (ty==0 and y==1) or (ty!=0 and (y % (2*ty)) != 0): |
---|
| 1357 | tileposy = 0 |
---|
| 1358 | else: |
---|
| 1359 | tileposy = self.tilesize |
---|
| 1360 | if tx: |
---|
| 1361 | tileposx = x % (2*tx) * self.tilesize |
---|
| 1362 | elif tx==0 and x==1: |
---|
| 1363 | tileposx = self.tilesize |
---|
| 1364 | else: |
---|
| 1365 | tileposx = 0 |
---|
| 1366 | dsquery.WriteRaster( tileposx, tileposy, self.tilesize, self.tilesize, |
---|
| 1367 | dsquerytile.ReadRaster(0,0,self.tilesize,self.tilesize), |
---|
| 1368 | band_list=list(range(1,tilebands+1))) |
---|
| 1369 | children.append( [x, y, tz+1] ) |
---|
| 1370 | |
---|
| 1371 | self.scale_query_to_tile(dsquery, dstile, tilefilename) |
---|
| 1372 | # Write a copy of tile to png/jpg |
---|
| 1373 | if self.options.resampling != 'antialias': |
---|
| 1374 | # Write a copy of tile to png/jpg |
---|
| 1375 | self.out_drv.CreateCopy(tilefilename, dstile, strict=0) |
---|
| 1376 | |
---|
| 1377 | if self.options.verbose: |
---|
| 1378 | print("\tbuild from zoom", tz+1," tiles:", (2*tx, 2*ty), (2*tx+1, 2*ty),(2*tx, 2*ty+1), (2*tx+1, 2*ty+1)) |
---|
| 1379 | |
---|
| 1380 | # Create a KML file for this tile. |
---|
| 1381 | if self.kml: |
---|
| 1382 | f = open( os.path.join(self.output, '%d/%d/%d.kml' % (tz, tx, ty)), 'w') |
---|
| 1383 | f.write( self.generate_kml( tx, ty, tz, children ) ) |
---|
| 1384 | f.close() |
---|
| 1385 | |
---|
| 1386 | if not self.options.verbose: |
---|
| 1387 | self.progressbar( ti / float(tcount) ) |
---|
| 1388 | |
---|
| 1389 | |
---|
| 1390 | # ------------------------------------------------------------------------- |
---|
| 1391 | def geo_query(self, ds, ulx, uly, lrx, lry, querysize = 0): |
---|
| 1392 | """For given dataset and query in cartographic coordinates |
---|
| 1393 | returns parameters for ReadRaster() in raster coordinates and |
---|
| 1394 | x/y shifts (for border tiles). If the querysize is not given, the |
---|
| 1395 | extent is returned in the native resolution of dataset ds.""" |
---|
| 1396 | |
---|
| 1397 | geotran = ds.GetGeoTransform() |
---|
| 1398 | rx= int((ulx - geotran[0]) / geotran[1] + 0.001) |
---|
| 1399 | ry= int((uly - geotran[3]) / geotran[5] + 0.001) |
---|
| 1400 | rxsize= int((lrx - ulx) / geotran[1] + 0.5) |
---|
| 1401 | rysize= int((lry - uly) / geotran[5] + 0.5) |
---|
| 1402 | |
---|
| 1403 | if not querysize: |
---|
| 1404 | wxsize, wysize = rxsize, rysize |
---|
| 1405 | else: |
---|
| 1406 | wxsize, wysize = querysize, querysize |
---|
| 1407 | |
---|
| 1408 | # Coordinates should not go out of the bounds of the raster |
---|
| 1409 | wx = 0 |
---|
| 1410 | if rx < 0: |
---|
| 1411 | rxshift = abs(rx) |
---|
| 1412 | wx = int( wxsize * (float(rxshift) / rxsize) ) |
---|
| 1413 | wxsize = wxsize - wx |
---|
| 1414 | rxsize = rxsize - int( rxsize * (float(rxshift) / rxsize) ) |
---|
| 1415 | rx = 0 |
---|
| 1416 | if rx+rxsize > ds.RasterXSize: |
---|
| 1417 | wxsize = int( wxsize * (float(ds.RasterXSize - rx) / rxsize) ) |
---|
| 1418 | rxsize = ds.RasterXSize - rx |
---|
| 1419 | |
---|
| 1420 | wy = 0 |
---|
| 1421 | if ry < 0: |
---|
| 1422 | ryshift = abs(ry) |
---|
| 1423 | wy = int( wysize * (float(ryshift) / rysize) ) |
---|
| 1424 | wysize = wysize - wy |
---|
| 1425 | rysize = rysize - int( rysize * (float(ryshift) / rysize) ) |
---|
| 1426 | ry = 0 |
---|
| 1427 | if ry+rysize > ds.RasterYSize: |
---|
| 1428 | wysize = int( wysize * (float(ds.RasterYSize - ry) / rysize) ) |
---|
| 1429 | rysize = ds.RasterYSize - ry |
---|
| 1430 | |
---|
| 1431 | return (rx, ry, rxsize, rysize), (wx, wy, wxsize, wysize) |
---|
| 1432 | |
---|
| 1433 | # ------------------------------------------------------------------------- |
---|
| 1434 | def scale_query_to_tile(self, dsquery, dstile, tilefilename=''): |
---|
| 1435 | """Scales down query dataset to the tile dataset""" |
---|
| 1436 | |
---|
| 1437 | querysize = dsquery.RasterXSize |
---|
| 1438 | tilesize = dstile.RasterXSize |
---|
| 1439 | tilebands = dstile.RasterCount |
---|
| 1440 | |
---|
| 1441 | if self.options.resampling == 'average': |
---|
| 1442 | |
---|
| 1443 | # Function: gdal.RegenerateOverview() |
---|
| 1444 | for i in range(1,tilebands+1): |
---|
| 1445 | # Black border around NODATA |
---|
| 1446 | #if i != 4: |
---|
| 1447 | # dsquery.GetRasterBand(i).SetNoDataValue(0) |
---|
| 1448 | res = gdal.RegenerateOverview( dsquery.GetRasterBand(i), |
---|
| 1449 | dstile.GetRasterBand(i), 'average' ) |
---|
| 1450 | if res != 0: |
---|
| 1451 | self.error("RegenerateOverview() failed on %s, error %d" % (tilefilename, res)) |
---|
| 1452 | |
---|
| 1453 | elif self.options.resampling == 'antialias': |
---|
| 1454 | |
---|
| 1455 | # Scaling by PIL (Python Imaging Library) - improved Lanczos |
---|
| 1456 | array = numpy.zeros((querysize, querysize, tilebands), numpy.uint8) |
---|
| 1457 | for i in range(tilebands): |
---|
| 1458 | array[:,:,i] = gdalarray.BandReadAsArray(dsquery.GetRasterBand(i+1), 0, 0, querysize, querysize) |
---|
| 1459 | im = Image.fromarray(array, 'RGBA') # Always four bands |
---|
| 1460 | im1 = im.resize((tilesize,tilesize), Image.ANTIALIAS) |
---|
| 1461 | if os.path.exists(tilefilename): |
---|
| 1462 | im0 = Image.open(tilefilename) |
---|
| 1463 | im1 = Image.composite(im1, im0, im1) |
---|
| 1464 | im1.save(tilefilename,self.tiledriver) |
---|
| 1465 | |
---|
| 1466 | else: |
---|
| 1467 | |
---|
| 1468 | # Other algorithms are implemented by gdal.ReprojectImage(). |
---|
| 1469 | dsquery.SetGeoTransform( (0.0, tilesize / float(querysize), 0.0, 0.0, 0.0, tilesize / float(querysize)) ) |
---|
| 1470 | dstile.SetGeoTransform( (0.0, 1.0, 0.0, 0.0, 0.0, 1.0) ) |
---|
| 1471 | |
---|
| 1472 | res = gdal.ReprojectImage(dsquery, dstile, None, None, self.resampling) |
---|
| 1473 | if res != 0: |
---|
| 1474 | self.error("ReprojectImage() failed on %s, error %d" % (tilefilename, res)) |
---|
| 1475 | |
---|
| 1476 | # ------------------------------------------------------------------------- |
---|
| 1477 | def generate_tilemapresource(self): |
---|
| 1478 | """ |
---|
| 1479 | Template for tilemapresource.xml. Returns filled string. Expected variables: |
---|
| 1480 | title, north, south, east, west, isepsg4326, projection, publishurl, |
---|
| 1481 | zoompixels, tilesize, tileformat, profile |
---|
| 1482 | """ |
---|
| 1483 | |
---|
| 1484 | args = {} |
---|
| 1485 | args['title'] = self.options.title |
---|
| 1486 | args['south'], args['west'], args['north'], args['east'] = self.swne |
---|
| 1487 | args['tilesize'] = self.tilesize |
---|
| 1488 | args['tileformat'] = self.tileext |
---|
| 1489 | args['publishurl'] = self.options.url |
---|
| 1490 | args['profile'] = self.options.profile |
---|
| 1491 | |
---|
| 1492 | if self.options.profile == 'mercator': |
---|
| 1493 | args['srs'] = "EPSG:900913" |
---|
| 1494 | elif self.options.profile == 'geodetic': |
---|
| 1495 | args['srs'] = "EPSG:4326" |
---|
| 1496 | elif self.options.s_srs: |
---|
| 1497 | args['srs'] = self.options.s_srs |
---|
| 1498 | elif self.out_srs: |
---|
| 1499 | args['srs'] = self.out_srs.ExportToWkt() |
---|
| 1500 | else: |
---|
| 1501 | args['srs'] = "" |
---|
| 1502 | |
---|
| 1503 | s = """<?xml version="1.0" encoding="utf-8"?> |
---|
| 1504 | <TileMap version="1.0.0" tilemapservice="http://tms.osgeo.org/1.0.0"> |
---|
| 1505 | <Title>%(title)s</Title> |
---|
| 1506 | <Abstract></Abstract> |
---|
| 1507 | <SRS>%(srs)s</SRS> |
---|
| 1508 | <BoundingBox minx="%(south).14f" miny="%(west).14f" maxx="%(north).14f" maxy="%(east).14f"/> |
---|
| 1509 | <Origin x="%(south).14f" y="%(west).14f"/> |
---|
| 1510 | <TileFormat width="%(tilesize)d" height="%(tilesize)d" mime-type="image/%(tileformat)s" extension="%(tileformat)s"/> |
---|
| 1511 | <TileSets profile="%(profile)s"> |
---|
| 1512 | """ % args |
---|
| 1513 | for z in range(self.tminz, self.tmaxz+1): |
---|
| 1514 | if self.options.profile == 'raster': |
---|
| 1515 | s += """ <TileSet href="%s%d" units-per-pixel="%.14f" order="%d"/>\n""" % (args['publishurl'], z, (2**(self.nativezoom-z) * self.out_gt[1]), z) |
---|
| 1516 | elif self.options.profile == 'mercator': |
---|
| 1517 | s += """ <TileSet href="%s%d" units-per-pixel="%.14f" order="%d"/>\n""" % (args['publishurl'], z, 156543.0339/2**z, z) |
---|
| 1518 | elif self.options.profile == 'geodetic': |
---|
| 1519 | s += """ <TileSet href="%s%d" units-per-pixel="%.14f" order="%d"/>\n""" % (args['publishurl'], z, 0.703125/2**z, z) |
---|
| 1520 | s += """ </TileSets> |
---|
| 1521 | </TileMap> |
---|
| 1522 | """ |
---|
| 1523 | return s |
---|
| 1524 | |
---|
| 1525 | # ------------------------------------------------------------------------- |
---|
| 1526 | def generate_kml(self, tx, ty, tz, children = [], **args ): |
---|
| 1527 | """ |
---|
| 1528 | Template for the KML. Returns filled string. |
---|
| 1529 | """ |
---|
| 1530 | args['tx'], args['ty'], args['tz'] = tx, ty, tz |
---|
| 1531 | args['tileformat'] = self.tileext |
---|
| 1532 | if 'tilesize' not in args: |
---|
| 1533 | args['tilesize'] = self.tilesize |
---|
| 1534 | |
---|
| 1535 | if 'minlodpixels' not in args: |
---|
| 1536 | args['minlodpixels'] = int( args['tilesize'] / 2 ) # / 2.56) # default 128 |
---|
| 1537 | if 'maxlodpixels' not in args: |
---|
| 1538 | args['maxlodpixels'] = int( args['tilesize'] * 8 ) # 1.7) # default 2048 (used to be -1) |
---|
| 1539 | if children == []: |
---|
| 1540 | args['maxlodpixels'] = -1 |
---|
| 1541 | |
---|
| 1542 | if tx==None: |
---|
| 1543 | tilekml = False |
---|
| 1544 | args['title'] = self.options.title |
---|
| 1545 | else: |
---|
| 1546 | tilekml = True |
---|
| 1547 | args['title'] = "%d/%d/%d.kml" % (tz, tx, ty) |
---|
| 1548 | args['south'], args['west'], args['north'], args['east'] = self.tileswne(tx, ty, tz) |
---|
| 1549 | |
---|
| 1550 | if tx == 0: |
---|
| 1551 | args['drawOrder'] = 2 * tz + 1 |
---|
| 1552 | elif tx != None: |
---|
| 1553 | args['drawOrder'] = 2 * tz |
---|
| 1554 | else: |
---|
| 1555 | args['drawOrder'] = 0 |
---|
| 1556 | |
---|
| 1557 | url = self.options.url |
---|
| 1558 | if not url: |
---|
| 1559 | if tilekml: |
---|
| 1560 | url = "../../" |
---|
| 1561 | else: |
---|
| 1562 | url = "" |
---|
| 1563 | |
---|
| 1564 | s = """<?xml version="1.0" encoding="utf-8"?> |
---|
| 1565 | <kml xmlns="http://earth.google.com/kml/2.1"> |
---|
| 1566 | <Document> |
---|
| 1567 | <Name>%(title)s</Name> |
---|
| 1568 | <Description></Description> |
---|
| 1569 | <Style> |
---|
| 1570 | <ListStyle id="hideChildren"> |
---|
| 1571 | <listItemType>checkHideChildren</listItemType> |
---|
| 1572 | </ListStyle> |
---|
| 1573 | </Style>""" % args |
---|
| 1574 | if tilekml: |
---|
| 1575 | s += """ |
---|
| 1576 | <Region> |
---|
| 1577 | <Lod> |
---|
| 1578 | <minLodPixels>%(minlodpixels)d</minLodPixels> |
---|
| 1579 | <maxLodPixels>%(maxlodpixels)d</maxLodPixels> |
---|
| 1580 | </Lod> |
---|
| 1581 | <LatLonAltBox> |
---|
| 1582 | <north>%(north).14f</north> |
---|
| 1583 | <south>%(south).14f</south> |
---|
| 1584 | <east>%(east).14f</east> |
---|
| 1585 | <west>%(west).14f</west> |
---|
| 1586 | </LatLonAltBox> |
---|
| 1587 | </Region> |
---|
| 1588 | <GroundOverlay> |
---|
| 1589 | <drawOrder>%(drawOrder)d</drawOrder> |
---|
| 1590 | <Icon> |
---|
| 1591 | <href>%(ty)d.%(tileformat)s</href> |
---|
| 1592 | </Icon> |
---|
| 1593 | <LatLonBox> |
---|
| 1594 | <north>%(north).14f</north> |
---|
| 1595 | <south>%(south).14f</south> |
---|
| 1596 | <east>%(east).14f</east> |
---|
| 1597 | <west>%(west).14f</west> |
---|
| 1598 | </LatLonBox> |
---|
| 1599 | </GroundOverlay> |
---|
| 1600 | """ % args |
---|
| 1601 | |
---|
| 1602 | for cx, cy, cz in children: |
---|
| 1603 | csouth, cwest, cnorth, ceast = self.tileswne(cx, cy, cz) |
---|
| 1604 | s += """ |
---|
| 1605 | <NetworkLink> |
---|
| 1606 | <name>%d/%d/%d.%s</name> |
---|
| 1607 | <Region> |
---|
| 1608 | <Lod> |
---|
| 1609 | <minLodPixels>%d</minLodPixels> |
---|
| 1610 | <maxLodPixels>-1</maxLodPixels> |
---|
| 1611 | </Lod> |
---|
| 1612 | <LatLonAltBox> |
---|
| 1613 | <north>%.14f</north> |
---|
| 1614 | <south>%.14f</south> |
---|
| 1615 | <east>%.14f</east> |
---|
| 1616 | <west>%.14f</west> |
---|
| 1617 | </LatLonAltBox> |
---|
| 1618 | </Region> |
---|
| 1619 | <Link> |
---|
| 1620 | <href>%s%d/%d/%d.kml</href> |
---|
| 1621 | <viewRefreshMode>onRegion</viewRefreshMode> |
---|
| 1622 | <viewFormat/> |
---|
| 1623 | </Link> |
---|
| 1624 | </NetworkLink> |
---|
| 1625 | """ % (cz, cx, cy, args['tileformat'], args['minlodpixels'], cnorth, csouth, ceast, cwest, url, cz, cx, cy) |
---|
| 1626 | |
---|
| 1627 | s += """ </Document> |
---|
| 1628 | </kml> |
---|
| 1629 | """ |
---|
| 1630 | return s |
---|
| 1631 | |
---|
| 1632 | # ------------------------------------------------------------------------- |
---|
| 1633 | def generate_googlemaps(self): |
---|
| 1634 | """ |
---|
| 1635 | Template for googlemaps.html implementing Overlay of tiles for 'mercator' profile. |
---|
| 1636 | It returns filled string. Expected variables: |
---|
| 1637 | title, googlemapskey, north, south, east, west, minzoom, maxzoom, tilesize, tileformat, publishurl |
---|
| 1638 | """ |
---|
| 1639 | args = {} |
---|
| 1640 | args['title'] = self.options.title |
---|
| 1641 | args['googlemapskey'] = self.options.googlekey |
---|
| 1642 | args['south'], args['west'], args['north'], args['east'] = self.swne |
---|
| 1643 | args['minzoom'] = self.tminz |
---|
| 1644 | args['maxzoom'] = self.tmaxz |
---|
| 1645 | args['tilesize'] = self.tilesize |
---|
| 1646 | args['tileformat'] = self.tileext |
---|
| 1647 | args['publishurl'] = self.options.url |
---|
| 1648 | args['copyright'] = self.options.copyright |
---|
| 1649 | |
---|
| 1650 | s = """<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> |
---|
| 1651 | <html xmlns="http://www.w3.org/1999/xhtml" xmlns:v="urn:schemas-microsoft-com:vml"> |
---|
| 1652 | <head> |
---|
| 1653 | <title>%(title)s</title> |
---|
| 1654 | <meta http-equiv="content-type" content="text/html; charset=utf-8"/> |
---|
| 1655 | <meta http-equiv='imagetoolbar' content='no'/> |
---|
| 1656 | <style type="text/css"> v\:* {behavior:url(#default#VML);} |
---|
| 1657 | html, body { overflow: hidden; padding: 0; height: 100%%; width: 100%%; font-family: 'Lucida Grande',Geneva,Arial,Verdana,sans-serif; } |
---|
| 1658 | body { margin: 10px; background: #fff; } |
---|
| 1659 | h1 { margin: 0; padding: 6px; border:0; font-size: 20pt; } |
---|
| 1660 | #header { height: 43px; padding: 0; background-color: #eee; border: 1px solid #888; } |
---|
| 1661 | #subheader { height: 12px; text-align: right; font-size: 10px; color: #555;} |
---|
| 1662 | #map { height: 95%%; border: 1px solid #888; } |
---|
| 1663 | </style> |
---|
| 1664 | <script src='http://maps.google.com/maps?file=api&v=2&key=%(googlemapskey)s' type='text/javascript'></script> |
---|
| 1665 | <script type="text/javascript"> |
---|
| 1666 | //<![CDATA[ |
---|
| 1667 | |
---|
| 1668 | /* |
---|
| 1669 | * Constants for given map |
---|
| 1670 | * TODO: read it from tilemapresource.xml |
---|
| 1671 | */ |
---|
| 1672 | |
---|
| 1673 | var mapBounds = new GLatLngBounds(new GLatLng(%(south)s, %(west)s), new GLatLng(%(north)s, %(east)s)); |
---|
| 1674 | var mapMinZoom = %(minzoom)s; |
---|
| 1675 | var mapMaxZoom = %(maxzoom)s; |
---|
| 1676 | |
---|
| 1677 | var opacity = 0.75; |
---|
| 1678 | var map; |
---|
| 1679 | var hybridOverlay; |
---|
| 1680 | |
---|
| 1681 | /* |
---|
| 1682 | * Create a Custom Opacity GControl |
---|
| 1683 | * http://www.maptiler.org/google-maps-overlay-opacity-control/ |
---|
| 1684 | */ |
---|
| 1685 | |
---|
| 1686 | var CTransparencyLENGTH = 58; |
---|
| 1687 | // maximum width that the knob can move (slide width minus knob width) |
---|
| 1688 | |
---|
| 1689 | function CTransparencyControl( overlay ) { |
---|
| 1690 | this.overlay = overlay; |
---|
| 1691 | this.opacity = overlay.getTileLayer().getOpacity(); |
---|
| 1692 | } |
---|
| 1693 | CTransparencyControl.prototype = new GControl(); |
---|
| 1694 | |
---|
| 1695 | // This function positions the slider to match the specified opacity |
---|
| 1696 | CTransparencyControl.prototype.setSlider = function(pos) { |
---|
| 1697 | var left = Math.round((CTransparencyLENGTH*pos)); |
---|
| 1698 | this.slide.left = left; |
---|
| 1699 | this.knob.style.left = left+"px"; |
---|
| 1700 | this.knob.style.top = "0px"; |
---|
| 1701 | } |
---|
| 1702 | |
---|
| 1703 | // This function reads the slider and sets the overlay opacity level |
---|
| 1704 | CTransparencyControl.prototype.setOpacity = function() { |
---|
| 1705 | // set the global variable |
---|
| 1706 | opacity = this.slide.left/CTransparencyLENGTH; |
---|
| 1707 | this.map.clearOverlays(); |
---|
| 1708 | this.map.addOverlay(this.overlay, { zPriority: 0 }); |
---|
| 1709 | if (this.map.getCurrentMapType() == G_HYBRID_MAP) { |
---|
| 1710 | this.map.addOverlay(hybridOverlay); |
---|
| 1711 | } |
---|
| 1712 | } |
---|
| 1713 | |
---|
| 1714 | // This gets called by the API when addControl(new CTransparencyControl()) |
---|
| 1715 | CTransparencyControl.prototype.initialize = function(map) { |
---|
| 1716 | var that=this; |
---|
| 1717 | this.map = map; |
---|
| 1718 | |
---|
| 1719 | // Is this MSIE, if so we need to use AlphaImageLoader |
---|
| 1720 | var agent = navigator.userAgent.toLowerCase(); |
---|
| 1721 | if ((agent.indexOf("msie") > -1) && (agent.indexOf("opera") < 1)){this.ie = true} else {this.ie = false} |
---|
| 1722 | |
---|
| 1723 | // create the background graphic as a <div> containing an image |
---|
| 1724 | var container = document.createElement("div"); |
---|
| 1725 | container.style.width="70px"; |
---|
| 1726 | container.style.height="21px"; |
---|
| 1727 | |
---|
| 1728 | // Handle transparent PNG files in MSIE |
---|
| 1729 | if (this.ie) { |
---|
| 1730 | var loader = "filter:progid:DXImageTransform.Microsoft.AlphaImageLoader(src='http://www.maptiler.org/img/opacity-slider.png', sizingMethod='crop');"; |
---|
| 1731 | container.innerHTML = '<div style="height:21px; width:70px; ' +loader+ '" ></div>'; |
---|
| 1732 | } else { |
---|
| 1733 | container.innerHTML = '<div style="height:21px; width:70px; background-image: url(http://www.maptiler.org/img/opacity-slider.png)" ></div>'; |
---|
| 1734 | } |
---|
| 1735 | |
---|
| 1736 | // create the knob as a GDraggableObject |
---|
| 1737 | // Handle transparent PNG files in MSIE |
---|
| 1738 | if (this.ie) { |
---|
| 1739 | var loader = "progid:DXImageTransform.Microsoft.AlphaImageLoader(src='http://www.maptiler.org/img/opacity-slider.png', sizingMethod='crop');"; |
---|
| 1740 | this.knob = document.createElement("div"); |
---|
| 1741 | this.knob.style.height="21px"; |
---|
| 1742 | this.knob.style.width="13px"; |
---|
| 1743 | this.knob.style.overflow="hidden"; |
---|
| 1744 | this.knob_img = document.createElement("div"); |
---|
| 1745 | this.knob_img.style.height="21px"; |
---|
| 1746 | this.knob_img.style.width="83px"; |
---|
| 1747 | this.knob_img.style.filter=loader; |
---|
| 1748 | this.knob_img.style.position="relative"; |
---|
| 1749 | this.knob_img.style.left="-70px"; |
---|
| 1750 | this.knob.appendChild(this.knob_img); |
---|
| 1751 | } else { |
---|
| 1752 | this.knob = document.createElement("div"); |
---|
| 1753 | this.knob.style.height="21px"; |
---|
| 1754 | this.knob.style.width="13px"; |
---|
| 1755 | this.knob.style.backgroundImage="url(http://www.maptiler.org/img/opacity-slider.png)"; |
---|
| 1756 | this.knob.style.backgroundPosition="-70px 0px"; |
---|
| 1757 | } |
---|
| 1758 | container.appendChild(this.knob); |
---|
| 1759 | this.slide=new GDraggableObject(this.knob, {container:container}); |
---|
| 1760 | this.slide.setDraggableCursor('pointer'); |
---|
| 1761 | this.slide.setDraggingCursor('pointer'); |
---|
| 1762 | this.container = container; |
---|
| 1763 | |
---|
| 1764 | // attach the control to the map |
---|
| 1765 | map.getContainer().appendChild(container); |
---|
| 1766 | |
---|
| 1767 | // init slider |
---|
| 1768 | this.setSlider(this.opacity); |
---|
| 1769 | |
---|
| 1770 | // Listen for the slider being moved and set the opacity |
---|
| 1771 | GEvent.addListener(this.slide, "dragend", function() {that.setOpacity()}); |
---|
| 1772 | //GEvent.addListener(this.container, "click", function( x, y ) { alert(x, y) }); |
---|
| 1773 | |
---|
| 1774 | return container; |
---|
| 1775 | } |
---|
| 1776 | |
---|
| 1777 | // Set the default position for the control |
---|
| 1778 | CTransparencyControl.prototype.getDefaultPosition = function() { |
---|
| 1779 | return new GControlPosition(G_ANCHOR_TOP_RIGHT, new GSize(7, 47)); |
---|
| 1780 | } |
---|
| 1781 | |
---|
| 1782 | /* |
---|
| 1783 | * Full-screen Window Resize |
---|
| 1784 | */ |
---|
| 1785 | |
---|
| 1786 | function getWindowHeight() { |
---|
| 1787 | if (self.innerHeight) return self.innerHeight; |
---|
| 1788 | if (document.documentElement && document.documentElement.clientHeight) |
---|
| 1789 | return document.documentElement.clientHeight; |
---|
| 1790 | if (document.body) return document.body.clientHeight; |
---|
| 1791 | return 0; |
---|
| 1792 | } |
---|
| 1793 | |
---|
| 1794 | function getWindowWidth() { |
---|
| 1795 | if (self.innerWidth) return self.innerWidth; |
---|
| 1796 | if (document.documentElement && document.documentElement.clientWidth) |
---|
| 1797 | return document.documentElement.clientWidth; |
---|
| 1798 | if (document.body) return document.body.clientWidth; |
---|
| 1799 | return 0; |
---|
| 1800 | } |
---|
| 1801 | |
---|
| 1802 | function resize() { |
---|
| 1803 | var map = document.getElementById("map"); |
---|
| 1804 | var header = document.getElementById("header"); |
---|
| 1805 | var subheader = document.getElementById("subheader"); |
---|
| 1806 | map.style.height = (getWindowHeight()-80) + "px"; |
---|
| 1807 | map.style.width = (getWindowWidth()-20) + "px"; |
---|
| 1808 | header.style.width = (getWindowWidth()-20) + "px"; |
---|
| 1809 | subheader.style.width = (getWindowWidth()-20) + "px"; |
---|
| 1810 | // map.checkResize(); |
---|
| 1811 | } |
---|
| 1812 | |
---|
| 1813 | |
---|
| 1814 | /* |
---|
| 1815 | * Main load function: |
---|
| 1816 | */ |
---|
| 1817 | |
---|
| 1818 | function load() { |
---|
| 1819 | |
---|
| 1820 | if (GBrowserIsCompatible()) { |
---|
| 1821 | |
---|
| 1822 | // Bug in the Google Maps: Copyright for Overlay is not correctly displayed |
---|
| 1823 | var gcr = GMapType.prototype.getCopyrights; |
---|
| 1824 | GMapType.prototype.getCopyrights = function(bounds,zoom) { |
---|
| 1825 | return ["%(copyright)s"].concat(gcr.call(this,bounds,zoom)); |
---|
| 1826 | } |
---|
| 1827 | |
---|
| 1828 | map = new GMap2( document.getElementById("map"), { backgroundColor: '#fff' } ); |
---|
| 1829 | |
---|
| 1830 | map.addMapType(G_PHYSICAL_MAP); |
---|
| 1831 | map.setMapType(G_PHYSICAL_MAP); |
---|
| 1832 | |
---|
| 1833 | map.setCenter( mapBounds.getCenter(), map.getBoundsZoomLevel( mapBounds )); |
---|
| 1834 | |
---|
| 1835 | hybridOverlay = new GTileLayerOverlay( G_HYBRID_MAP.getTileLayers()[1] ); |
---|
| 1836 | GEvent.addListener(map, "maptypechanged", function() { |
---|
| 1837 | if (map.getCurrentMapType() == G_HYBRID_MAP) { |
---|
| 1838 | map.addOverlay(hybridOverlay); |
---|
| 1839 | } else { |
---|
| 1840 | map.removeOverlay(hybridOverlay); |
---|
| 1841 | } |
---|
| 1842 | } ); |
---|
| 1843 | |
---|
| 1844 | var tilelayer = new GTileLayer(GCopyrightCollection(''), mapMinZoom, mapMaxZoom); |
---|
| 1845 | var mercator = new GMercatorProjection(mapMaxZoom+1); |
---|
| 1846 | tilelayer.getTileUrl = function(tile,zoom) { |
---|
| 1847 | if ((zoom < mapMinZoom) || (zoom > mapMaxZoom)) { |
---|
| 1848 | return "http://www.maptiler.org/img/none.png"; |
---|
| 1849 | } |
---|
| 1850 | var ymax = 1 << zoom; |
---|
| 1851 | var y = ymax - tile.y -1; |
---|
| 1852 | var tileBounds = new GLatLngBounds( |
---|
| 1853 | mercator.fromPixelToLatLng( new GPoint( (tile.x)*256, (tile.y+1)*256 ) , zoom ), |
---|
| 1854 | mercator.fromPixelToLatLng( new GPoint( (tile.x+1)*256, (tile.y)*256 ) , zoom ) |
---|
| 1855 | ); |
---|
| 1856 | if (mapBounds.intersects(tileBounds)) { |
---|
| 1857 | return zoom+"/"+tile.x+"/"+y+".png"; |
---|
| 1858 | } else { |
---|
| 1859 | return "http://www.maptiler.org/img/none.png"; |
---|
| 1860 | } |
---|
| 1861 | } |
---|
| 1862 | // IE 7-: support for PNG alpha channel |
---|
| 1863 | // Unfortunately, the opacity for whole overlay is then not changeable, either or... |
---|
| 1864 | tilelayer.isPng = function() { return true;}; |
---|
| 1865 | tilelayer.getOpacity = function() { return opacity; } |
---|
| 1866 | |
---|
| 1867 | overlay = new GTileLayerOverlay( tilelayer ); |
---|
| 1868 | map.addOverlay(overlay); |
---|
| 1869 | |
---|
| 1870 | map.addControl(new GLargeMapControl()); |
---|
| 1871 | map.addControl(new GHierarchicalMapTypeControl()); |
---|
| 1872 | map.addControl(new CTransparencyControl( overlay )); |
---|
| 1873 | """ % args |
---|
| 1874 | if self.kml: |
---|
| 1875 | s += """ |
---|
| 1876 | map.addMapType(G_SATELLITE_3D_MAP); |
---|
| 1877 | map.getEarthInstance(getEarthInstanceCB); |
---|
| 1878 | """ |
---|
| 1879 | s += """ |
---|
| 1880 | |
---|
| 1881 | map.enableContinuousZoom(); |
---|
| 1882 | map.enableScrollWheelZoom(); |
---|
| 1883 | |
---|
| 1884 | map.setMapType(G_HYBRID_MAP); |
---|
| 1885 | } |
---|
| 1886 | resize(); |
---|
| 1887 | } |
---|
| 1888 | """ |
---|
| 1889 | if self.kml: |
---|
| 1890 | s += """ |
---|
| 1891 | function getEarthInstanceCB(object) { |
---|
| 1892 | var ge = object; |
---|
| 1893 | |
---|
| 1894 | if (ge) { |
---|
| 1895 | var url = document.location.toString(); |
---|
| 1896 | url = url.substr(0,url.lastIndexOf('/'))+'/doc.kml'; |
---|
| 1897 | var link = ge.createLink(""); |
---|
| 1898 | if ("%(publishurl)s") { link.setHref("%(publishurl)s/doc.kml") } |
---|
| 1899 | else { link.setHref(url) }; |
---|
| 1900 | var networkLink = ge.createNetworkLink(""); |
---|
| 1901 | networkLink.setName("TMS Map Overlay"); |
---|
| 1902 | networkLink.setFlyToView(true); |
---|
| 1903 | networkLink.setLink(link); |
---|
| 1904 | ge.getFeatures().appendChild(networkLink); |
---|
| 1905 | } else { |
---|
| 1906 | // alert("You should open a KML in Google Earth"); |
---|
| 1907 | // add div with the link to generated KML... - maybe JavaScript redirect to the URL of KML? |
---|
| 1908 | } |
---|
| 1909 | } |
---|
| 1910 | """ % args |
---|
| 1911 | s += """ |
---|
| 1912 | onresize=function(){ resize(); }; |
---|
| 1913 | |
---|
| 1914 | //]]> |
---|
| 1915 | </script> |
---|
| 1916 | </head> |
---|
| 1917 | <body onload="load()"> |
---|
| 1918 | <div id="header"><h1>%(title)s</h1></div> |
---|
| 1919 | <div id="subheader">Generated by <a href="http://www.maptiler.org/">MapTiler</a>/<a href="http://www.klokan.cz/projects/gdal2tiles/">GDAL2Tiles</a>, Copyright © 2008 <a href="http://www.klokan.cz/">Klokan Petr Pridal</a>, <a href="http://www.gdal.org/">GDAL</a> & <a href="http://www.osgeo.org/">OSGeo</a> <a href="http://code.google.com/soc/">GSoC</a> |
---|
| 1920 | <!-- PLEASE, LET THIS NOTE ABOUT AUTHOR AND PROJECT SOMEWHERE ON YOUR WEBSITE, OR AT LEAST IN THE COMMENT IN HTML. THANK YOU --> |
---|
| 1921 | </div> |
---|
| 1922 | <div id="map"></div> |
---|
| 1923 | </body> |
---|
| 1924 | </html> |
---|
| 1925 | """ % args |
---|
| 1926 | |
---|
| 1927 | return s |
---|
| 1928 | |
---|
| 1929 | |
---|
| 1930 | # ------------------------------------------------------------------------- |
---|
| 1931 | def generate_openlayers( self ): |
---|
| 1932 | """ |
---|
| 1933 | Template for openlayers.html implementing overlay of available Spherical Mercator layers. |
---|
| 1934 | |
---|
| 1935 | It returns filled string. Expected variables: |
---|
| 1936 | title, googlemapskey, yahooappid, north, south, east, west, minzoom, maxzoom, tilesize, tileformat, publishurl |
---|
| 1937 | """ |
---|
| 1938 | |
---|
| 1939 | args = {} |
---|
| 1940 | args['title'] = self.options.title |
---|
| 1941 | args['googlemapskey'] = self.options.googlekey |
---|
| 1942 | args['yahooappid'] = self.options.yahookey |
---|
| 1943 | args['south'], args['west'], args['north'], args['east'] = self.swne |
---|
| 1944 | args['minzoom'] = self.tminz |
---|
| 1945 | args['maxzoom'] = self.tmaxz |
---|
| 1946 | args['tilesize'] = self.tilesize |
---|
| 1947 | args['tileformat'] = self.tileext |
---|
| 1948 | args['publishurl'] = self.options.url |
---|
| 1949 | args['copyright'] = self.options.copyright |
---|
| 1950 | if self.options.profile == 'raster': |
---|
| 1951 | args['rasterzoomlevels'] = self.tmaxz+1 |
---|
| 1952 | args['rastermaxresolution'] = 2**(self.nativezoom) * self.out_gt[1] |
---|
| 1953 | |
---|
| 1954 | s = """<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> |
---|
| 1955 | <html xmlns="http://www.w3.org/1999/xhtml" |
---|
| 1956 | <head> |
---|
| 1957 | <title>%(title)s</title> |
---|
| 1958 | <meta http-equiv='imagetoolbar' content='no'/> |
---|
| 1959 | <style type="text/css"> v\:* {behavior:url(#default#VML);} |
---|
| 1960 | html, body { overflow: hidden; padding: 0; height: 100%%; width: 100%%; font-family: 'Lucida Grande',Geneva,Arial,Verdana,sans-serif; } |
---|
| 1961 | body { margin: 10px; background: #fff; } |
---|
| 1962 | h1 { margin: 0; padding: 6px; border:0; font-size: 20pt; } |
---|
| 1963 | #header { height: 43px; padding: 0; background-color: #eee; border: 1px solid #888; } |
---|
| 1964 | #subheader { height: 12px; text-align: right; font-size: 10px; color: #555;} |
---|
| 1965 | #map { height: 95%%; border: 1px solid #888; } |
---|
| 1966 | </style>""" % args |
---|
| 1967 | |
---|
| 1968 | if self.options.profile == 'mercator': |
---|
| 1969 | s += """ |
---|
| 1970 | <script src='http://dev.virtualearth.net/mapcontrol/mapcontrol.ashx?v=6.1'></script> |
---|
| 1971 | <script src='http://maps.google.com/maps?file=api&v=2&key=%(googlemapskey)s' type='text/javascript'></script> |
---|
| 1972 | <script src="http://api.maps.yahoo.com/ajaxymap?v=3.0&appid=%(yahooappid)s"></script>""" % args |
---|
| 1973 | |
---|
| 1974 | s += """ |
---|
| 1975 | <script src="http://www.openlayers.org/api/2.7/OpenLayers.js" type="text/javascript"></script> |
---|
| 1976 | <script type="text/javascript"> |
---|
| 1977 | var map; |
---|
| 1978 | var mapBounds = new OpenLayers.Bounds( %(west)s, %(south)s, %(east)s, %(north)s); |
---|
| 1979 | var mapMinZoom = %(minzoom)s; |
---|
| 1980 | var mapMaxZoom = %(maxzoom)s; |
---|
| 1981 | |
---|
| 1982 | // avoid pink tiles |
---|
| 1983 | OpenLayers.IMAGE_RELOAD_ATTEMPTS = 3; |
---|
| 1984 | OpenLayers.Util.onImageLoadErrorColor = "transparent"; |
---|
| 1985 | |
---|
| 1986 | function init(){""" % args |
---|
| 1987 | |
---|
| 1988 | if self.options.profile == 'mercator': |
---|
| 1989 | s += """ |
---|
| 1990 | var options = { |
---|
| 1991 | controls: [], |
---|
| 1992 | projection: new OpenLayers.Projection("EPSG:900913"), |
---|
| 1993 | displayProjection: new OpenLayers.Projection("EPSG:4326"), |
---|
| 1994 | units: "m", |
---|
| 1995 | maxResolution: 156543.0339, |
---|
| 1996 | maxExtent: new OpenLayers.Bounds(-20037508, -20037508, 20037508, 20037508.34) |
---|
| 1997 | }; |
---|
| 1998 | map = new OpenLayers.Map('map', options); |
---|
| 1999 | |
---|
| 2000 | // create Google Mercator layers |
---|
| 2001 | var gmap = new OpenLayers.Layer.Google("Google Streets", |
---|
| 2002 | { sphericalMercator: true, numZoomLevels: 20} ); |
---|
| 2003 | var gsat = new OpenLayers.Layer.Google("Google Satellite", |
---|
| 2004 | {type: G_SATELLITE_MAP, sphericalMercator: true, numZoomLevels: 20} ); |
---|
| 2005 | var ghyb = new OpenLayers.Layer.Google("Google Hybrid", |
---|
| 2006 | {type: G_HYBRID_MAP, sphericalMercator: true, numZoomLevels: 20}); |
---|
| 2007 | var gter = new OpenLayers.Layer.Google("Google Terrain", |
---|
| 2008 | {type: G_PHYSICAL_MAP, sphericalMercator: true, numZoomLevels: 20 }); |
---|
| 2009 | |
---|
| 2010 | // create Virtual Earth layers |
---|
| 2011 | OpenLayers.Layer.VirtualEarth.prototype.MAX_ZOOM_LEVEL=19; |
---|
| 2012 | OpenLayers.Layer.VirtualEarth.prototype.RESOLUTIONS=OpenLayers.Layer.Google.prototype.RESOLUTIONS |
---|
| 2013 | var veroad = new OpenLayers.Layer.VirtualEarth("Virtual Earth Roads", |
---|
| 2014 | {'type': VEMapStyle.Road, 'sphericalMercator': true, numZoomLevels: 20}); |
---|
| 2015 | var veaer = new OpenLayers.Layer.VirtualEarth("Virtual Earth Aerial", |
---|
| 2016 | {'type': VEMapStyle.Aerial, 'sphericalMercator': true, numZoomLevels: 20 }); |
---|
| 2017 | var vehyb = new OpenLayers.Layer.VirtualEarth("Virtual Earth Hybrid", |
---|
| 2018 | {'type': VEMapStyle.Hybrid, 'sphericalMercator': true}); |
---|
| 2019 | |
---|
| 2020 | // create Yahoo layer |
---|
| 2021 | var yahoo = new OpenLayers.Layer.Yahoo("Yahoo Street", |
---|
| 2022 | {'sphericalMercator': true}); |
---|
| 2023 | var yahoosat = new OpenLayers.Layer.Yahoo("Yahoo Satellite", |
---|
| 2024 | {'type': YAHOO_MAP_SAT, 'sphericalMercator': true}); |
---|
| 2025 | var yahoohyb = new OpenLayers.Layer.Yahoo("Yahoo Hybrid", |
---|
| 2026 | {'type': YAHOO_MAP_HYB, 'sphericalMercator': true}); |
---|
| 2027 | |
---|
| 2028 | // create OSM/OAM layer |
---|
| 2029 | var osm = new OpenLayers.Layer.TMS( "OpenStreetMap", |
---|
| 2030 | "http://tile.openstreetmap.org/", |
---|
| 2031 | { type: 'png', getURL: osm_getTileURL, displayOutsideMaxExtent: true, attribution: '<a href="http://www.openstreetmap.org/">OpenStreetMap</a>'} ); |
---|
| 2032 | var oam = new OpenLayers.Layer.TMS( "OpenAerialMap", |
---|
| 2033 | "http://tile.openaerialmap.org/tiles/1.0.0/openaerialmap-900913/", |
---|
| 2034 | { type: 'png', getURL: osm_getTileURL } ); |
---|
| 2035 | |
---|
| 2036 | // create TMS Overlay layer |
---|
| 2037 | var tmsoverlay = new OpenLayers.Layer.TMS( "TMS Overlay", "", |
---|
| 2038 | { // url: '', serviceVersion: '.', layername: '.', |
---|
| 2039 | type: 'png', getURL: overlay_getTileURL, alpha: true, |
---|
| 2040 | isBaseLayer: false |
---|
| 2041 | }); |
---|
| 2042 | if (OpenLayers.Util.alphaHack() == false) { tmsoverlay.setOpacity(0.7); } |
---|
| 2043 | |
---|
| 2044 | map.addLayers([gmap, gsat, ghyb, gter, veroad, veaer, vehyb, |
---|
| 2045 | yahoo, yahoosat, yahoohyb, osm, oam, |
---|
| 2046 | tmsoverlay]); |
---|
| 2047 | |
---|
| 2048 | var switcherControl = new OpenLayers.Control.LayerSwitcher(); |
---|
| 2049 | map.addControl(switcherControl); |
---|
| 2050 | switcherControl.maximizeControl(); |
---|
| 2051 | |
---|
| 2052 | map.zoomToExtent( mapBounds.transform(map.displayProjection, map.projection ) ); |
---|
| 2053 | """ % args |
---|
| 2054 | |
---|
| 2055 | elif self.options.profile == 'geodetic': |
---|
| 2056 | s += """ |
---|
| 2057 | var options = { |
---|
| 2058 | controls: [], |
---|
| 2059 | projection: new OpenLayers.Projection("EPSG:4326"), |
---|
| 2060 | maxResolution: 0.703125, |
---|
| 2061 | maxExtent: new OpenLayers.Bounds(-180, -90, 180, 90) |
---|
| 2062 | }; |
---|
| 2063 | map = new OpenLayers.Map('map', options); |
---|
| 2064 | |
---|
| 2065 | layer = new OpenLayers.Layer.WMS( "Blue Marble", |
---|
| 2066 | "http://labs.metacarta.com/wms-c/Basic.py?", {layers: 'satellite' } ); |
---|
| 2067 | map.addLayer(layer); |
---|
| 2068 | wms = new OpenLayers.Layer.WMS( "VMap0", |
---|
| 2069 | "http://labs.metacarta.com/wms-c/Basic.py?", {layers: 'basic', format: 'image/png' } ); |
---|
| 2070 | map.addLayer(wms); |
---|
| 2071 | |
---|
| 2072 | var tmsoverlay = new OpenLayers.Layer.TMS( "TMS Overlay", "", |
---|
| 2073 | { |
---|
| 2074 | serviceVersion: '.', layername: '.', alpha: true, |
---|
| 2075 | type: 'png', getURL: overlay_getTileURL, |
---|
| 2076 | isBaseLayer: false |
---|
| 2077 | }); |
---|
| 2078 | map.addLayer(tmsoverlay); |
---|
| 2079 | if (OpenLayers.Util.alphaHack() == false) { tmsoverlay.setOpacity(0.7); } |
---|
| 2080 | |
---|
| 2081 | var switcherControl = new OpenLayers.Control.LayerSwitcher(); |
---|
| 2082 | map.addControl(switcherControl); |
---|
| 2083 | switcherControl.maximizeControl(); |
---|
| 2084 | |
---|
| 2085 | map.zoomToExtent( mapBounds ); |
---|
| 2086 | """ |
---|
| 2087 | |
---|
| 2088 | elif self.options.profile == 'raster': |
---|
| 2089 | s += """ |
---|
| 2090 | var options = { |
---|
| 2091 | controls: [], |
---|
| 2092 | maxExtent: new OpenLayers.Bounds( %(west)s, %(south)s, %(east)s, %(north)s ), |
---|
| 2093 | maxResolution: %(rastermaxresolution)f, |
---|
| 2094 | numZoomLevels: %(rasterzoomlevels)d |
---|
| 2095 | }; |
---|
| 2096 | map = new OpenLayers.Map('map', options); |
---|
| 2097 | |
---|
| 2098 | var layer = new OpenLayers.Layer.TMS( "TMS Layer","", |
---|
| 2099 | { url: '', serviceVersion: '.', layername: '.', alpha: true, |
---|
| 2100 | type: 'png', getURL: overlay_getTileURL |
---|
| 2101 | }); |
---|
| 2102 | map.addLayer(layer); |
---|
| 2103 | map.zoomToExtent( mapBounds ); |
---|
| 2104 | """ % args |
---|
| 2105 | |
---|
| 2106 | |
---|
| 2107 | s += """ |
---|
| 2108 | map.addControl(new OpenLayers.Control.PanZoomBar()); |
---|
| 2109 | map.addControl(new OpenLayers.Control.MousePosition()); |
---|
| 2110 | map.addControl(new OpenLayers.Control.MouseDefaults()); |
---|
| 2111 | map.addControl(new OpenLayers.Control.KeyboardDefaults()); |
---|
| 2112 | } |
---|
| 2113 | """ % args |
---|
| 2114 | |
---|
| 2115 | if self.options.profile == 'mercator': |
---|
| 2116 | s += """ |
---|
| 2117 | function osm_getTileURL(bounds) { |
---|
| 2118 | var res = this.map.getResolution(); |
---|
| 2119 | var x = Math.round((bounds.left - this.maxExtent.left) / (res * this.tileSize.w)); |
---|
| 2120 | var y = Math.round((this.maxExtent.top - bounds.top) / (res * this.tileSize.h)); |
---|
| 2121 | var z = this.map.getZoom(); |
---|
| 2122 | var limit = Math.pow(2, z); |
---|
| 2123 | |
---|
| 2124 | if (y < 0 || y >= limit) { |
---|
| 2125 | return "http://www.maptiler.org/img/none.png"; |
---|
| 2126 | } else { |
---|
| 2127 | x = ((x %% limit) + limit) %% limit; |
---|
| 2128 | return this.url + z + "/" + x + "/" + y + "." + this.type; |
---|
| 2129 | } |
---|
| 2130 | } |
---|
| 2131 | |
---|
| 2132 | function overlay_getTileURL(bounds) { |
---|
| 2133 | var res = this.map.getResolution(); |
---|
| 2134 | var x = Math.round((bounds.left - this.maxExtent.left) / (res * this.tileSize.w)); |
---|
| 2135 | var y = Math.round((bounds.bottom - this.tileOrigin.lat) / (res * this.tileSize.h)); |
---|
| 2136 | var z = this.map.getZoom(); |
---|
| 2137 | if (this.map.baseLayer.name == 'Virtual Earth Roads' || this.map.baseLayer.name == 'Virtual Earth Aerial' || this.map.baseLayer.name == 'Virtual Earth Hybrid') { |
---|
| 2138 | z = z + 1; |
---|
| 2139 | } |
---|
| 2140 | if (mapBounds.intersectsBounds( bounds ) && z >= mapMinZoom && z <= mapMaxZoom ) { |
---|
| 2141 | //console.log( this.url + z + "/" + x + "/" + y + "." + this.type); |
---|
| 2142 | return this.url + z + "/" + x + "/" + y + "." + this.type; |
---|
| 2143 | } else { |
---|
| 2144 | return "http://www.maptiler.org/img/none.png"; |
---|
| 2145 | } |
---|
| 2146 | } |
---|
| 2147 | """ % args |
---|
| 2148 | |
---|
| 2149 | elif self.options.profile == 'geodetic': |
---|
| 2150 | s += """ |
---|
| 2151 | function overlay_getTileURL(bounds) { |
---|
| 2152 | bounds = this.adjustBounds(bounds); |
---|
| 2153 | var res = this.map.getResolution(); |
---|
| 2154 | var x = Math.round((bounds.left - this.tileOrigin.lon) / (res * this.tileSize.w)); |
---|
| 2155 | var y = Math.round((bounds.bottom - this.tileOrigin.lat) / (res * this.tileSize.h)); |
---|
| 2156 | var z = this.map.getZoom(); |
---|
| 2157 | var path = this.serviceVersion + "/" + this.layername + "/" + z + "/" + x + "/" + y + "." + this.type; |
---|
| 2158 | var url = this.url; |
---|
| 2159 | if (mapBounds.intersectsBounds( bounds ) && z >= mapMinZoom && z <= mapMaxZoom) { |
---|
| 2160 | // console.log( this.url + z + "/" + x + "/" + y + "." + this.type); |
---|
| 2161 | return this.url + z + "/" + x + "/" + y + "." + this.type; |
---|
| 2162 | } else { |
---|
| 2163 | return "http://www.maptiler.org/img/none.png"; |
---|
| 2164 | } |
---|
| 2165 | } |
---|
| 2166 | """ % args |
---|
| 2167 | |
---|
| 2168 | elif self.options.profile == 'raster': |
---|
| 2169 | s += """ |
---|
| 2170 | function overlay_getTileURL(bounds) { |
---|
| 2171 | var res = this.map.getResolution(); |
---|
| 2172 | var x = Math.round((bounds.left - this.maxExtent.left) / (res * this.tileSize.w)); |
---|
| 2173 | var y = Math.round((bounds.bottom - this.maxExtent.bottom) / (res * this.tileSize.h)); |
---|
| 2174 | var z = this.map.getZoom(); |
---|
| 2175 | if (x >= 0 && y >= 0) { |
---|
| 2176 | return this.url + z + "/" + x + "/" + y + "." + this.type; |
---|
| 2177 | } else { |
---|
| 2178 | return "http://www.maptiler.org/img/none.png"; |
---|
| 2179 | } |
---|
| 2180 | } |
---|
| 2181 | """ % args |
---|
| 2182 | |
---|
| 2183 | s += """ |
---|
| 2184 | function getWindowHeight() { |
---|
| 2185 | if (self.innerHeight) return self.innerHeight; |
---|
| 2186 | if (document.documentElement && document.documentElement.clientHeight) |
---|
| 2187 | return document.documentElement.clientHeight; |
---|
| 2188 | if (document.body) return document.body.clientHeight; |
---|
| 2189 | return 0; |
---|
| 2190 | } |
---|
| 2191 | |
---|
| 2192 | function getWindowWidth() { |
---|
| 2193 | if (self.innerWidth) return self.innerWidth; |
---|
| 2194 | if (document.documentElement && document.documentElement.clientWidth) |
---|
| 2195 | return document.documentElement.clientWidth; |
---|
| 2196 | if (document.body) return document.body.clientWidth; |
---|
| 2197 | return 0; |
---|
| 2198 | } |
---|
| 2199 | |
---|
| 2200 | function resize() { |
---|
| 2201 | var map = document.getElementById("map"); |
---|
| 2202 | var header = document.getElementById("header"); |
---|
| 2203 | var subheader = document.getElementById("subheader"); |
---|
| 2204 | map.style.height = (getWindowHeight()-80) + "px"; |
---|
| 2205 | map.style.width = (getWindowWidth()-20) + "px"; |
---|
| 2206 | header.style.width = (getWindowWidth()-20) + "px"; |
---|
| 2207 | subheader.style.width = (getWindowWidth()-20) + "px"; |
---|
| 2208 | if (map.updateSize) { map.updateSize(); }; |
---|
| 2209 | } |
---|
| 2210 | |
---|
| 2211 | onresize=function(){ resize(); }; |
---|
| 2212 | |
---|
| 2213 | </script> |
---|
| 2214 | </head> |
---|
| 2215 | <body onload="init()"> |
---|
| 2216 | <div id="header"><h1>%(title)s</h1></div> |
---|
| 2217 | <div id="subheader">Generated by <a href="http://www.maptiler.org/">MapTiler</a>/<a href="http://www.klokan.cz/projects/gdal2tiles/">GDAL2Tiles</a>, Copyright © 2008 <a href="http://www.klokan.cz/">Klokan Petr Pridal</a>, <a href="http://www.gdal.org/">GDAL</a> & <a href="http://www.osgeo.org/">OSGeo</a> <a href="http://code.google.com/soc/">GSoC</a> |
---|
| 2218 | <!-- PLEASE, LET THIS NOTE ABOUT AUTHOR AND PROJECT SOMEWHERE ON YOUR WEBSITE, OR AT LEAST IN THE COMMENT IN HTML. THANK YOU --> |
---|
| 2219 | </div> |
---|
| 2220 | <div id="map"></div> |
---|
| 2221 | <script type="text/javascript" >resize()</script> |
---|
| 2222 | </body> |
---|
| 2223 | </html>""" % args |
---|
| 2224 | |
---|
| 2225 | return s |
---|
| 2226 | |
---|
| 2227 | # ============================================================================= |
---|
| 2228 | # ============================================================================= |
---|
| 2229 | # ============================================================================= |
---|
| 2230 | |
---|
| 2231 | if __name__=='__main__': |
---|
| 2232 | argv = gdal.GeneralCmdLineProcessor( sys.argv ) |
---|
| 2233 | if argv: |
---|
| 2234 | gdal2tiles = GDAL2Tiles( argv[1:] ) |
---|
| 2235 | gdal2tiles.process() |
---|