I recently came across a loop_zonal_stats function by @ustroetz found here: Issue Trying to create Zonal Statistics using Gdal and Python.. I am using a (single) polygon shapefile and a raster file to compute the zonal statistics. The polygon shapefile (http://www.epa.gov/wed/pages/ecoregions/na_eco.htm#Downloads) is a list of North American ecoregions. I first opened the file in QGIS and saved it in the EPSG:4326 - WGS 84 coordinate reference system. The raster file (http://worldclim.org/current) contains precipitation data and is an ESRI grid file with 30 arc-second resolution.
One issue I noticed with the code is that NoData values are included in the calculations (a "bug" that I noticed in QGIS as well). In order to fix this, I added the following lines of code to the zonal_stats function:
nodata = banddataraster.GetNoDataValue()
if nodata is not None:
zoneraster = np.ma.masked_equal(zoneraster, nodata)
After making this change, I am successfully able to exclude the NoData values from the calculation. However, I noticed something strange that I don't understand. Many of the polygons in my shapefile have the same attribute (and none of them overlap). I used QGIS to convert the shapefile of polygons to a shapefile of multipolygons (QGIS SinglePart to MultiPart) and then applied the loop_zonal_stats function to compute the sum and count for this multipolygon shapefile. However, the results from the original polygon shapefile are inconsistent with those from the multipolygon shapefile. For example, for a particular region with "NA_L3CODE" = 10.1.2, I get a count of 352 + 6873 + 317400 = 324625 for the original polygon shapefile and a count of 369288 for the multipolygon shapefile. Might anyone know what is going on?