How to create a circle vector layer with 12 sectors with Python and to save average value of underlying raster in each sector?

Question

Hi I am writing a code with Python 2.7 using Arcpy 10.1 and I guess I am doing stuck with a silly problem that lies under my eyes. I have created a grid of points (with "Create fishnet" function) over a region where I have also raster data (for instance land cover or DEM). For each point I need to associate a circle in vector format (with center in each point) divided in "n" sectors of 1 km radius for example. The goal is to calculate the mean value of all pixels of the land cover or of the DEM situated in each sector and to save them in the attribute table of the circles.

I found this very helpful post (How to create a circle vector layer with 12 sectors with Python/pyqgis?) but I am somehow stuck where I have to save the data in shapefile instead of geojson and to go further with the estimate of the mean values for each sector.

Here is part of the code take from the previous code:

# helper function to calculate point from relative polar coordinates (degrees)
def polar_point(origin_point, angle,  distance):
    return [origin_point.x + math.sin(math.radians(angle)) * distance, origin_point.y +math.cos(math.radians(angle)) * distance]

features = []

for x in xrange(0,int(sectors)):
    segment_vertices = []

    # first the center and first point
    segment_vertices.append(polar_point(center, x,y))
    segment_vertices.append(polar_point(center, start + x*sector_width,radius))

    # then the sector outline points
    for z in xrange(1, steps_per_sector):
        segment_vertices.append((polar_point(center, start + x * sector_width + z * step_angle_width,radius)))

    # then again the center point to finish the polygon
    segment_vertices.append(polar_point(center, start + x *  sector_width+sector_width,radius))
    segment_vertices.append(polar_point(center, x,y))

    # create feature
    features.append(geojson.Feature(geometry=Polygon(segment_vertices)))

# prepare geojson feature collection
res = geojson.FeatureCollection(features = features)

# write to file
f = open('sector.shp', 'w')
f.write(geojson.dumps(res))
f.close() }

Any piece of advice out there please? Is any other approach maybe more practical/feasible? I was thinking also about "Euclidean direction" and "zonal Statistics", but it deals with many raster data and the calculation might be very long in case of large regions to analyze.

Best regards.