bin/wkt2pop

   1 #!/usr/bin/python
   2
   3 """
   4 Find the total population contained within a geometric object.
   5 """
   6
   7 """
   8 Our input is an OGC Well-Known Text[1] string. This string is used as
   9 part of a database query that finds the population contained within
  10 (i.e. 'underneath') the geometric object corresponding to the WKT
  11 string.
  12
  13 [1] http://en.wikipedia.org/wiki/Well-known_text
  14 """
  15
  16 import sys
  17 import os
  18 import site
  19 import pgdb
  20 from optparse import OptionParser
  21
  22 # Basically, add '../src' to our path.
  23 # Needed for the imports that follow.
  24 site.addsitedir(os.path.dirname(os.path.abspath(sys.argv[0])) + '/../src')
  25
  26 import Configuration.Defaults
  27 import ExitCodes
  28
  29 usage = '%prog [options] <well-known text representation>'
  30
  31 # -h (help) Conflicts with -h HOSTNAME
  32 parser = OptionParser(usage=usage, add_help_option = False)
  33
  34 # Use this module's docstring as the description.
  35 parser.description = __doc__.strip()
  36
  37 parser.add_option('-h',
  38                   '--host',
  39                   help='The hostname/address where the database is located.',
  40                   default=Configuration.Defaults.DATABASE_HOST)
  41
  42 parser.add_option('-d',
  43                   '--database',
  44                   help='The database in which the population data are stored.',
  45                   default=Configuration.Defaults.DATABASE_NAME)
  46
  47 parser.add_option('-U',
  48                   '--username',
  49                   help='The username who has access to the database.',
  50                   default=Configuration.Defaults.DATABASE_USERNAME)
  51
  52 parser.add_option('-s',
  53                   '--srid',
  54                   type="int",
  55                   help="SRID of the input geometry. Defaults to %s." % Configuration.Defaults.SRID,
  56                   default=Configuration.Defaults.SRID)
  57
  58
  59 (options, args) = parser.parse_args()
  60
  61 if len(args) < 1:
  62     print "\nERROR: You must supply a geometric object in Well-Known Text format.\n"
  63     parser.print_help()
  64     print '' # Print a newline.
  65     raise SystemExit(ExitCodes.NOT_ENOUGH_ARGS)
  66
  67
  68 conn = pgdb.connect(host=options.host,
  69                     database=options.database,
  70                     user=options.username)
  71
  72 cursor = conn.cursor()
  73
  74 # We're ready to build our query, one step at a time. Firsy, we store
  75 # the Text->Geom conversion in a variable; this just makes the query a
  76 # little easier to read.
  77 geometric_object = "ST_GeomFromText('%s', %d)" % (args[0], options.srid)
  78
  79 # We want to compute the population "under" the geometric object. We
  80 # can compute the percentage of a block that is covered by taking the
  81 # area of (the intersection of the object and the block) divided by
  82 # the total area of the block.
  83 #
  84 # Once we know the percentage covered, we just multiply that value by
  85 # the total population in the block to find the population that is
  86 # covered. The sum of these values over all blocks is our final
  87 # result.
  88 #
  89 query = """
  90 SELECT SUM(sf1_blocks.pop100 *
  91              ( ST_Area(ST_Intersection(%s, tiger.the_geom))
  92              / ST_Area(tiger.the_geom) )
  93           ) AS covered_population
  94 """ % geometric_object
  95
  96
  97 # Join our two block tables, so that we have both the demographic
  98 # and geometric data.
  99 query += \
 100 """
 101 FROM (sf1_blocks INNER JOIN tiger
 102       ON sf1_blocks.tiger_blkidfp00 = tiger.blkidfp00)
 103 """
 104
 105
 106 # We only need to calculate the covered population for the blocks
 107 # that actually intersect our object.
 108 query += \
 109 """
 110 WHERE (ST_Intersects(%s, tiger.the_geom))
 111 """  % geometric_object
 112
 113
 114 # And we only take the first result, since they're all going to be the
 115 # same (our query returns the sum once for each block).
 116 query += \
 117 """
 118 LIMIT 1
 119 """
 120
 121 cursor.execute(query)
 122 rows = cursor.fetchall()
 123
 124 if len(rows) > 0:
 125     population = rows[0][0]
 126     print population
 127 else:
 128     print 'Error: No rows returned.'
 129     raise SystemExit(ExitCodes.NO_RESULTS)
 130
 131 conn.close()