Identify clusters of oft-repeated fire pixels with a date range to create non-fire mask (polygon)

fireloc_cluster.py

#!/usr/bin/env python

import os, sys
from optparse import OptionParser
from datetime import datetime
import logging
import psycopg2, psycopg2.extras
import math
from collections import defaultdict
#
# adapted from http://stackoverflow.com/questions/10108368/detecting-geographic-clusters
#
desc_text = """Identify clusters of oft-repeated fire pixels with a date range to create non-fire mask (polygon)"""

usage_text = """usage: %prog [options] start_date end_date
   start/end dates in yyyy-mm-dd format"""

class ActiveFire:
    def __init__(self, dsn):
        self.dsn = dsn
        self.conn = psycopg2.connect(dsn)
        #self.cur = self.conn.cursor(cursor_factory=psycopg2.extras.DictCursor)
        self.cur = self.conn.cursor()
    def __del__(self):
        if hasattr(self, 'cur'):
            self.cur.close()
        if hasattr(self, 'conn'):
            self.conn.close()
    def execute(self, sql, paramdict=None, commit=None):
        self.cur.execute(sql, paramdict)
        if commit is True:
            self.conn.commit()
        return self.cur
    def commit(self):
        self.conn.commit()
        

to_rad = math.pi / 180.0   # convert lat or lng to radians
threshhold_dist_m=800      # adjust to your needs
threshhold_locations=None  # minimum # of locations needed in a cluster

def dist(lat1,lng1,lat2,lng2):
    '''Compute the distance in metres between 2 coordinates'''
    global to_rad
    earth_radius_m = 6371000

    dLat = (lat2-lat1) * to_rad
    dLon = (lng2-lng1) * to_rad
    lat1_rad = lat1 * to_rad
    lat2_rad = lat2 * to_rad

    a = math.sin(dLat/2) * math.sin(dLat/2) + math.sin(dLon/2) * math.sin(dLon/2) * math.cos(lat1_rad) * math.cos(lat2_rad)
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a)); 
    dist = earth_radius_m * c
    return dist

def sitesDist(site1,site2): 
    #just a helper to shorted list comprehension below 
    return dist(site1[0],site1[1], site2[0], site2[1])

def main(argv=None):
    if argv is None:
        argv = sys.argv

    debuglevelD = {
        'debug': logging.DEBUG,
        'info': logging.INFO,
        'warning': logging.WARNING,
        'error': logging.ERROR,
        'critical': logging.CRITICAL,
    }
    
    defvals = {
        'dateformat': '%Y-%m-%d',
        'start_hour': 8,
        'end_hour': 16,
        'sql0': '''SELECT ST_X(the_geom), ST_Y(the_geom) FROM af_modis_2012 WHERE obs_time_stamp >= %(start_date)s AND obs_time_stamp < %(end_date)s AND date_part('hour', obs_time_stamp) > %(start_hour)s AND date_part('hour', obs_time_stamp) < %(end_hour)s''',
        'sql1': '''INSERT INTO af_mask(start_date, end_date, the_geom) VALUES(%(start_date)s, %(end_date)s, ST_ConvexHull(ST_GeomFromText(%(wkt)s, 4326)))''',
    }
    
    parser = OptionParser(usage=usage_text, description=desc_text)
    parser.add_option("--dsn", dest="dsn", type="string", \
        help="Database connection string", metavar="'host=HOST port=PORT dbname=DATABASE user=DBUSER password=PASSWORD'")
    parser.add_option("--select", dest="sql0", type="string", \
        help="Statement that returns a list of (longitude,latitude) tuples (%s)"%defvals['sql0'], metavar="SQL")
    parser.add_option("--start-hour", dest="start_hour", type="int", \
        help="Filter for fire after start hour (%s)"%defvals['start_hour'], metavar="INT")
    parser.add_option("--end-hour", dest="end_hour", type="int", \
        help="Filter for fire before end hour (%s)"%defvals['end_hour'], metavar="INT")        
    parser.add_option("--insert", dest="sql1", type="string", \
        help="Statement that inserts MULTIPOINT WKT into database (%s)"%defvals['sql1'], metavar="SQL")
    parser.add_option("-l", "--loglevel", dest="loglevel", type="string", \
        help="Verbosity %s"%debuglevelD.keys(), metavar='LOGLEVEL')
    parser.set_defaults(**defvals)
    (options, args) = parser.parse_args()

    if options.loglevel:
        if options.loglevel not in debuglevelD: raise AssertionError("Log level must be one of: %s"%debuglevelD.keys())
        dbglvl = debuglevelD[options.loglevel]
    else:
        dbglvl = logging.WARNING
    logger = logging.getLogger()
    logger.setLevel(dbglvl)
    ch = logging.StreamHandler()
    ch.setFormatter( logging.Formatter('%(asctime)s %(lineno)d %(name)s %(funcName)s - %(levelname)s - %(message)s') )
    ch.setLevel(dbglvl)
    logger.addHandler(ch)

    if len(args) != 2:
        parser.error("Requires start and end dates in %s format" % options.dateformat)
    
    try:
        sdate = datetime.strptime(args[0], options.dateformat)
        edate = datetime.strptime(args[1], options.dateformat)
    except Exception, err:
        parser.error("Invalid input date: %s" % err)

    numdays = ((edate - sdate).days - 1)
    if numdays < 10:
        parser.error("End date must be greater than start date by at least 10 days.")
    if numdays < 30:
        threshhold_locations = numdays / 2
    else:
        threshhold_locations = math.log(numdays) * numdays / 10
    logger.info("Threshhold locations: %d" % threshhold_locations)
        
    afobj = ActiveFire(options.dsn)
    cur = afobj.execute(options.sql0, dict(
            start_date=sdate.strftime('%Y-%m-%d'),
            end_date=edate.strftime('%Y-%m-%d'),
            start_hour=options.start_hour,
            end_hour=options.end_hour,
        ))
    logger.debug("query returned %d rows" % cur.rowcount)

    sites_dict = {}
    sites = defaultdict(tuple)
    for row in cur:
        lng = float(row[0])
        lat = float(row[1])
        sites[(lng,lat)] = (lng,lat)

    for site in sites: 
        #for each site put it in a dictionary with its value being an array of neighbors 
        sites_dict[site] = [x for x in sites if x != site and sitesDist(site,x) < threshhold_dist_m]

    results = {}
    for site in sites: 
        j = len(sites_dict[site])
        if j >= threshhold_locations:
            #coord = bounding_box( site, sites_dict[site] )
            #results[coord] = coord
            results[site] = sites_dict[site]
    
    count = 0
    for site, lnglats in results.items():
        wkt = "MULTIPOINT(%s)" % (",".join([ "%s %s"%(x) for x in lnglats ]))
        logger.debug(wkt)
        if options.sql1:
            cur = afobj.execute(options.sql1, dict(
                    start_date=sdate.strftime('%Y-%m-%d'),
                    end_date=edate.strftime('%Y-%m-%d'),
                    wkt=wkt,
                ))
        count += 1
    if options.sql1:
        afobj.commit()
    logger.info("Found %d clusters" % count)

if __name__ == "__main__":
    sys.exit(main())