Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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'm a total noob when it comes to GIS, so maybe it's very obvious (and the code may be totally the wrong strategy). Any ideas, tutorials or links would be appreciated.

share|improve this question
    
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 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 at 16:00

1 Answer 1

up vote 3 down vote accepted

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]
....
share|improve this answer
    
Thank you. I didn't know about those scripting options with Qgis. –  Massagran Apr 22 at 5:11

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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