Given a set of latitude-longitude coordinates representing a route, this module will create a set of bounding boxes.

boxer.py

#	module boxer.py
#
#	Given a set of latitude-longitude coordinates representing a route, this module will
#	create a set of bounding boxes, accessible via the boxes property.
#	

import math
from pprint import pprint

class coordinate:
	"""Class representing a latitude/longitude coordinate."""
	def __init__(self, lat, lng):
		self.lat = float(lat)
		self.lng = float(lng)
	
	def __str__(self):
		return str(self.__dict__)
	
	def __eq__(self, other):
		return ((self.lat == other.lat) and (self.lng == other.lng))

	
class point:
	"""Class representing a single point (in mercator space)."""
	def __init__(self, x_coord, y_coord):
		self.x = float(x_coord)
		self.y = float(y_coord)
		
	def __str__(self):
		return str(self.__dict__)
	
	def __eq__(self, other):
		return (self.x == other.x) and (self.y == other.y)

class size:
	"""Class representing a size in terms of width and height."""
	def __init__(self, wi, ht):
		self.width = float(math.fabs(wi))
		self.height = float(math.fabs(ht))

	def __str__(self):
		return str(self.__dict__)

class rect:
	"""Class representing a rectangle in terms of an origin point and size."""
	def __init__(self, pt, sz):
		self.origin = point(pt.x, pt.y)
		self.size = size(sz.width, sz.height)
		
		ulx = self.origin.x - self.size.width/2.0
		uly = self.origin.y + self.size.height/2.0
		self.upperLeft = point(ulx, uly)
		
		lrx = self.origin.x + self.size.width/2.0
		lry = self.origin.y - self.size.height/2.0
		self.lowerRight = point(lrx, lry)

	def __str__(self):
		return "{{origin: {0}  size: {1}}}".format(str(self.origin), str(self.size))


	#	Method to return the union of two rectangles.
	def union(self, rect2):
		minX = min(self.upperLeft.x, rect2.upperLeft.x)
		maxX = max(self.lowerRight.x, rect2.lowerRight.x)
		minY = min(self.lowerRight.y, rect2.lowerRight.y)
		maxY = max(self.upperLeft.y, rect2.upperLeft.y)
		return rect(point((minX + maxX) / 2.0, (minY + maxY) / 2.0), 
			size(math.fabs(maxX-minX), math.fabs(maxY-minY)))

class boundingBox:
	"""Utility class to represent a bounding box in terms of upper left
		and lower right coordinates."""
	def __init__(self, ul, lr):
		self.upperLeft = ul
		self.lowerRight = lr
	
	def __str__(self):
 		return "{{ul: {0} lr: {1}}}".format(str(self.upperLeft), str(self.lowerRight))
 		
 		
class grid:
	"""Class representing a grid of cell objects; its methods help build up a 
	list of marked cells.
	Internally, grid uses a set to keep track of marked cells.
	Each marked cell is represented as a tuple of the form (x-index, y-index)."""
	def __init__(self, boundingRect, edgeLength):
		self.boundingRect = boundingRect
		self.edgeLength = edgeLength
		self.upperLeftCell = self.cellForMapPoint(boundingRect.upperLeft)
		self.lowerRightCell = self.cellForMapPoint(boundingRect.lowerRight)
		self.xIndices = range(self.upperLeftCell[0], self.lowerRightCell[0]+1) 
		self.yIndices = range(self.lowerRightCell[1], self.upperLeftCell[1]+1)
		self.markedCells = set()
		
	def __str__(self):
		return str(self.__dict__)

	def markCell(self, cell ):
		self.markedCells.add( cell )
	
	def unmarkCell(self, cell ):
		self.markedCells.remove( cell )
		
	def markCellAndNeighbors(self, cell ):
		minX = int(max(min(self.xIndices), cell[0]-1))
		maxX = int(min(max(self.xIndices), cell[0]+2))
		minY = int(max(min(self.yIndices), cell[1]-1))
		maxY = int(min(max(self.yIndices), cell[1]+2))
		for xi in range(minX, maxX):
			for yi in range(minY, maxY):
				self.markCell( (xi, yi) )
	
	def cellMarked(self, cell):
		return cell in self.markedCells
		
	def cellForMapPoint(self, pt):
		"""Returns a tuple of the form (x-index, y-index) corresponding to the cell
		for the given map (Mercator) point"""
		normX = pt.x - self.boundingRect.origin.x
		normY = pt.y - self.boundingRect.origin.y
		xIdx = int(round( normX / self.edgeLength ))
		yIdx = int(round( normY / self.edgeLength ))
		return (xIdx, yIdx)
	
	def mapPointForCell(self, cell ):
		x = cell[0] * self.edgeLength + self.boundingRect.origin.x
		y = cell[1] * self.edgeLength + self.boundingRect.origin.y
		return point(x, y)
		
	def mapRectForCell(self, cell ):
		return rect( self.mapPointForCell( cell ), 
			size(self.edgeLength, self.edgeLength) )
		
	def cellsAreAdjacent(self, cell1, cell2):
		return (((abs(cell1[0] - cell2[0]) == 1) and (cell1[1] == cell2[1])) 
			or ((cell1[0] == cell2[0]) and (abs(cell1[1] - cell2[1]) == 1)))	
	
# Mercator translation functions
EARTH_RADIUS = 6378137
def deg2rad(deg):
	return deg * math.pi / 180.0
def rad2deg(rad):
	return rad * 180.0 / math.pi
	
def x2lng(x):
	return rad2deg( x/EARTH_RADIUS )
	
def lng2x(lng):
	return deg2rad(lng) * EARTH_RADIUS

def y2lat(y):
	return rad2deg( (2.0*math.atan(math.exp(y / EARTH_RADIUS))-math.pi/2.0) )

def lat2y(lat):
	return math.log( math.tan ( math.pi / 4.0 + deg2rad(lat) / 2.0 ) ) * EARTH_RADIUS
#  return 180.0/math.pi*math.log(math.tan(math.pi/4.0+a*(math.pi/180.0)/2.0))


def coord2MercPoint(coord):
	return point(lng2x(coord.lng), lat2y(coord.lat))

def mercPoint2Coord(pt):
	return coordinate( y2lat(pt.y), x2lng(pt.x))

class RouteBoxer:
	"""Class to decompose a route of points, given in lat/lng format, into a list 
	of bounding boxes."""
	EDGE_LENGTH = 2000	# Edge length in meters
	
	def __init__(self, path, upperLeft, lowerRight):
		self.path = path
		self.upperLeftCoord = upperLeft
		self.lowerRightCoord = lowerRight
		self.upperLeftMapPoint = coord2MercPoint(upperLeft)
		self.lowerRightMapPoint = coord2MercPoint(lowerRight)
		deltaX = math.fabs(self.lowerRightMapPoint.x - self.upperLeftMapPoint.x)
		deltaY = math.fabs(self.upperLeftMapPoint.y - self.lowerRightMapPoint.y)
		originX = self.upperLeftMapPoint.x + deltaX / 2.0
		originY = self.lowerRightMapPoint.y + deltaY / 2.0
		boundingBox = rect(point(originX, originY), size(deltaX, deltaY))
		self.grid = grid(boundingBox, RouteBoxer.EDGE_LENGTH)
		self.boxesX = []
		self.boxesY = []
		self.buildGrid()
		self.mergeIntersectingCells()
	
	def buildGrid(self):
		#self.grid = grid(boundingBox, RouteBoxer.EDGE_LENGTH)
		lastMapPoint = coord2MercPoint(self.path[0])
		lastCell = self.grid.cellForMapPoint(lastMapPoint)
		self.grid.markCellAndNeighbors(lastCell)
		
		for loc in self.path:
			curMapPoint = coord2MercPoint(loc)
			curCell = self.grid.cellForMapPoint(curMapPoint)
			# TODO:  test for failure to find cell
			# If the current cell is the same cell as the last one, skip it.
			if curCell == lastCell:
				continue
			self.grid.markCellAndNeighbors(curCell)
			
			# If the cell is next to the last cell, continue.
			if self.grid.cellsAreAdjacent(curCell, lastCell):
				continue
			# Otherwise, the cells are at some distance from each other.
			# If the cells share the same X or Y index, mark the intervening ones.
			if curCell[0] == lastCell[0]:
				lowerY = min(lastCell[1], curCell[1])
				upperY = max(lastCell[1], curCell[1])
				for y in range (lowerY, upperY + 1):
					self.grid.markCellAndNeighbors( (curCell[0], y) )
			
			elif curCell[1] == lastCell[1]:
				lowerX = min(lastCell[0], curCell[0])
				upperX = max(lastCell[0], curCell[0])
				for x in range(lowerX, upperX+1):
					self.grid.markCellAndNeighbors( (x, curCell[1]) )
			else:
			# The cells lie on a line not sharing an X or Y index.
			# Calculate the slope and length of the line connecting the last and
			# current map points.
				deltaX = curMapPoint.x - lastMapPoint.x
				deltaY = curMapPoint.y - lastMapPoint.y
				length = math.sqrt(pow(deltaY, 2) + pow(deltaX, 2))
				theta = math.atan2(deltaY, deltaX)
			# Calculate the number of mapPointRange segments are in the line
			# connecting the last and current map point.  Iterate over them and 
			# mark the enclosing cells. This is equivalent to walking the connecting
			# line in segments of mapPointRange length and marking the enclosing
			# cells. Since the segment used is shorter than the sides of the cells,
			# we're guaranteed to hit almost every intervening cell. The
			# markCellAndNeighbors function takes care of the rest.
				numSegments = int(math.floor(length / RouteBoxer.EDGE_LENGTH))
				if numSegments > 1:
					for i in range(numSegments):
						nextX = lastMapPoint.x + i * RouteBoxer.EDGE_LENGTH * math.cos(theta)
						nextY = lastMapPoint.y + i * RouteBoxer.EDGE_LENGTH * math.sin(theta)
						nextMapPoint = point(nextX, nextY)
						nextCell = self.grid.cellForMapPoint(nextMapPoint)
						self.grid.markCellAndNeighbors(nextCell)
			
			lastMapPoint = curMapPoint
			lastCell = curCell			
	
	def mergeIntersectingCells(self):
		self.boxesX = []
		self.boxesY = []
		curMapRect= None
		for y in self.grid.yIndices:
			for x in self.grid.xIndices:
				if self.grid.cellMarked( (x,y) ):
					if curMapRect is None:
						curMapRect = self.grid.mapRectForCell( (x,y) )
					else:
						newRect = self.grid.mapRectForCell( (x,y) )
						curMapRect = curMapRect.union(newRect)	
				else:
					if curMapRect:
						self.boxesX.append(curMapRect)
						curMapRect = None
			if curMapRect:
				self.boxesX.append(curMapRect)
				curMapRect = None
		
		curMapRect = None
		for x in self.grid.xIndices:
			for y in self.grid.yIndices:
				if self.grid.cellMarked( (x,y) ):
					if curMapRect is None:
						curMapRect = self.grid.mapRectForCell( (x,y) )
					else:
						curMapRect = curMapRect.union(self.grid.mapRectForCell( (x,y) ))
				else:
					if curMapRect:
						self.boxesY.append(curMapRect)
						curMapRect = None
			if curMapRect:
				self.boxesY.append(curMapRect)
				curMapRect = None
		
	def boxes(self):
		if len(self.boxesX) < len(self.boxesY):
			return self.boxesX
		else:
			return self.boxesY
	
	def boxCoords(self):
		coords = []
		for box in self.boxes():
			ul = mercPoint2Coord(box.upperLeft)
			lr = mercPoint2Coord(box.lowerRight)
			coords.append(boundingBox(ul, lr))

		return coords	


# The remaining code is for testing.		

shapePoints =  ["40.730426",
                "-73.98634300000001",
                "40.729778",
                "-73.986808",
                "40.729778",
                "-73.986808",
                "40.728782",
                "-73.984436",
                "40.728782",
                "-73.984436",
                "40.747833",
                "-73.97049699999999",
                "40.747833",
                "-73.97049699999999",
                "40.751506",
                "-73.967803",
                "40.751506",
                "-73.967803",
                "40.752311",
                "-73.96700199999999",
                "40.752311",
                "-73.96700199999999",
                "40.753406",
                "-73.96333300000001",
                "40.753406",
                "-73.96333300000001",
                "40.774616",
                "-73.94306899999999",
                "40.782993",
                "-73.94375599999999",
                "40.797519",
                "-73.929176",
                "40.797519",
                "-73.929176",
                "40.799316",
                "-73.92913",
                "40.799316",
                "-73.92913",
                "40.8031",
                "-73.929748",
                "40.8031",
                "-73.929748",
                "40.804843",
                "-73.926208",
                "40.804843",
                "-73.926208",
                "40.805072",
                "-73.923301",
                "40.805072",
                "-73.923301",
                "40.804153",
                "-73.912933",
                "40.804153",
                "-73.912933",
                "40.804519",
                "-73.91154400000001",
                "40.804519",
                "-73.91154400000001",
                "40.822376",
                "-73.887001",
                "40.829086",
                "-73.835975",
                "40.829086",
                "-73.835975",
                "40.835872",
                "-73.824676",
                "40.859527",
                "-73.82719400000001",
                "40.859527",
                "-73.82719400000001",
                "40.860515",
                "-73.827316",
                "40.860515",
                "-73.827316",
                "40.863605",
                "-73.825912",
                "40.863605",
                "-73.825912",
                "40.871265",
                "-73.818969",
                "40.883472",
                "-73.816078",
                "40.883472",
                "-73.816078",
                "40.886619",
                "-73.813957",
                "40.887519",
                "-73.802474",
                "40.906391",
                "-73.79098500000001",
                "40.941104",
                "-73.75151",
                "40.959793",
                "-73.74033300000001",
                "40.965164",
                "-73.72895800000001",
                "40.97406",
                "-73.72264800000001",
                "40.975032",
                "-73.702934",
                "40.986881",
                "-73.681915",
                "40.988426",
                "-73.666122",
                "41.017787",
                "-73.63499400000001",
                "41.027153",
                "-73.60393500000001",
                "41.035533",
                "-73.59448999999999",
                "41.04454",
                "-73.568077",
                "41.044246",
                "-73.553276",
                "41.04747",
                "-73.54207599999999",
                "41.04747",
                "-73.54207599999999",
                "41.048072",
                "-73.539299",
                "41.048072",
                "-73.539299",
                "41.053428",
                "-73.53936"]



def doit():
	#	Bounds of testing rectangle.
	lr = coordinate(40.728782, -73.539299)
	ul = coordinate(41.053428, -73.986808)
	
	# Build list of coordinates from the shape points
	shapeCoords = []
	n = len(shapePoints)
	print "Analyzing", n, "points"
	for idx in range(len(shapePoints)/2):
		lat = shapePoints[idx*2]
		lng = shapePoints[idx*2 + 1]
		coord = coordinate(lat, lng)
		shapeCoords.append(coord)

	print shapeCoords[0]
	rb = RouteBoxer(shapeCoords, ul, lr)
	#import pdb; pdb.set_trace()
	rb.buildGrid()
	print "Marked Cells:", pprint(rb.grid.markedCells)
	rb.mergeIntersectingCells()
	print "\nBoxesX:", ', '.join(map(str, rb.boxesX))
	print "\nBoxesY:", ', '.join(map(str, rb.boxesY))
	boxes = rb.boxes()
	if rb.boxesX == boxes:
		print "Least boxes is boxesX with {0} boxes".format(len(rb.boxesX))
	else:
		print "Least boxes is boxesY with {0} boxes".format(len(rb.boxesY))
	
	print "\nBox coords:", ', '.join(map(str, rb.boxCoords()))
	

if __name__ == "__main__":
	doit()