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