Skip to content

Instantly share code, notes, and snippets.

@JoppeSchwartz
Last active September 8, 2022 18:05
Show Gist options
  • Save JoppeSchwartz/61db22a48de29414db80 to your computer and use it in GitHub Desktop.
Save JoppeSchwartz/61db22a48de29414db80 to your computer and use it in GitHub Desktop.
Given a set of latitude-longitude coordinates representing a route, this module will create a set of bounding boxes.
# 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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment