from struct import unpack import dbfUtils, math XY_POINT_RECORD_LENGTH = 16 db = [] def loadShapefile(file_name): global db shp_bounding_box = [] shp_type = 0 file_name = file_name records = [] # open dbf file and get records as a list dbf_file = file_name[0:-4] + '.dbf' dbf = open(dbf_file, 'rb') db = list(dbfUtils.dbfreader(dbf)) dbf.close() fp = open(file_name, 'rb') # get basic shapefile configuration fp.seek(32) shp_type = readAndUnpack('i', fp.read(4)) shp_bounding_box = readBoundingBox(fp) # fetch Records fp.seek(100) while True: shp_record = createRecord(fp) if shp_record == False: break records.append(shp_record) return records 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 == '': print 'doner' return False content_length = readAndUnpack('>L', fp.read(4)) record_shape_type = readAndUnpack(' maxarea: maxarea = ringArea biggest = ring #now get the true centroid tempPoint = {'x':0, 'y':0} if biggest[points][0] != biggest[points][len(biggest[points])-1]: print "mug", biggest[points][0], biggest[points][len(biggest[points])-1] for i in range(0, len(biggest[points])-1): j = (i + 1) % (len(biggest[points])-1) tempPoint['x'] -= (biggest[points][i]['x'] + biggest[points][j]['x']) * ((biggest[points][i]['x'] * biggest[points][j]['y']) - (biggest[points][j]['x'] * biggest[points][i]['y'])) tempPoint['y'] -= (biggest[points][i]['y'] + biggest[points][j]['y']) * ((biggest[points][i]['x'] * biggest[points][j]['y']) - (biggest[points][j]['x'] * biggest[points][i]['y'])) tempPoint['x'] = tempPoint['x'] / ((6) * maxarea) tempPoint['y'] = tempPoint['y'] / ((6) * maxarea) feature['shp_data']['truecentroid'] = tempPoint def getArea(ring, points): #returns the area of a polygon #needs to be spherical area, but isn't area = 0 for i in range(0,len(ring[points])-1): j = (i + 1) % (len(ring[points])-1) area += ring[points][i]['x'] * ring[points][j]['y'] area -= ring[points][i]['y'] * ring[points][j]['x'] return math.fabs(area/2) def getNeighbors(records): #for each feature for i in range(len(records)): #print i, records[i]['dbf_data']['ADMIN_NAME'] if not 'neighbors' in records[i]['shp_data']: records[i]['shp_data']['neighbors'] = [] #for each other feature for j in range(i+1, len(records)): numcommon = 0 #first check to see if the bounding boxes overlap if overlap(records[i], records[j]): #if so, check every single point in this feature to see if it matches a point in the other feature #for each part: for part in records[i]['shp_data']['parts']: #for each point: for point in part['points']: for otherPart in records[j]['shp_data']['parts']: if point in otherPart['points']: numcommon += 1 if numcommon == 2: if not 'neighbors' in records[j]['shp_data']: records[j]['shp_data']['neighbors'] = [] records[i]['shp_data']['neighbors'].append(j) records[j]['shp_data']['neighbors'].append(i) #now break out to the next j break if numcommon == 2: break if numcommon == 2: break def projectShapefile(records, whatProjection, lonCenter=0, latCenter=0): print 'projecting to ', whatProjection for feature in records: for part in feature['shp_data']['parts']: part['projectedPoints'] = [] for point in part['points']: tempPoint = projectPoint(point, whatProjection, lonCenter, latCenter) part['projectedPoints'].append(tempPoint) def projectPoint(fromPoint, whatProjection, lonCenter, latCenter): latRadians = fromPoint['y'] * math.pi/180 if latRadians > 1.5: latRadians = 1.5 if latRadians < -1.5: latRadians = -1.5 lonRadians = fromPoint['x'] * math.pi/180 lonCenter = lonCenter * math.pi/180 latCenter = latCenter * math.pi/180 newPoint = {} if whatProjection == "MERCATOR": newPoint['x'] = (180/math.pi) * (lonRadians - lonCenter) newPoint['y'] = (180/math.pi) * math.log(math.tan(latRadians) + (1/math.cos(latRadians))) if newPoint['y'] > 200: newPoint['y'] = 200 if newPoint['y'] < -200: newPoint['y'] = 200 return newPoint if whatProjection == "EQUALAREA": newPoint['x'] = 0 newPoint['y'] = 0 return newPoint def overlap(feature1, feature2): if (feature1['shp_data']['xmax'] > feature2['shp_data']['xmin'] and feature1['shp_data']['ymax'] > feature2['shp_data']['ymin'] and feature1['shp_data']['xmin'] < feature2['shp_data']['xmax'] and feature1['shp_data']['ymin'] < feature2['shp_data']['ymax']): return True else: return False