/django/contrib/gis/maps/google/zoom.py

  1from django.contrib.gis.geos import GEOSGeometry, LinearRing, Polygon, Point
  2from django.contrib.gis.maps.google.gmap import GoogleMapException
  3from math import pi, sin, cos, log, exp, atan
  4
  5# Constants used for degree to radian conversion, and vice-versa.
  6DTOR = pi / 180.
  7RTOD = 180. / pi
  8
  9class GoogleZoom(object):
 10    """
 11    GoogleZoom is a utility for performing operations related to the zoom
 12    levels on Google Maps.
 13
 14    This class is inspired by the OpenStreetMap Mapnik tile generation routine
 15    `generate_tiles.py`, and the article "How Big Is the World" (Hack #16) in
 16    "Google Maps Hacks" by Rich Gibson and Schuyler Erle.
 17
 18    `generate_tiles.py` may be found at:
 19      http://trac.openstreetmap.org/browser/applications/rendering/mapnik/generate_tiles.py
 20
 21    "Google Maps Hacks" may be found at http://safari.oreilly.com/0596101619
 22    """
 23    
 24    def __init__(self, num_zoom=19, tilesize=256):
 25        "Initializes the Google Zoom object."
 26        # Google's tilesize is 256x256, square tiles are assumed.
 27        self._tilesize = tilesize
 28        
 29        # The number of zoom levels
 30        self._nzoom = num_zoom
 31
 32        # Initializing arrays to hold the parameters for each one of the 
 33        # zoom levels.
 34        self._degpp = [] # Degrees per pixel
 35        self._radpp = [] # Radians per pixel
 36        self._npix  = [] # 1/2 the number of pixels for a tile at the given zoom level
 37        
 38        # Incrementing through the zoom levels and populating the parameter arrays.
 39        z = tilesize # The number of pixels per zoom level.
 40        for i in xrange(num_zoom):
 41            # Getting the degrees and radians per pixel, and the 1/2 the number of
 42            # for every zoom level.
 43            self._degpp.append(z / 360.) # degrees per pixel
 44            self._radpp.append(z / (2 * pi)) # radians per pixl
 45            self._npix.append(z / 2) # number of pixels to center of tile
 46
 47            # Multiplying `z` by 2 for the next iteration.
 48            z *= 2
 49
 50    def __len__(self):
 51        "Returns the number of zoom levels."
 52        return self._nzoom
 53
 54    def get_lon_lat(self, lonlat):
 55        "Unpacks longitude, latitude from GEOS Points and 2-tuples."
 56        if isinstance(lonlat, Point):
 57            lon, lat = lonlat.coords
 58        else:
 59            lon, lat = lonlat
 60        return lon, lat
 61
 62    def lonlat_to_pixel(self, lonlat, zoom):
 63        "Converts a longitude, latitude coordinate pair for the given zoom level."
 64        # Setting up, unpacking the longitude, latitude values and getting the
 65        # number of pixels for the given zoom level.
 66        lon, lat = self.get_lon_lat(lonlat)
 67        npix = self._npix[zoom]
 68
 69        # Calculating the pixel x coordinate by multiplying the longitude value
 70        # with with the number of degrees/pixel at the given zoom level.
 71        px_x = round(npix + (lon * self._degpp[zoom]))
 72
 73        # Creating the factor, and ensuring that 1 or -1 is not passed in as the 
 74        # base to the logarithm.  Here's why:
 75        #  if fac = -1, we'll get log(0) which is undefined; 
 76        #  if fac =  1, our logarithm base will be divided by 0, also undefined.
 77        fac = min(max(sin(DTOR * lat), -0.9999), 0.9999)
 78
 79        # Calculating the pixel y coordinate.
 80        px_y = round(npix + (0.5 * log((1 + fac)/(1 - fac)) * (-1.0 * self._radpp[zoom])))
 81
 82        # Returning the pixel x, y to the caller of the function.
 83        return (px_x, px_y)
 84
 85    def pixel_to_lonlat(self, px, zoom):
 86        "Converts a pixel to a longitude, latitude pair at the given zoom level."
 87        if len(px) != 2:
 88            raise TypeError('Pixel should be a sequence of two elements.')
 89
 90        # Getting the number of pixels for the given zoom level.
 91        npix = self._npix[zoom]
 92
 93        # Calculating the longitude value, using the degrees per pixel.
 94        lon = (px[0] - npix) / self._degpp[zoom]
 95
 96        # Calculating the latitude value.
 97        lat = RTOD * ( 2 * atan(exp((px[1] - npix)/ (-1.0 * self._radpp[zoom]))) - 0.5 * pi)
 98
 99        # Returning the longitude, latitude coordinate pair.
100        return (lon, lat)
101    
102    def tile(self, lonlat, zoom):
103        """
104        Returns a Polygon  corresponding to the region represented by a fictional
105        Google Tile for the given longitude/latitude pair and zoom level. This
106        tile is used to determine the size of a tile at the given point.
107        """
108        # The given lonlat is the center of the tile.
109        delta = self._tilesize / 2
110
111        # Getting the pixel coordinates corresponding to the
112        # the longitude/latitude.
113        px = self.lonlat_to_pixel(lonlat, zoom)
114
115        # Getting the lower-left and upper-right lat/lon coordinates
116        # for the bounding box of the tile.
117        ll = self.pixel_to_lonlat((px[0]-delta, px[1]-delta), zoom)
118        ur = self.pixel_to_lonlat((px[0]+delta, px[1]+delta), zoom)
119
120        # Constructing the Polygon, representing the tile and returning.
121        return Polygon(LinearRing(ll, (ll[0], ur[1]), ur, (ur[0], ll[1]), ll), srid=4326)
122        
123    def get_zoom(self, geom):
124        "Returns the optimal Zoom level for the given geometry."
125        # Checking the input type.
126        if not isinstance(geom, GEOSGeometry) or geom.srid != 4326:
127            raise TypeError('get_zoom() expects a GEOS Geometry with an SRID of 4326.')
128
129        # Getting the envelope for the geometry, and its associated width, height
130        # and centroid.
131        env = geom.envelope
132        env_w, env_h = self.get_width_height(env.extent)
133        center = env.centroid
134
135        for z in xrange(self._nzoom):
136            # Getting the tile at the zoom level.
137            tile_w, tile_h = self.get_width_height(self.tile(center, z).extent)
138
139            # When we span more than one tile, this is an approximately good
140            # zoom level.
141            if (env_w > tile_w) or (env_h > tile_h):
142                if z == 0: 
143                    raise GoogleMapException('Geometry width and height should not exceed that of the Earth.')
144                return z-1
145        
146        # Otherwise, we've zoomed in to the max.
147        return self._nzoom-1
148
149    def get_width_height(self, extent):
150        """
151        Returns the width and height for the given extent.
152        """
153        # Getting the lower-left, upper-left, and upper-right
154        # coordinates from the extent.
155        ll = Point(extent[:2])
156        ul = Point(extent[0], extent[3])
157        ur = Point(extent[2:])
158        # Calculating the width and height.
159        height = ll.distance(ul)
160        width  = ul.distance(ur)
161        return width, height