The here's the next challenge: take a geographic area of interest, and find urls for all the tiles which covering that area (at zoom levels of interest).
My area of interest will be the fine City of New Orleans
'''
The City of New Orleans
xmin: -90.283741
ymin: 29.890626
xmax: -89.912952
ymax: 30.057766
'''
Disclaimer: Persisting map tiles locally breaks the terms of service for many popular basemap providers like Google, Bing, and ArcGIS Online. I'm not a lawyer, nor do I aspire to be one, but if you start storing tiles offline without permission, then you do so at your own risk.
Helper Methods for Finding Tiles
There are a bunch of functions you will need at your disposal to accomplish this task. These functions can be found in Joe Schwartz's Awesome White Paper on the Bing Maps Tiling System. This is one of the best explanations of map tile systems I've found on the web. Joe goes over the mercator projection, along with explanations of ground resolution, map scale, and pixel coordinates. Really really great stuff and I urge everybody to really dig into it.
At the end of the article, Joe provides sample code which contains the major functions you will need to dominate map tile systems. While I've grown to really like C# development, I almost always translate interesting code back into warm and friendly Python. Here's Joe's TileSystem class translated into Python with a couple extra methods for additional abstraction:
'''
Mostly direct port of awesome article by Joe Schwartz - http://msdn.microsoft.com/en-us/library/bb259689.aspx
'''
from __future__ import division
import math
class TileUtils(object):
earth_radius = 6378137
min_lat = -85.05112878
max_lat = 85.05112878
min_lng = -180
max_lng = 180
def clipValue(self, value, minValue, maxValue):
'''
Makes sure that value is within a specific range.
If not, then the lower or upper bounds is returned
'''
return min(max(value, minValue), maxValue)
def getMapDimensionsByZoomLevel(self, zoomLevel):
'''
Returns the width/height in pixels of the entire map
based on the zoom level.
'''
return 256 << zoomLevel
def getGroundResolution(self, latitude, level):
'''
returns the ground resolution for based on latitude and zoom level.
'''
latitude = self.clipValue(latitude, self.min_lat, self.max_lat);
mapSize = self.getMapDimensionsByZoomLevel(level)
return math.cos(latitude * math.pi / 180) * 2 * math.pi * self.earth_radius / mapSize
def getMapScale(self, latitude, level, dpi=96):
'''
returns the map scale on the dpi of the screen
'''
dpm = dpi / 0.0254 #convert to dots per meter
return self.getGroundResolution(latitude, level) * dpm
def convertLatLngToPixelXY(self, lat, lng, level):
'''
returns the x and y values of the pixel corresponding to a latitude and longitude.
'''
mapSize = self.getMapDimensionsByZoomLevel(level)
lat = self.clipValue(lat, self.min_lat, self.max_lat)
lng = self.clipValue(lng, self.min_lng, self.max_lng)
x = (lng + 180) / 360
sinlat = math.sin(lat * math.pi / 180)
y = 0.5 - math.log((1 + sinlat) / (1 - sinlat)) / (4 * math.pi)
pixelX = int(self.clipValue(x * mapSize + 0.5, 0, mapSize - 1))
pixelY = int(self.clipValue(y * mapSize + 0.5, 0, mapSize - 1))
return (pixelX, pixelY)
def convertPixelXYToLngLat(self, pixelX, pixelY, level):
'''
converts a pixel x, y to a latitude and longitude.
'''
mapSize = self.getMapDimensionsByZoomLevel(level)
x = (self.clipValue(pixelX, 0, mapSize - 1) / mapSize) - 0.5
y = 0.5 - (self.clipValue(pixelY, 0, mapSize - 1) / mapSize)
lat = 90 - 360 * math.atan(math.exp(-y * 2 * math.pi)) / math.pi
lng = 360 * x
return (lng, lat)
def convertPixelXYToTileXY(self, pixelX, pixelY):
'''
Converts pixel XY coordinates into tile XY coordinates of the tile containing
'''
return(int(pixelX / 256), int(pixelY / 256))
def convertTileXYToPixelXY(self, tileX, tileY):
'''
Converts tile XY coordinates into pixel XY coordinates of the upper-left pixel
'''
return(tileX * 256, tileY * 256)
def tileXYZToQuadKey(self, x, y, z):
'''
Computes quadKey value based on tile x, y and z values.
'''
quadKey = ''
for i in range(z, 0, -1):
digit = 0
mask = 1 << (i - 1)
if(x & mask) != 0:
digit += 1
if(y & mask) != 0:
digit += 2
quadKey += str(digit)
return quadKey
def quadKeyToTileXYZ(self, quadKey):
'''
Computes tile x, y and z values based on quadKey.
'''
tileX = 0
tileY = 0
tileZ = len(quadKey)
for i in range(tileZ, 0, -1):
mask = 1 << (i - 1)
value = quadKey[tileZ - i]
if value == '0':
continue
elif value == '1':
tileX |= mask
elif value == '2':
tileY |= mask
elif value == '3':
tileX |= mask
tileY |= mask
else:
raise Exception('Invalid QuadKey')
return (tileX, tileY, tileZ)
def convertLngLatToTileXY(self, lng, lat, level):
pixelX, pixelY = self.convertLatLngToPixelXY(lat, lng, level)
return self.convertPixelXYToTileXY(pixelX, pixelY)
def getTileOrigin(self, tileX, tileY, level):
'''
Returns the upper-left hand corner lat/lng for a tile
'''
pixelX, pixelY = self.convertTileXYToPixelXY(tileX, tileY)
lng, lat = self.convertPixelXYToLngLat(pixelX, pixelY, level)
return (lat, lng)
The function convertLngLatToTileXY really sums up the class for our purposes. This function will help get the tile x/y values for the sides of the bounding box. Using these functions, we can write an additional helper class to generate all the tile urls within our bounding box based off of a template tile url like the ones mentioned in my last post:
class TileFinder(object):
def __init__(self, tileUtils, tileTemplate):
'''
Tile template is a template url which use {{x}}, {{y}}, {{z}} which will be replaced by
their respective values.
'''
self.tileUtils = tileUtils
self.tileTemplate = tileTemplate
def getTileUrlsByLatLngExtent(self, xmin, ymin, xmax, ymax, level):
'''
Returns a list of tile urls by extent
'''
#Upper-Left Tile
tileXMin, tileYMin = self.tileUtils.convertLngLatToTileXY(xmin, ymax, level)
#Lower-Right Tile
tileXMax, tileYMax = self.tileUtils.convertLngLatToTileXY(xmax, ymin, level)
tileUrls = []
for y in range(tileYMax, tileYMin - 1, -1):
for x in range(tileXMin, tileXMax + 1, 1):
tileUrls.append(self.createTileUrl(x, y, level))
return tileUrls
def createTileUrl(self, x, y, z):
'''
returns new tile url based on template
'''
return self.tileTemplate.replace('{{x}}', str(x)).replace('{{y}}', str(y)).replace('{{z}}', str(z))
Check out the constructor here. It has two required arguments, one of which is a tile template. This tile template uses {{x}}, {{y}}, and {{z}} values to indicate where substitutions should take place. These substitutions happen in createTileUrl which is called for each x,y pair generated in getTileUrlsByLatLngExtent.
So now we have our helper classes for generating a list of tile urls based on an extent. The last value to complete the trifecta is the zoom level. I chose level 11 since my area of interest is fairly small (city level). Here is a snippet for running the code:
'''
The City of New Orleans
xmin: -90.283741
ymin: 29.890626
xmax: -89.912952
ymax: 30.057766
'''
import urllib2
#constructing the template url for arcgis online
template = 'http://services.arcgisonline.com/ArcGIS/rest/services/World_Topo_Map/MapServer/tile/{{z}}/{{y}}/{{x}}.png'
tileUtils = TileUtils()
tileFinder = TileFinder(tileUtils, template)
tile_urls = tileFinder.getTileUrlsByLatLngExtent(-90.283741, 29.890626, -89.912952, 30.057766, 11)
for t in tile_urls: urllib2.urlopen(t).read()
So that's about all there is to it. Below you can see the six tile urls that were produced in the code snippet above. A special thanks to Stephen Ansari (@sansari22) who really got me started on learning tile schemas and provided me sample code to get going.
 |  |  |
 |  |  |
