Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I have a shapefile with polygons. And I have a global raster file. I want to overlay the shapefile's polygons onto the raster grid and calculate the mean raster value for each polygon.

How can I do this using GDAL, writing the results to the shapefile?

share|improve this question
Is GDAL the only tool you'd like to use? – Simbamangu Mar 18 '12 at 17:15
@Simbamangu no, basically anything is fine, and it would be great if it were in Python – andreash Apr 16 '12 at 15:25

In R you can do

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.

share|improve this answer
that's a really nice and easy one - thanks! – andreash Jun 20 '12 at 16:26
up vote 6 down vote accepted

Following advice I got on the gdal-dev mailing list, I used StarSpan:

starspan --vector V --raster R1 R2 ... --stats mystats.csv avg mode

The results are saved to CSV format. At that time, that was already enough for me, but it should be possible to somehow forge a Shapefile from that info.

share|improve this answer
StarSpan seems to have moved to GitHub. Get it here. – Richard Nov 27 '14 at 5:05

Load your shapefile and your raster in PostGIS 2.0 and do:

SELECT (ST_SummaryStats(ST_Clip(rast, geom))).*
FROM rastertable, geomtable
share|improve this answer

I don't think GDAL is the best tool for this, but you can use gdal_rasterize to "clear" all the values outside the polygon.

Something like:

gdal_translate -a_nodata 0 original.tif work.tif
gdal_rasterize -burn 0 -b 1 -i work.tif yourpolygon.shp -l yourpolygon
gdalinfo -stats work.tif
rm work.tif

The gdal_rasterize program modifies the file, so we make a copy to work on. We also mark some particular value (zero in this case) to be nodata. The "-burn 0 -b 1" means burn a value of zero into band 1 of the target file (work.tif). The "-i" means invert rasterization so we burn values outside the polygon instead of inside it. The gdalinfo command with -stats reports on band statistics. I believe it will exclude the nodata value (which we marked earlier with -a_nodata).

share|improve this answer

Transform the shape file in raster by gdal_rasterize and the use the code in to calculate zonal statistic for each polygon. You can run if you want obtain a tif with your rasters statistic. Enjoy the code Ciao Giuseppe

share|improve this answer
Do you happen to have a copy of the code you refer to? Unfortunately the link to the Python file is dead. – ustroetz May 21 '13 at 17:51

The following script allows you to do the task with GDAL:

# Calculates statistics (mean) on values of a raster within the zones of an polygon shapefile

import gdal, ogr, osr, numpy

def zonal_stats(input_value_raster, input_zone_polygon):

    # Open data
    raster = gdal.Open(input_value_raster)
    driver = ogr.GetDriverByName('ESRI Shapefile')
    shp = driver.Open(input_zone_polygon)
    lyr = shp.GetLayer()

    # get raster georeference info
    transform = raster.GetGeoTransform()
    xOrigin = transform[0]
    yOrigin = transform[3]
    pixelWidth = transform[1]
    pixelHeight = transform[5]

    # reproject geometry to same projection as raster
    sourceSR = lyr.GetSpatialRef()
    targetSR = osr.SpatialReference()
    coordTrans = osr.CoordinateTransformation(sourceSR,targetSR)
    feat = lyr.GetNextFeature()
    geom = feat.GetGeometryRef()

    # Get extent of geometry
    ring = geom.GetGeometryRef(0)
    numpoints = ring.GetPointCount()
    pointsX = []; pointsY = []
    for p in range(numpoints):
            lon, lat, z = ring.GetPoint(p)
    xmin = min(pointsX)
    xmax = max(pointsX)
    ymin = min(pointsY)
    ymax = max(pointsY)

    # Specify offset and rows and columns to read
    xoff = int((xmin - xOrigin)/pixelWidth)
    yoff = int((yOrigin - ymax)/pixelWidth)
    xcount = int((xmax - xmin)/pixelWidth)+1
    ycount = int((ymax - ymin)/pixelWidth)+1

    # create memory target raster
    target_ds = gdal.GetDriverByName('MEM').Create('', xcount, ycount, gdal.GDT_Byte)
        xmin, pixelWidth, 0,
        ymax, 0, pixelHeight,

    # create for target raster the same projection as for the value raster
    raster_srs = osr.SpatialReference()

    # rasterize zone polygon to raster
    gdal.RasterizeLayer(target_ds, [1], lyr, burn_values=[1])

    # read raster as arrays
    banddataraster = raster.GetRasterBand(1)
    dataraster = banddataraster.ReadAsArray(xoff, yoff, xcount, ycount).astype(numpy.float)

    bandmask = target_ds.GetRasterBand(1)
    datamask = bandmask.ReadAsArray(0, 0, xcount, ycount).astype(numpy.float)

    # mask zone of raster
    zoneraster =,  numpy.logical_not(datamask))

    # calculate mean of zonal raster
    return numpy.mean(zoneraster)
share|improve this answer
This is really useful, thank you so much – Doon_Bogan Nov 19 '13 at 11:42

This is not possible using GDAL. You could use other free tools however, eg saga gis:

saga_cmd shapes_grid "Grid Values to Shapes" -GRIDS=grid.sgrd -POLYGONS=in.shp -SHAPES=out.shp-NODATA -TYPE=1
share|improve this answer
I went with this approach, though the function name is actually "Grid Statistics for Polygons". – bananafish Jan 18 '13 at 4:46

Your Answer


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.