#!/usr/bin/python

"""
Find the total population contained within a geometric object.
"""

"""
Our input is an OGC Well-Known Text[1] string. This string is used as
part of a database query that finds the population contained within
(i.e. 'underneath') the geometric object corresponding to the WKT
string.

[1] http://en.wikipedia.org/wiki/Well-known_text
"""

import sys
import os
import site
import pgdb
from optparse import OptionParser

# Basically, add '../src' to our path.
# Needed for the imports that follow.
site.addsitedir(os.path.dirname(os.path.abspath(sys.argv[0])) + '/../src')

import Configuration.Defaults
import ExitCodes

usage = '%prog [options] <well-known text representation>'

# -h (help) Conflicts with -h HOSTNAME
parser = OptionParser(usage=usage, add_help_option = False)

# Use this module's docstring as the description.
parser.description = __doc__.strip()

parser.add_option('-h',
                  '--host',
                  help='The hostname/address where the database is located.',
                  default=Configuration.Defaults.DATABASE_HOST)

parser.add_option('-d',
                  '--database',
                  help='The database in which the population data are stored.',
                  default=Configuration.Defaults.DATABASE_NAME)

parser.add_option('-U',
                  '--username',
                  help='The username who has access to the database.',
                  default=Configuration.Defaults.DATABASE_USERNAME)

parser.add_option('-s',
                  '--srid',
                  help="SRID of the input geometry. Defaults to %s." % Configuration.Defaults.SRID,
                  default=Configuration.Defaults.SRID)


(options, args) = parser.parse_args()

if len(args) < 1:
    print "\nERROR: You must supply a geometric object in Well-Known Text format.\n"
    parser.print_help()
    print '' # Print a newline.
    raise SystemExit(ExitCodes.NOT_ENOUGH_ARGS)


conn = pgdb.connect(host=options.host,
                    database=options.database,
                    user=options.username)

cursor = conn.cursor()

# We're ready to build our query, one step at a time. Firsy, we store
# the Text->Geom conversion in a variable; this just makes the query a
# little easier to read.
geometric_object = "ST_GeomFromText('%s', %s)" % (args[0], options.srid)

# We want to compute the population "under" the geometric object. We
# can compute the percentage of a block that is covered by taking the
# area of (the intersection of the object and the block) divided by
# the total area of the block.
#
# Once we know the percentage covered, we just multiply that value by
# the total population in the block to find the population that is
# covered. The sum of these values over all blocks is our final
# result.
#
query = """
SELECT SUM(sf1_blocks.pop100 *
             ( ST_Area(ST_Intersection(%s, tiger.the_geom))
             / ST_Area(tiger.the_geom) )
          ) AS covered_population
""" % geometric_object


# Join our two block tables, so that we have both the demographic
# and geometric data.
query += \
"""          
FROM (sf1_blocks INNER JOIN tiger
      ON sf1_blocks.tiger_blkidfp00 = tiger.blkidfp00)
"""


# We only need to calculate the covered population for the blocks
# that actually intersect our object.
query += \
"""
WHERE (ST_Intersects(%s, tiger.the_geom))
"""  % geometric_object


# And we only take the first result, since they're all going to be the
# same (our query returns the sum once for each block).
query += \
"""
LIMIT 1
"""

cursor.execute(query)
rows = cursor.fetchall()

if len(rows) > 0:
    population = rows[0][0]
    print population
else:
    print 'Error: No rows returned.'
    raise SystemExit(ExitCodes.NO_RESULTS)
    
conn.close()