In R you can do
library(raster)
library(rgdal)
r <- raster('raster_filename')
p <- readOGR('shp_path', 'shp_file')
e <- extract(r, p, fun=mean)
e is a vector with the mean of the raster cell values for each polygon.
QgsGeometry.fromMultiPolygon() expects a list of polygons. A polygon, in turn, is a list of rings: An outer ring, and a number of inner rings. A ring, in turn, is a list of QgsPoint objects.
A valid MultiPolygon could thus look like this:
[ # The list of polygons
[ # The list of rings for the first polygon
[ # The list of QgsCoordinates
QgsPoint(0,0), QgsPoint(1,0), QgsPoint(0,1)
]
],
[ # The list of rings for the second polygon
[ # The list of QgsCoordinates for the first polygon
QgsPoint(2,2), QgsPoint(3,2), QgsPoint(2,3)
]
]
]
You can transform your list of coordinates into that structure using
[[[QgsPoint(point[0],point[1]) for point in polygon ] for polygon in polygon_points ]]
So to create the feature, you could use
feature.setGeometry(QgsGeometry.fromMultiPolygon([[[QgsPoint(point[0],point[1]) for point in polygon ] for polygon in polygon_points ]]))
Best Answer
See the PyQGIS Cookbook for advise how to add attributes to vector layers:
http://www.qgis.org/pyqgis-cookbook/vector.html#adding-and-removing-fields
However the easier way in your situation would be to do a simple spatial join to add your point values to the polygons.