From: Michael Orlitzky <michael@orlitzky.com>
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] <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()