zoom.py 6.5 KB

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