source: src/django_gheat/gheat/globalmaptiles.py@ 9029

Last change on this file since 9029 was 9028, checked in by dennisw, 14 years ago

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