From: Michael Orlitzky Date: Thu, 24 Sep 2009 04:56:46 +0000 (-0400) Subject: Added the wkt2pop script. X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=8f026c4ceb6aed31c1b1ac0f306c1569eba000bd;p=dead%2Fcensus-tools.git Added the wkt2pop script. --- diff --git a/bin/wkt2pop b/bin/wkt2pop new file mode 100755 index 0000000..af35595 --- /dev/null +++ b/bin/wkt2pop @@ -0,0 +1,130 @@ +#!/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] ' + +# -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()