Identify Census blocks inside urban areas (using pyshp or similar)

Question

I am trying to find the Census blocks that are inside specific urban areas or urban clusters.

The Census Bureau provides a list of urban areas/urban clusters with the following information:

| Column Name | Column Description | | ------------- | ----------------------------------------------------------------- | | UACE | Five digit code that uniquely identifies 2010 urban area entities | | NAME | 2010 urban area name | | POP | 2010 population of the urban area | | HU | 2010 housing unit count for the urban area | | AREALAND | Land area for the urban area in square meters | | AREALANDSQMI | Land area for the record converted to square miles | | AREAWATER | Water area for the record in square meters | | AREAWATERSQMI | Water area for the record converted to square miles | | POPDEN | Population density | | LSADC | Legal/Statistical area description code for | | | the urban area (75 for urbanized area, 76 for urban cluster) |

The first urban area in this file is Abbeville, LA (UACE=00037). As an example, I would like to find all the Census Blocks that are in this urban area.

As far as I can tell, there is not an API or other method for getting a list of Census Blocks from a specified urban area, so I am trying to find alternatives.

I am an experienced programmer but have no GIS experience. Theoretically, it seems like it should be possible to take the shapefile for urban areas and the shapefile for Census Blocks in Vermilion Parish (where Abbeville is located), overlay the two, and find all the Census Blocks that are covered by the Abbeville urban area polygon.

I need to do this a number of times so I was hoping to do this programmatically with pyshp or similar. But I am open to any solution that I can automate and isn't buying ArcGIS.

Answer 1

I have written a set of scripts to solve this problem and have published the output in case anyone else needs it. See the bottom of this answer for details.

GeospatialPython.com pointed me in the right direction. I was able to throw together a quick Python script to do this. I'm sure there are issues here because I'm not very familiar with pyqgis, but it does appear to work (here is how to verify it through the QGIS GUI.

from qgis.core import *
# Note: add /Applications/QGIS.app/Contents/Resources/python/ to the PYTHONPATH environment variable for this to import.

# supply path to where is your qgis installed
QgsApplication.setPrefixPath("/Applications/QGIS.app/Contents/MacOS", True)

# load providers
QgsApplication.initQgis()

# If you have trouble loading layers, see this: https://gis.stackexchange.com/questions/59069/creating-qgis-layers-in-python-console-vs-stand-alone-application

urban_areas_layer = QgsVectorLayer("data/census/urban_areas/tl_2010_us_uac10.shp", "urban_areas", "ogr")
if not urban_areas_layer.isValid():
  print "Layer failed to load!"

block_layer = QgsVectorLayer("data/census/blocks/tl_2010_22113_tabblock10/tl_2010_22113_tabblock10.shp", "blocks", "ogr")
if not block_layer.isValid():
  print "Layer failed to load!"

# Gets the names of the fields in a layer
# https://gis.stackexchange.com/questions/76364/how-to-get-field-names-in-pyqgis-2-0
field_names = [field.name() for field in urban_areas_layer.pendingFields() ]

field_names_blocks = [field.name() for field in block_layer.pendingFields() ]

for urban_area in urban_areas_layer.getFeatures():
  attributes =  dict(zip(field_names, urban_area.attributes()))
  if attributes['UACE10'] == '00037':
    block_count = 0
    print attributes

    for block in block_layer.getFeatures():
      if block.geometry().touches(urban_area.geometry()) == True:
        block_attributes =  dict(zip(field_names_blocks, block.attributes()))
        block_count += 1
        # print block_attributes['GEOID10']
    print "%d block(s) found" % block_count

QgsApplication.exitQgis()

Edit 1:

Here is Python code that fully automates the process of downloading the county shapefiles needed for each urban area. The output should be a CSV file listing each Census Block that is in an urban area.

There are some prerequisites for this code:

• qgis Python package installed (it is part of the QGIS GUI install)
• pandas package
• Create these nested folders off your working directory for county shapefiles to go: data/census/blocks/
• Download tl_2010_us_uac10.shp to data/census/blocks/urban_areas/
• Download ua_county_rel_10.txt to data/census/blocks/urban_areas/

import pandas as pd
import os # for running shell commands
from qgis.core import * # Note: add /Applications/QGIS.app/Contents/Resources/python/ to the PYTHONPATH environment variable for this to import.


# supply path to where is your qgis installed
QgsApplication.setPrefixPath("/Applications/QGIS.app/Contents/MacOS", True)

# load providers
QgsApplication.initQgis()

# If you have trouble loading layers, see this: https://gis.stackexchange.com/questions/59069/creating-qgis-layers-in-python-console-vs-stand-alone-application

urban_areas_layer = QgsVectorLayer("data/census/urban_areas/tl_2010_us_uac10.shp", "urban_areas", "ogr")
if not urban_areas_layer.isValid():
  print "Layer failed to load!"



# Valid state numerical codes
# From https://www.census.gov/geo/reference/ansi_statetables.html
valid_states = [1, 2, 4, 5, 6, 8, 9, 10, 11, 12, 13, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 44, 45, 46, 47, 48, 49, 50, 51, 53, 54, 55, 56]

# Relationships between urban areas and counties
# http://www2.census.gov/geo/docs/maps-data/data/rel/ua_county_rel_10.txt
urban_areas_counties = pd.read_csv('data/census/urban_areas/urban_area_county.txt', quotechar='"')


# Gets the names of the fields in a layer
# https://gis.stackexchange.com/questions/76364/how-to-get-field-names-in-pyqgis-2-0
field_names = [field.name() for field in urban_areas_layer.pendingFields() ]


# Set up output dataframe
columns = ['ua_name', 'ua_id', 'state', 'county', 'tract', 'block', 'block_geoid']
output_df = pd.DataFrame(columns=columns)

for urban_area in urban_areas_layer.getFeatures():
  block_count = 0
  attributes =  dict(zip(field_names, urban_area.attributes()))
  print "Starting urban area %s..." % attributes['NAME10']

  # Look up counties that this UA covers
  counties = urban_areas_counties[(urban_areas_counties.UA == int(float(attributes['UACE10'])))]

  # Skip any counties that aren't in the US
  counties = counties[counties['STATE'].isin(valid_states)]

  # Check to see if shapefiles are already downloaded
  for county in counties.iterrows():
    print "    Starting county %s" % county[1]['CNAME']

    # Figure out filename for shapefile for county
    zero_pad_state = '%02d' % int(float(county[1]['STATE']))
    zero_pad_county = '%03d' % int(float(county[1]['COUNTY']))
    filename = 'tl_2010_%s%s_tabblock10' % (zero_pad_state, zero_pad_county)

    # Download shapefile for county if it hasn't been downloaded already
    if not os.path.exists("data/census/blocks/%s/%s.shp" % (filename, filename)):
      cmd = "cd data/census/blocks && rm -rf %s && mkdir %s && cd %s && curl -O ftp://ftp2.census.gov//geo/tiger/TIGER2010/TABBLOCK/2010/%s.zip && unzip %s.zip" % (filename, filename, filename, filename, filename)
      os.system(cmd)

    # Load in block layer for county
    block_layer = QgsVectorLayer("data/census/blocks/%s/%s.shp" % (filename, filename), "blocks", "ogr")
    if not block_layer.isValid():
      print "Layer %s.shp to load!" % filename
    else:
      # Get field names for block
      field_names_blocks = [field.name() for field in block_layer.pendingFields() ]

      for block in block_layer.getFeatures():
        if block.geometry().within(urban_area.geometry()) == True:
          block_attributes =  dict(zip(field_names_blocks, block.attributes()))
          block_count += 1
          # print "        Block found: %s" % block_attributes['GEOID10']

          # Need to dump the following attributes:
          #   STATEFP10
          #   COUNTYFP10
          #   TRACTCE10
          #   BLOCKCE10
          #   GEOID10
          #
          # Along with this, save the urban area info:
          #   attributes['NAME10']
          #   attributes['UACE10']
          #
          # columns = ['ua_name', 'ua_id', 'state', 'county', 'tract', 'block', 'block_geoid']

          data = [{
            'ua_name': attributes['NAME10'],
            'ua_id': attributes['UACE10'],
            'state': block_attributes['STATEFP10'],
            'county': block_attributes['COUNTYFP10'],
            'tract': block_attributes['TRACTCE10'],
            'block': block_attributes['BLOCKCE10'],
            'block_geoid': block_attributes['GEOID10']
          }]
          output_df = output_df.append(data)


          # code.interact(local=locals())
          # sys.exit()
  print "%d block(s) found for %s" % (block_count, attributes['NAME10'])

output_df.to_csv('data/census/urban_area_blocks.csv', index = False)

# cleanup
QgsApplication.exitQgis()

---

I have released all the code for doing this on GitHub, and have published the output (a single CSV file with a row for each block that is in an urban area/urban cluster, along with that block's population) here. These scripts/data are not verified so use at your own risk. Pull requests welcome.

Answer 2

I'm the author of PyShp. What you are trying to do is very standard GIS stuff. Download QGIS. It is basically the open-source ArcGIS. It bundles all of the best of breed open source GIS tools into a nice GUI as well as a fantastic Python automation framework. There are dozens of free tutorials and videos for performing different tasks. There also lots of good books. Drag your shapefiles (.shp) onto the canvas and prototype it in the GUI and then batch automate using the Processing a Toolbox. In QGIS speak (and the underlying OGR library) you want to see what blocks an urban area "contains". You can run your automated script in QGIS or you can turn it into a command line program using OGR which has the contains() method if you need it to run outside of QGIS. Regardless you'll be using Python.