I am trying to use QGIS to extract raster statistics for shapefiles (I have shapefiles of species distribution and I want to extract environmental data from within each species' range). This is possible using the Zonal Statistics plugin however I have about 300 shapefiles I need to get data from. I've tried merging the shapefiles into one but this results in massive files that cause QGIS to crash (using shapefiles >15mB causes problems).

Is there a way to automate the procedure using a Python script? I found a post that details how to run ZonalStats using Python (How to calculate raster statistics for polygons? second answer) but can't get it to work. I know absolutely nothing about Python and so I'm not sure if there's something that needs to be done before attempting to run the script or if the problem is something else. A similar script that will run on its own without user input would be ideal.

I'm running QGIS 1.8 on a Linux machine.

Any help would be much appreciated.

  • "using shapefiles >15mB causes problems". Are you sure? I've successfully load at qgis (both in windows and linux) shapefiles which were over 100 mbytes
    – nickves
    Commented Dec 11, 2012 at 23:38
  • If you don't know anything about python scripting your next option is using GRASS, which you can very easily write a script which iterates over any number of shapefiles. Also I suggest to invest some time to learn python, because it is a really powerful tool for scientific analyses
    – nickves
    Commented Dec 11, 2012 at 23:44
  • Thanks for your comments Nick. Shapefiles >15mB cause problems for the Zonal Stats plugin - it takes hours to finish and then returns results for only 15 or so polygons. Running shapefiles one at a time (without merging into larger files) works fine however. I had a brief play with GRASS and found it a bit daunting but I'll look into it. If you could offer some guidance on exactly what it is I need to do in GRASS that would be a big help.
    – Thomas
    Commented Dec 12, 2012 at 15:50
  • @Thomas - Your issue is likely to be more about the underlying raster size than the shapefiles. Any process which does a cell-by-cell look-up on a large raster will be hammered at high resolutions. Doubling the resolition, quadruples the raster data so things can quickly get out of hand. SciPy handles the data differently. I recently processed over a hundred vectors averaging about 15mB deriving values from rasters averaging nearly 400mB each. SciPy took between a few seconds to about five minutes per file. That time saving alone makes it worth getting to grips with Python and SciPy. Commented Jan 4, 2013 at 21:32
  • In general, how long does this process take? I've been using very large raster files (TIF's around 700 mb) and a vector layer with about 200 polygons; it's been several hours and the process still hasn't completed. Sorry if this marked as an answer -- I didn't want to open a new thread to address the same topic, though. Commented Dec 18, 2013 at 20:06

3 Answers 3


Yes you can do this in a stand alone script. To do this you efficiently, you will need SciPy. SciPy is a cousin of NumPy and together they extend Python with a lot of extremely powerful array handling capabilities. In SciPy you will use the 'ndimage' module.

Your procedure is to make a raster of your zones. Then import your original raster and the zones raster each into a NumPy array. Make a third numpy array of the unique values in your zones array to use as an index for the ndimage function. Finally use one of the ndimage statistics to extract your result, e.g.

myIndex = numpy.unique(zonesArray)    
minResult = ndimage.minimum(myRasterArray, labels=zonesArray, index=indexArray)

(obviously I have not included all the code to open a raster or vector file etc. but see below for help with that).

I often write my results to a a simple csv text file and can then do a table join to pull the statistics back into my zone polygons. If you have many rasters, then you can automate the process into a batch process.

While you can do this with GDAL and indeed, will need GDAL to read the original raster and polygon data, I would not recommend that you attempt it without SciPy because the process will have glacial alacrity! See this excellent series on coding GDAL with Python. Note tutorial 6 in the series as this is especially relevant. However, it doesn't mention SciPy as it is a bit old (and where it mentions 'Numeric' understand that this has been replaced by 'NumPy'). ndimage is very fast and much better than manually iterating over a raster in a bitwise fashion or even iterating over an array cell by cell. The speed difference of using SciPy vs. NumPy alone can be seconds compared to hours for the same data (from personal experience). As a bonus, it takes less code too :)

There is a lot packed into this message and it may appear daunting if you know nothing of Python. In which case, start with the tutorial series I mention above, which assumes zero coding experience. Then after the first couple of lessons, skip forward to the raster material. Once you know how to open a raster and get it into a Numpy (aka 'numeric') array, read the SciPy documentation on ndimage.

  • Thanks for your answer. You're right, the stuff you mention does sound a bit daunting. I will keep this in mind as an option but I was really after a way to do this in QGIS. QGIS works well and I'm 90% of the way there, I just need to learn a bit of python.
    – Thomas
    Commented Dec 12, 2012 at 15:57
  • The tutorial series will really help with learning Python. Not only will it give you an introduction to Python, but it will also put it in the context of using the GDAL/OGR Python API - which is a double win. Commented Dec 13, 2012 at 8:21

In general: My QGIS (1.9 - Master ( ;-) ) on Debian Linux) works quite well with big shapefiles (hell, i have some that are over 300mb big). Maybe this is some problem with your conversion/merging process and if so you should switch to other tools (search for it on this website).

Regarding your data:

  • You didn't specify the nature of your "shapefiles". If you have a point layer, you could just use the Point sampling tool in QGIS (search for it via the plugin downloader)
  • If you have a shapefile with polygons and you're somehow experienced with statistical analyses in R, then see this little example of mine. Use the R-Code to extract any values from your raster.
  • If you strictly want a python-code than you should look into the gdal api and code it.
  • You are also encouraged to use the other very good tools available out there (even for linux). Both GRASS and SAGA have a function, which does exactly what you want. And if you're clever you could even use those functions directly in QGIS via the Sextante toolbox.
  • You're working with qgis master (pre 1.9 at the moment). QGIS 1.9 hasn't been released yet. I speculate it will be out in a couple of months.
    – nickves
    Commented Dec 11, 2012 at 23:40
  • Thanks for your answer Curlew. I am able to load large shapefiles, the problem is with the Zonal Stats plugin - shapefiles >15mB take hours and results for only a small fraction of polygons are returned (my shapefiles are polygons). I am familiar with R however until now have been unable to find an easy way to get R do what I want. I will look into your example, thanks. GRASS is a bit scary but may be worth a look. I've tried doing this in SAGA already - the results it returns are inconsistent i.e. running the same thing twice gives a different result each time.
    – Thomas
    Commented Dec 12, 2012 at 16:12
  • You might want to try the QGIS plugin i coded, LecoS. Starting with version 1.4 it is capable of extracting raster values per polygon-overlays using numpy and scipy. More here (tinyurl.com/cy6hs9q)
    – Curlew
    Commented Jan 5, 2013 at 13:55

The solution I was after has now been answered in another question I posted (How to use QGIS Zonal Stats Plugin from Python Console?). The actual solution I used however was in R, as suggested here by Curlew.

The code I used was as follows:


# extract raster data from overlying polygon
rast <- raster('raster.tif')
poly <- readOGR(dsn="polygons", layer="polygon")
ext.poly <- extract(rast, poly, fun = mean, na.rm=TRUE, df=TRUE)

# and then append the result to the .dbf file of the original polygon
poly.dbf <- read.dbf("polygons/polygon.dbf", as.is = TRUE)
poly.dbf$result <- ext.poly$extraction
write.dbf(poly.dbf, "polygon_with_extraction.dbf")

Thanks everyone for the help.

  • 1
    This seems quite a lot slower than QGIS Zonal Statistics in my case.
    – thadk
    Commented Apr 29, 2015 at 20:16
  • In my case too. Especially because rgdal has substantial problems to handle big shapefiles (100MB +). Commented Jan 25, 2016 at 9:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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