1 from django.contrib.gis.geos import GEOSGeometry, LinearRing, Polygon, Point

     2 from django.contrib.gis.maps.google.gmap import GoogleMapException

     3 from math import pi, sin, cos, log, exp, atan

     4 

     5 # Constants used for degree to radian conversion, and vice-versa.

     6 DTOR = pi / 180.

     7 RTOD = 180. / pi

     8 

     9 class 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