from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.core import *
from qgis.gui import *
from qgis.utils import *
from qgis.analysis import *
from qgis.networkanalysis import *
import sys

def run_script(iface):
    print sys.version
    
    layer = iface.activeLayer()
    lid = layer.getLayerID()
    print "Layer: %s" % lid
    lct = layer.featureCount()
    print "Features: %d" % lct
    prov = layer.dataProvider()
    feat = QgsFeature()
    attr = prov.attributeIndexes()
    amap = prov.fieldNameMap()
    imap = {v:k for k, v in amap.items()}
    
    for (k,attr) in imap.iteritems():
        print "Attribute %d: %s" % (k,imap[k])
    prov.select(attr)
    maxfeat = 10
    curfeat = 0
    rad = 500
    field = "COVERAGE"
    fix = prov.fieldNameIndex(field)
    while prov.nextFeature(feat):
        curfeat += 1
        if curfeat == maxfeat:
            break
        geom = feat.geometry()
        cent = geom.centroid().asPoint()
        x = cent.x()
        y = cent.y()
        rect = QgsRectangle(x-rad,y-rad,x+rad,y+rad)
        #freq = QgsFeatureRequest().setFilterRect(rect)
        #wind = layer.getFeatures(freq)
        layer = iface.activeLayer()
        wind = iface.activeLayer().getFeatures(QgsFeatureRequest().setFilterRect(rect))
        total = 0
        while wind.nextFeature(wfeat):
            attrs = wfeat.attributeMap()
            total += attrs[fix]
        print "Total in window: %d" % total
##        print "%d, %d" % (x,y)
        
        
##    print type(layer)
##        print type(geom)
##        print str(geom)
##        print "Feature ID %d@: "% feat.id() 
##        if 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)
##        attrs = feat.attributeMap()
##        for (k,attr) in attrs.iteritems():
##            print "%d: %s" % (k, attr.toString())