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