Alternative raster stats implementations

Alternatives:

ArcPY

>>> import arcpy
>>> arcpy.sa.ZonalStatisticsAsTable("fcid", "FID", "dem_aea2_feet.tif", "E:\workspace\stats1.dbf")
<geoprocessing server result object object at 0x059F1FC8>

34 seconds outputs dbf format, return values not accessible as python variables

QGIS

    from qgis.analysis import QgsZonalStatistics

    #specify polygon shapefile vector
    polygonLayer = QgsVectorLayer('F:/temp/zonalstat/zonePoly.shp', 'zonepolygons', "ogr") 

    # specify raster filename
    rasterFilePath = 'F:/temp/zonalstat/raster1.tif'

    # usage - QgsZonalStatistics (QgsVectorLayer *polygonLayer, const QString &rasterFile, const QString &attributePrefix="", int rasterBand=1)
    zoneStat = QgsZonalStatistics (polygonLayer, rasterFilePath, 'pre-', 1)
    zoneStat.calculateStatistics(None)

Slower (in gui took 111 sec) limited stats, variables not accessible, alters source dataset

R

library(raster)
library(rgdal)
library(foreign)

# 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")

Starspan

Grass

Postgis 2.0+