Extracting shapefile attributes from raster colors using PyQGIS?

I have two layers:

• a shapefile with states and counties (shp)


• a raster file (with colors and some text and borders, but not a perfect map) (tiff)

I want to use the colors in the raster file to associate values to one of the attributes of the shapefile. I've started to do so manually with qgis. I overlayed the shapefile and the raster, and I read the color and fill in the value in the attribute editor:

I was wondering if this can be done programmatically. The pseudocode would be something like:

for county in polygons_in_shapefile:

    center_of_polygon = county.center()
    color = raster_color(center_of_polygon.x, center_of_polygon.y)
    if color == blue:
       county.value = 1
    if color == red:
       county.value = 2
    ... etc ...

I am totally new when it comes to GIS, so maybe it's very obvious (and the code may be totally the wrong strategy).

Answer

The can use raster.dataProvider().identify() (Using Raster Layers):

raster.dataProvider().identify(QgsPoint(x,y),QgsRaster.IdentifyFormatValue).results()

The result is a Python dictionary with the band number as key (R,G,B), and the values

# number of bands (raster)
print raster.bandCount()
3
# iterating over the layer
for county in polygons_in_shapefile.getFeatures():
    # geometry of county

    geom  = county.geometry()
    # centroid of each county
    x,y = geom.centroid().asPoint()
    # raster value 
    print raster.dataProvider().identify(QgsPoint(x,y), QgsRaster.IdentifyFormatValue).results()

# 1 = Red, 2 = Blue, 3 = Green
{1: 189.0, 2: 120.0, 3: 63.0}
{1: 215.0, 2: 204.0, 3: 148.0}
....

And you can create an "universal" function:

def Val_raster(point,raster):
    return raster.dataProvider().identify(point, QgsRaster.IdentifyFormatValue).results().values()

for county in polygons_in_shapefile.getFeatures():
     point = county.geometry().centroid().asPoint()
     print Val_raster(point,raster)
[189.0, 120.0, 63.0]
[215.0, 204.0, 148.0]

....

I there is only one band (DEM for example, the result is the z value):

print DEM.bandCount()
1
for county in polygons_in_shapefile.getFeatures():
     point = county.geometry().centroid().asPoint()
     print Val_raster(point,DEM)
[343.0]
[352.0]

....