Ushahidi Clustering - Python

gistfile1.py

#!/usr/bin/env python
import math
import sys
import MySQLdb

MAP_ZOOM     = 21
MAP_OFFSET   = 268435456    # half the Earth's circumference at zoom level 21
MAP_RADIUS   = MAP_OFFSET / math.pi

def longitude_to_x(longitude):
    """
    Returns the x pixel coordinate for a given longitude value.
    Longitude values are in decimal degrees format.
    """
    x = round(MAP_OFFSET + MAP_RADIUS * longitude * math.pi / 180)

    return x

def latitude_to_y(latitude):
    """
    Returns the y pixel coordinate for a given latitude value.
    Latitude values are in decimal degrees format.
    """
    y = round(MAP_OFFSET - MAP_RADIUS * \
              math.log((1 + math.sin(latitude * math.pi / 180)) / \
                       (1 - math.sin(latitude * math.pi / 180)) \
                      ) / 2 \
             )

    return y

def plane_distance(latitude1, longitude1, latitude2, longitude2, zoom):
    """
    Given a pair of lat/long coordinates and a map zoom level, returns
    the distance between the two points in pixels
    """
    x1 = longitude_to_x(longitude1)
    y1 = latitude_to_y(latitude1)

    x2 = longitude_to_x(longitude2)
    y2 = latitude_to_y(latitude2)

    distance = int(math.sqrt((x1 - x2)**2) + math.sqrt((y1 - y2)**2))

    return int(distance) >> (int(MAP_ZOOM) - int(zoom))

def cluster(points, cluster_distance, zoom):
    """
    Groups points that are less than cluster_distance pixels apart at
    a given zoom level into a cluster.
    """
    distance = (10000000 >> int(zoom)) / 100000
    
    clusters = []

    while len(points) > 0:
        point1 = points.pop()

        cluster = []

        for point2 in points[:]:

            pixel_distance = plane_distance(point1['latitude'],
                                            point1['longitude'],
                                            point2['latitude'],
                                            point2['longitude'],
                                            zoom)
                                            
#            pixel_distance = abs(point1['longitude']-point2['longitude']) + abs(point1['latitude']-point2['latitude'])
            if pixel_distance < distance:
                points.remove(point2)
                cluster.append(point2)

        # add the first point to the cluster
        if len(cluster) > 0:
            cluster.append(point1)
            clusters.append(cluster)
        else:
            clusters.append([point1])

    return clusters
    
# Zoom
try:
    sys.argv[1]
    zoom = sys.argv[1]
except NameError:
    zoom = 8
    
# Category
try:
    sys.argv[2]
    category_id = sys.argv[2]
except NameError:
    category_id = 0    

# Start Date
try:
    sys.argv[3]
    start_date = sys.argv[3]
except NameError:
    start_date = 0

# Start Date
try:
    sys.argv[4]
    end_date = sys.argv[4]
except NameError:
    end_date = 0
    
# Create the filter
db_filter = ""
if (category_id != 0):
    db_filter += " AND ( c.id="+category_id+" OR c.parent_id="+category_id+") "
if (start_date != 0):
    db_filter += " AND i.incident_date >= '"+start_date+"'"
if (end_date != 0):
    db_filter += " AND i.incident_date <= '"+end_date+"'"


# Connect To Ushahidi Database
db = MySQLdb.connect(host="localhost", user="root", passwd="voodoo", db="ushahidi")

# create a database cursor
cursor = db.cursor()

# execute SQL select statement
cursor.execute("SELECT DISTINCT i.id, i.incident_title, l.`latitude`, l.`longitude` FROM `incident` AS i INNER JOIN `location` AS l ON (l.`id` = i.`location_id`) INNER JOIN `incident_category` AS ic ON (i.`id` = ic.`incident_id`) INNER JOIN `category` AS c ON (ic.`category_id` = c.`id`) WHERE i.incident_active=1 "+db_filter+" ORDER BY i.`id` ASC")
rows = cursor.fetchall()
db.close()

points = [{'latitude': float(row[2]),  
    'longitude': float(row[3])} for row in rows]
    
point1 = points.pop()
point2 = points.pop()

clusters = cluster(points, 20, zoom)

print clusters

You can get division by 0 when latitude_to_y(90.0)