Some overengineered code for taking a bounding box and computing its constituent geohashes

pyrolite.py

import math
import geohash as ghash

WGS84_RADIUS = 6370997.0 # meters
METERS_PER_DEGREE = WGS84_RADIUS * 2 * math.pi / 360.0

class Point(list):
    def __init__(self, latitude, longitude):
        list.__init__(self, (latitude, longitude))

    @property
    def latitude(self):
        return self[0]

    @property
    def longitude(self):
        return self[1]

    def __repr__(self):
        return "Point(%.4f %.4f)" % self

    def distance_to(self, other):
        if abs(self.latitude) == 90 and self.latitude == self.longitude:
            return 0.0

        # Haversine distance
        lat1, lon1, lat2, lon2 = [
            n * math.pi / 180 for n in self.latitude, self.longitude, other.latitude, other.longitude]
        d_lat = (lat2 - lat1) / 2.0
        d_lon = (lon2 - lon1) / 2.0
        a = math.sin(d_lat) ** 2 + math.cos(lat1) * math.cos(lat2) * math.sin(d_lon) ** 2
        c = 2 * math.atan2(math.sqrt(a), math.sqrt(1.0 - a))
        return WGS84_RADIUS * c # return meters

    @property
    def hash(self):
        return ghash.encode(self.latitude, self.longitude)


class Envelope(object):
    __slots__ = ("south", "west", "north", "east", "depth", "geohash")

    def __init__(self, southwest=(-90.0, -180.0), northeast=(90.0, 180.0), depth=""):
        self.south, self.west = southwest
        self.north, self.east = northeast
        self.depth = depth
        self.geohash = None

    def __repr__(self):
        return "Envelope(%.4f %.4f, %.4f %.4f, %s)" % (
            self.south, self.west, self.north, self.east, self.depth)

    @property
    def center(self):
        return Point((self.south + self.north) / 2,
                     (self.west + self.east) / 2)

    @property
    def hash(self):
        if self.geohash: return self.geohash
        center = self.center
        self.geohash = ghash.encode(center.latitude, center.longitude)[:int(len(self.depth)/5)]
        return self.geohash 

    @classmethod
    def from_hash(cls, hash):
        lat, lon, dlat, dlon = ghash.decode_exactly(hash)
        return cls((lat-dlat, lon-dlon), (lat+dlat, lon+dlon), "0"*(len(hash)*5))

    def contains(self, other):
        return (self.west <= other.east <= self.east and
                self.west <= other.west <= self.east and
                self.south <= other.south <= self.north and
                self.south <= other.north <= self.north)

    def overlaps(self, other):
       return (not self.east < other.west and
               not self.west > other.east and
               not self.north < other.south and
               not self.south > other.north)

    def children(self):
        mid_lat = (self.south + self.north) / 2.0
        mid_lon = (self.west + self.east) / 2.0
        return (
            Envelope((self.south, self.west), (mid_lat, mid_lon), self.depth+"00"),
            Envelope((self.south, mid_lon), (mid_lat, self.east), self.depth+"10"),
            Envelope((mid_lat, self.west), (self.north, mid_lon), self.depth+"01"),
            Envelope((mid_lat, mid_lon), (self.north, self.east), self.depth+"11")
        )

    def count_hashes(self, cells):
        return len(set(cell.depth[:int(len(cell.depth)/5)*5] for cell in cells))
       
    """This is actually where the magic happens: takes an arbitrary envelope and returns
       a set of n geohashes that minimally contain it."""
    def coverage(self, min_count=32):
        """Returns a list of uniformly sized regions that completely cover the envelope."""
        cells = [Envelope()]
        while self.count_hashes(cells) < min_count:
            new_cells = []
            for cell in cells:
                new_cells.extend([c for c in cell.children() if self.overlaps(c)])
            cells = new_cells
        return cells

    def snap_coverage(self, min_count=32):
        hashes = set()
        result = []
        for cell in self.coverage(min_count):
            if cell.hash in hashes: continue
            result.append(cell)
            hashes.add(cell.hash)
        return result

class Circle(Envelope):
    def __init__(self, origin, radius):
        assert isinstance(origin, Point)
        radius_lat = radius / METERS_PER_DEGREE
        if origin.latitude == 90:
            bottom_left = Point(90 - radius_lat, -180)
            top_right = Point(90, 180)
        elif origin.latitude == -90:
            bottom_left = Point(-90, -180)
            top_right = Point(-90 + radius_lat, 180)
        else:
            radius_lon = radius_lat / math.cos(math.radians(origin.latitude))
            bottom_left = Point(origin.latitude - radius_lat,
                                origin.longitude - radius_lon)
            top_right   = Point(origin.latitude + radius_lat,
                                origin.longitude + radius_lon)
        Envelope.__init__(self, bottom_left, top_right)
        self.origin = origin
        self.radius = radius

    @property
    def center(self):
        return self.origin

    def __repr__(self):
        return "Circle(%s %.4f)" % (self.origin, self.radius)

    def coverage(self, min_count=32):
        cells = Envelope.coverage(self, min_count)
        cells = [(cell.center.distance_to(self.center), cell) for cell in cells]
        cells.sort()
        return [cell for distance, cell in cells]

    sorted_coverage = Envelope.snap_coverage

if __name__ == "__main__":
    import sys
    lat1, lon1, lat2, lon2  = map(float, sys.argv[1:5])
    circle = Circle(Point(lat1, lon1), 1000.0)
    print [cell.hash for cell in circle.coverage()]