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                   help="SRID of the input geometry. Defaults to %s." % Configuration.Defaults.SRID,
  55                   default=Configuration.Defaults.SRID)
  56
  57
  58 (options, args) = parser.parse_args()
  59
  60 if len(args) < 1:
  61     print "\nERROR: You must supply a geometric object in Well-Known Text format.\n"
  62     parser.print_help()
  63     print '' # Print a newline.
  64     raise SystemExit(ExitCodes.NOT_ENOUGH_ARGS)
  65
  66
  67 conn = pgdb.connect(host=options.host,
  68                     database=options.database,
  69                     user=options.username)
  70
  71 cursor = conn.cursor()
  72
  73 # We're ready to build our query, one step at a time. Firsy, we store
  74 # the Text->Geom conversion in a variable; this just makes the query a
  75 # little easier to read.
  76 geometric_object = "ST_GeomFromText('%s', %s)" % (args[0], options.srid)
  77
  78 # We want to compute the population "under" the geometric object. We
  79 # can compute the percentage of a block that is covered by taking the
  80 # area of (the intersection of the object and the block) divided by
  81 # the total area of the block.
  82 #
  83 # Once we know the percentage covered, we just multiply that value by
  84 # the total population in the block to find the population that is
  85 # covered. The sum of these values over all blocks is our final
  86 # result.
  87 #
  88 query = """
  89 SELECT SUM(sf1_blocks.pop100 *
  90              ( ST_Area(ST_Intersection(%s, tiger.the_geom))
  91              / ST_Area(tiger.the_geom) )
  92           ) AS covered_population
  93 """ % geometric_object
  94
  95
  96 # Join our two block tables, so that we have both the demographic
  97 # and geometric data.
  98 query += \
  99 """
 100 FROM (sf1_blocks INNER JOIN tiger
 101       ON sf1_blocks.tiger_blkidfp00 = tiger.blkidfp00)
 102 """
 103
 104
 105 # We only need to calculate the covered population for the blocks
 106 # that actually intersect our object.
 107 query += \
 108 """
 109 WHERE (ST_Intersects(%s, tiger.the_geom))
 110 """  % geometric_object
 111
 112
 113 # And we only take the first result, since they're all going to be the
 114 # same (our query returns the sum once for each block).
 115 query += \
 116 """
 117 LIMIT 1
 118 """
 119
 120 cursor.execute(query)
 121 rows = cursor.fetchall()
 122
 123 if len(rows) > 0:
 124     population = rows[0][0]
 125     print population
 126 else:
 127     print 'Error: No rows returned.'
 128     raise SystemExit(ExitCodes.NO_RESULTS)
 129
 130 conn.close()