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:

enter image description here

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).

  • In the raster image are the pixel values consistent? For example, do all of the 'green' pixels have exactly the same values? N.
    – nhopton
    Apr 21 '14 at 8:15
  • @nhopton, yes all greens have the same color, all blues have the same color, etc. (however, some lines are blurry, and there are letters and numbers on top of the map which could lead to: color = black for some points)
    – Massagran
    Apr 21 '14 at 16:00

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


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()
# 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()
for county in polygons_in_shapefile.getFeatures():
     point = county.geometry().centroid().asPoint()
     print Val_raster(point,DEM)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.