text: #!/usr/bin/env python
  
  # shpUtils.py
  # Original version by Zachary Forest Johnson
  # http://indiemaps.com/blog/index.php/code/pyShapefile.txt
  # This version modified by Michael Geary
  
  from struct import unpack
  import dbfUtils
  XY_POINT_RECORD_LENGTH = 16
  db = []
  
  def loadShapefile( filename ):
  	# open dbf file and get features as a list
  	global db
  	dbfile = open( filename[0:-4] + '.dbf', 'rb' )
  	db = list( dbfUtils.dbfreader(dbfile) )
  	dbfile.close()
  	
  	fp = open( filename, 'rb' )
  	
  	# get basic shapefile configuration
  	fp.seek(32)
  	filetype = readAndUnpack('i', fp.read(4))
  	bounds = readBounds( fp )
  	
  	# fetch Records
  	fp.seek(100)
  	features = []
  	while True:
  		feature = createRecord(fp)
  		if feature == False: break
  		getPolyInfo( feature )
  		features.append( feature )
  	return { 'type': filetype, 'bounds': bounds, 'features': features }
  
  record_class = { 0:'RecordNull', 1:'RecordPoint', 8:'RecordMultiPoint', 3:'RecordPolyLine', 5:'RecordPolygon' }
  
  def createRecord(fp):
  	# read header
  	record_number = readAndUnpack('>L', fp.read(4))
  	if record_number == '': return False
  	content_length = readAndUnpack('>L', fp.read(4))
  	rectype = readAndUnpack('<L', fp.read(4))
  
  	shape = readRecordAny(fp,rectype)
  	shape['type'] = rectype
  	info = {}
  	names = db[0]
  	values = db[record_number+1]
  	for i in xrange(len(names)):
  		value = values[i]
  		if isinstance( value, str ):
  			value = value.strip()
  		info[ names[i] ] = value
  	
  	return { 'shape':shape, 'info':info }
  	
  # Reading defs
  
  def readRecordAny(fp, rectype):
  	if rectype==0:
  		return readRecordNull(fp)
  	elif rectype==1:
  		return readRecordPoint(fp)
  	elif rectype==8:
  		return readRecordMultiPoint(fp)
  	elif rectype==3 or rectype==5:
  		return readRecordPolyLine(fp)
  	else:
  		return False
  
  def readRecordNull(fp):
  	return {}
  
  point_count = 0
  def readRecordPoint(fp):
  	global point_count
  	point = [ readAndUnpack('d', fp.read(8)), readAndUnpack('d', fp.read(8)) ]
  	point_count += 1
  	return point
  
  
  def readRecordMultiPoint(fp):
  	shape = { 'bounds': readBounds(fp) }
  	points = shape['points'] = []
  	nPoints = readAndUnpack('i', fp.read(4))
  	for i in xrange(nPoints):
  		points.append(readRecordPoint(fp))
  	return shape
  
  def readRecordPolyLine(fp):
  	shape = { 'bounds': readBounds(fp) }
  	nParts = readAndUnpack('i', fp.read(4))
  	nPoints = readAndUnpack('i', fp.read(4))
  	if readAndUnpack('i', fp.read(4)):
  		print 'ERROR: First part offset must be 0'
  	counts = []; prev = 0
  	for i in xrange(nParts-1):
  		next = readAndUnpack('i', fp.read(4))
  		counts.append( next - prev )
  		prev = next
  	counts.append( nPoints - prev )
  	parts = shape['parts'] = []
  	for i in xrange(nParts):
  		part = {}
  		parts.append( part )
  		points = part['points'] = []
  		for j in xrange(counts[i]):
  			points.append(readRecordPoint(fp))
  	return shape
  
  # General defs
  
  def readBounds(fp):
  	return [
  		[ readAndUnpack('d',fp.read(8)), readAndUnpack('d',fp.read(8)) ],
  		[ readAndUnpack('d',fp.read(8)), readAndUnpack('d',fp.read(8)) ]
  	]
  
  def readAndUnpack(fieldtype, data):
  	if data=='': return data
  	return unpack(fieldtype, data)[0]
  
  def getPolyInfo( feature ):
  	nPoints = cx = cy = 0
  	shape = feature['shape']
  	shapetype = shape['type']
  	if shapetype == 3  or shapetype ==5:
  		for part in shape['parts']:
  			getPartInfo( part )
  
  def getPartInfo( part ):
  	points = part['points']
  	n = len(points)
  	area = cx = cy = 0
  	xmin = ymin = 360
  	xmax = ymax = -360
  	pt = points[n-1];  xx = pt[0];  yy = pt[1]
  	for pt in points:
  		x = pt[0];  y = pt[1]
  		# bounds
  		xmin = min( x, xmin );  ymin = min( y, ymin )
  		xmax = max( x, xmax );  ymax = max( y, ymax )
  		# area and centroid
  		a = xx * y - x * yy
  		area += a
  		cx += ( x + xx ) * a
  		cy += ( y + yy ) * a
  		# next
  		xx = x;  yy = y
  	area /= 2
  	if area:
  		centroid = [ cx / area / 6, cy / area / 6 ]
  	else:
  		centroid = None
  	part.update({
  		'area': abs(area),
  		'bounds': [ [ xmin, ymin ], [ xmax, ymax ] ],
  		'center': [ ( xmin + xmax ) / 2, ( ymin + ymax ) / 2 ],
  		'centroid': centroid,
  		'extent': [ abs( xmax - xmin ), abs( ymax - ymin ) ]
  	})
  
  def getBoundCenters(features):
  	for feature in features:
  		bounds = feature['shape']['bounds']
  		min = bounds[0];  max = bounds[1]
  		bounds['center'] = [
  			( min[0] + max[0] ) / 2,
  			( min[1] + max[1] ) / 2
  		]
  	
  def getMAT(features):
  	print 'feature not yet available'
  
  def dumpFeatureInfo( features ):
  	fields = []
  	rows = []
  	for feature in features:
  		info = feature['info']
  		if not len(fields):
  			for key in info:
  				fields.append( key )
  			rows.append( ','.join(fields) )
  		cols = []
  		for field in fields:
  			cols.append( str( feature['info'][field] ) )
  		rows.append( ','.join(cols) )
  	return '\r\n'.join(rows)
  