Antarctica, Winter 1982 ([info]ivorjawa) wrote,

decoding census geographic shape codes using QGIS and Python

Posted in response to this reddit discussion: http://www.reddit.com/r/programming/comments/7z8al/where_can_one_get_zip_code_boundary_polygon_data/


These are just some rough notes:

# click on layer

>>> l = iface.activeLayer()
>>> l
<qgis.core.QgsVectorLayer object at 0x1aa15300>

>>> l.name()
montana_poly

>>> l.source()
/Users/kujawa/Projects/state/montana_poly

# extent rectangle is in lat/lon

>>> e = l.extent()
>>> e
<qgis.core.QgsRectangle object at 0x1aa156f0>

>>> e.xMinimum()
-116.050003
>>> e.yMinimum()
44.358221
>>> e.yMaximum()
49.00139
>>> e.xMaximum()
-104.039563

>>> dp = l.dataProvider()
>>> dp
<qgis.core.QgsVectorDataProvider object at 0x1aa15738>

>>> feat = QgsFeature()
>>> dp.nextFeature(feat)
True

>>> geom = feat.geometry()

>>> x = geom.asPolygon()
>>> len(x)
1

#x is in array of arrays of tuples (all polygons?)
>>> poly0 = x[0]
#poly0 is an array of tuples

>>> v = poly0[0]
>>> v
(-116.049,49.0009)

#v is a vertix tuple of (lon, lat)





shell script to set up the environment:

#!/bin/sh
export PYTHONPATH="$PYTHONPATH:/Applications/qgis-1.0.0/qgis1.0.0.app/Contents/MacOS/share/qgis/python/"
export DYLD_FRAMEWORK_PATH="DYLD_FRAMEWORK_PATH:/Applications/qgis-1.0.0/qgis1.0.0.app/Contents/MacOS/lib/"
export DYLD_LIBRARY_PATH="DYLD_LIBRARY_PATH:/Applications/qgis-1.0.0/qgis1.0.0.app/Contents/MacOS/lib/"


python script to run it

#!/usr/bin/env python2.5
# state data from here: 
# http://www.census.gov/geo/www/cob/st2000.html

from qgis.core import *
QgsApplication.setPrefixPath("/Applications/qgis-1.0.0/qgis1.0.0.app/Contents/MacOS/", True)
QgsApplication.initQgis()

def go():
	source = "/Users/kujawa/Projects/state/st99_d00.shp"
	vlayer = QgsVectorLayer(source, "st99_d00_frak", "ogr")
	provider = vlayer.dataProvider()
	
	feat = QgsFeature()
	#allAttrs = provider.allAttributesList()
	
	# start data retreival: fetch geometry and all attributes for each feature
	#provider.select(allAttrs)

	fields = provider.fields()
	field_names = {}
	field_nums = []
	print "FIELDS:"
	for (k,attr) in fields.iteritems():
		field_nums.append(k)
		field_names[k] = attr.name()
	   	print "  %d: %s" % (k, attr.name())
	print


	provider.select(field_nums)
	
	# retreive every feature with its geometry and attributes
	while provider.nextFeature(feat):
		# fetch geometry
		geom = feat.geometry()
		print "Feature ID %d, type %s: " % (feat.id(), feat.typeName() ),

		# show some information about the feature
		if geom.type() == QGis.Point:
		  x = geom.asPoint()
		  print "Point: " + str(x)
		elif geom.type() == QGis.Line:
		  x = geom.asPolyline()
		  print "Line: %d points" % len(x)
		elif geom.type() == QGis.Polygon:
		  x = geom.asPolygon()
		  numPts = 0
		  for ring in x:
		    numPts += len(ring)
		  print "Polygon: %d rings with %d points" % (len(x), numPts)
		else:
		  print "Unknown"

		# fetch map of attributes
		attrs = feat.attributeMap()
		
		# attrs is a dictionary: key = field index, value = QgsFeatureAttribute
		# show all attributes and their values
		for (k,attr) in attrs.iteritems():
		   print "  %s: %s" % (field_names[k], attr.toString())	
		print
		
if __name__ == "__main__":
	go()


Tags: code

  • Post a new comment

    Error

    Your reply will be screened

    Your IP address will be recorded 

  • 0 comments
Create an Account
Forgot your login or password?
Facebook Twitter More login options
English • Español • Deutsch • Русский…