Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

gdal2tiles.py @ 23

Revision 23, 86.8 KB checked in by astrid.emde, 14 years ago (diff)

Line
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( txself.tileSize, tyself.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	txself.tileSizeres - 180,
394	tyself.tileSizeres - 90,
395	(tx+1)self.tileSizeres - 180,
396	(ty+1)self.tileSizeres - 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(2tz-1, tmaxx), min(2tz-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(2tz-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] + xself.tilesizepixelsizex
1042	east = west + self.tilesize*pixelsizex
1043	south = self.ominy + yself.tilesizepixelsizex
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('', 2self.tilesize, 2self.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(2ty,2ty+2):
1352	for x in range(2tx,2tx+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 % (2tx) 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:", (2tx, 2ty), (2tx+1, 2ty),(2tx, 2ty+1), (2tx+1, 2ty+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()

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: branches/0.5/win/PostNAS-0.5/bin/gdal2tiles.py @ 23

Download in other formats: