#!/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)

