# Issue Trying to create Zonal Statistics using Gdal and Python.

I am trying to create zonal statistics using Python and Gdal. I have a polygon shapefile and a raster file, and in order to do so, I am using a piece of code I found in StackExchange.

The raster and shapefile used can be found here

Here is the code I am using:

``````#!/usr/bin/python
#coding: utf-8

## Code from Stack Exchange: http://gis.stackexchange.com/questions/21888/how-to-overlay-shapefile-and-raster
## User: ustroetz

# 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()
targetSR.ImportFromWkt(raster.GetProjectionRef())
coordTrans = osr.CoordinateTransformation(sourceSR,targetSR)
feat = lyr.GetNextFeature()
geom = feat.GetGeometryRef()
geom.Transform(coordTrans)

# Get extent of geometry
ring = geom.GetGeometryRef(0)
numpoints = ring.GetPointCount()
pointsX = []; pointsY = []
for p in range(numpoints):
lon, lat, z = ring.GetPoint(p)
pointsX.append(lon)
pointsY.append(lat)
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)
target_ds.SetGeoTransform((
xmin, pixelWidth, 0,
ymax, 0, pixelHeight,
))

# create for target raster the same projection as for the value raster
raster_srs = osr.SpatialReference()
raster_srs.ImportFromWkt(raster.GetProjectionRef())
target_ds.SetProjection(raster_srs.ExportToWkt())

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

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

# calculate mean of zonal raster
return numpy.mean(zoneraster)
``````

However, when I implement it, the output (numpy.mean(zoneraster)) is just a number and not an array or a raster. What am I missing? Is there something wrong with the files imput? The code?

• The output is supposed to be only a number. The number is the mean of the raster of the area covered by your shapefile. What are you expecting to be in the array? You can find an updated version of my script here. Commented Nov 19, 2013 at 19:53
• I thought the code would create the mean value of my input raster for each polygon of my input shapefile. Therefore I was expecting to have, as an output something of the flavour: `[['polygon_label1', 'mean_raster1'],...,['polygon_labeln',mean_rastern']]` And by the way, thank you for posting that code! Commented Nov 20, 2013 at 10:55
• If you have many polygons, you need to rasterize them by some ID and then iterate over the unique IDs. You can accumulate the statistical output as a look-up table which you can join back to your zonal polygons. If speed becomes an issue as a result of the iteration use SciPy.ndimage instead of numpy. It is a couple of orders of magnitude faster in my experience. Commented Nov 20, 2013 at 11:10
• Do you know how I could loop over each polygon of the shapefile? Commented Nov 20, 2013 at 12:27
• I was looking for a script to do exactly this task, but, when trying it I've realized that averages are calculated over the whole bbox defined by each polygon in the shape layer. gdal.RasterizeLayer(target_ds, [1], lyr, burn_values=[1]) I don't really know if I am missunderstanding the code (I have seen the same approach in the oficial python gdal/ogr cookbook), or if it's in fact an error. However, any idea on how to parse just the feature instead of the whole layer in the mentione line? Many thanks in advance Commented Sep 9, 2014 at 14:59

With multiprocessing, for fastness! Has a little different output-formatting.

``````#!/usr/bin/python
import gdal, ogr, osr, numpy, sys
from multiprocessing import Pool

# Raster dataset
input_value_raster = sys.argv[1]

# Vector dataset(zones)
input_zone_polygon = sys.argv[2]

# Open data
rast = gdal.Open(input_value_raster)
shp = ogr.Open(input_zone_polygon)

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

def proc(fid):

# Each process needs its own pointer it seems.
shape = ogr.Open(input_zone_polygon)
lyr = shape.GetLayer()
feat = lyr.GetFeature(fid)
raster = gdal.Open(input_value_raster)

# Get extent of feat
geom = feat.GetGeometryRef()
if (geom.GetGeometryName() == 'MULTIPOLYGON'):
count = 0
pointsX = []; pointsY = []
for polygon in geom:
geomInner = geom.GetGeometryRef(count)
ring = geomInner.GetGeometryRef(0)
numpoints = ring.GetPointCount()
for p in range(numpoints):
lon, lat, z = ring.GetPoint(p)
pointsX.append(lon)
pointsY.append(lat)
count += 1
elif (geom.GetGeometryName() == 'POLYGON'):
ring = geom.GetGeometryRef(0)
numpoints = ring.GetPointCount()
pointsX = []; pointsY = []
for p in range(numpoints):
lon, lat, z = ring.GetPoint(p)
pointsX.append(lon)
pointsY.append(lat)

else:
sys.exit()

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)
target_ds.SetGeoTransform((
xmin, pixelWidth, 0,
ymax, 0, pixelHeight,
))

# Create for target raster the same projection as for the value raster
raster_srs = osr.SpatialReference()
raster_srs.ImportFromWkt(raster.GetProjectionRef())
target_ds.SetProjection(raster_srs.ExportToWkt())

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

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

# Calculate statistics of zonal raster
value = numpy.max(zoneraster)
#print value
return str(fid)+": "+str(value)

# Start the processes
layer = shp.GetLayer()
featList = range(layer.GetFeatureCount())
pool = Pool(processes=24)
print pool.map(proc,featList,8)

``````

If you want to get zonal statistics for several features in one shapefile, you have to loop over the `zonal_stats` function. You can write the results of the loop for example to a dictionary. Below is the modified `zonal_stats` function together with a loop, looping over the input shapefile. As an output you get a Dictionary containing for each Feature ID the mean of the covered raster.

For your sample dataset the dictionary for the first five features looks like this

``````{0: 114.57909872798288, 1: 21.889622561136882, 2: 287.79623686237102, 3: 35.350804240486838, 4: 19.63043032511781}
``````

The two functions:

``````import gdal, ogr, osr, numpy

def zonal_stats(feat, input_zone_polygon, input_value_raster):

# Open data
raster = gdal.Open(input_value_raster)
shp = ogr.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]

# Get extent of feat
geom = feat.GetGeometryRef()
if (geom.GetGeometryName() == 'MULTIPOLYGON'):
count = 0
pointsX = []; pointsY = []
for polygon in geom:
geomInner = geom.GetGeometryRef(count)
ring = geomInner.GetGeometryRef(0)
numpoints = ring.GetPointCount()
for p in range(numpoints):
lon, lat, z = ring.GetPoint(p)
pointsX.append(lon)
pointsY.append(lat)
count += 1
elif (geom.GetGeometryName() == 'POLYGON'):
ring = geom.GetGeometryRef(0)
numpoints = ring.GetPointCount()
pointsX = []; pointsY = []
for p in range(numpoints):
lon, lat, z = ring.GetPoint(p)
pointsX.append(lon)
pointsY.append(lat)

else:
sys.exit()

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)
target_ds.SetGeoTransform((
xmin, pixelWidth, 0,
ymax, 0, pixelHeight,
))

# Create for target raster the same projection as for the value raster
raster_srs = osr.SpatialReference()
raster_srs.ImportFromWkt(raster.GetProjectionRef())
target_ds.SetProjection(raster_srs.ExportToWkt())

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

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

# Calculate statistics of zonal raster
return numpy.mean(zoneraster)

def loop_zonal_stats(input_zone_polygon, input_value_raster):

shp = ogr.Open(input_zone_polygon)
lyr = shp.GetLayer()
featList = range(lyr.GetFeatureCount())
statDict = {}

for FID in featList:
feat = lyr.GetFeature(FID)
meanValue = zonal_stats(feat, input_zone_polygon, input_value_raster)
statDict[FID] = meanValue
return statDict

# Raster dataset
• I have another issue with the code in the section `#get extent of geometry`. My features can be multiple polygons. Here's my output for `geom.GetGeomRef(i)`, where `i' is each value of `geom.GetGeometryCount()`: `DATABASE: 7 LINEARRING (0.980468987882861 8.28964924042946, ..., 0.980468987882861 8.28964924042946) DATABASE: 8 POLYGON ((43.148330027378911 11.721389243497216,..., 11.72443924429683,43.148330027378911 11.721389243497216))` I modified the code to deal with these special cases (code written in next comment) but typically, it won't deal with polygons like "Database 8". Commented Nov 20, 2013 at 18:55
• My modified section is as follows: `# Get extent of geometry geonumbers=geom.GetGeometryCount() xmin_all=[] (...) ymax_all=[] for i in range(geonumbers): ring_i = geom.GetGeometryRef(i) numpoints = ring_i.GetPointCount() print ring_i pointsX_i = []; pointsY_i = [] for p in range(numpoints): lon, lat, z = ring_i.GetPoint(p) pointsX_i.append(lon) pointsY_i.append(lat) xmin_all.append(min(pointsX_i)) xmax_all.append(max(pointsX_i)) ymin_all.append(min(pointsY_i)) ymax_all.append(max(pointsY_i)) xmin=min(xmin_all) xmax=max(xmax_all) ymin=min(ymin_all) ymax=max(ymax_all) ` Commented Nov 20, 2013 at 18:57