imageryqgisrasterstatistics
i need to calculate raster statistics for RGBN (4 channels, multiband image) rasters in QGIS.
In Arcgis if you use "calculate statistics" in ArcCatalog, statistics are calculated for every image band.
In Qgis if you use "Raster layer statistics" you get a log file with a summary of image stats, so you cannot extract statistics for each band.
Do you know how this is achieved in Qgis?
I was struggling to do exactly the same thing, but for various reasons I'm committed to using QGIS. I tried using v.rast.stats using the GRASS plugin and also via the Sextante plugin. The latter approach failed, because it seems to attach the stats to a temporary vector layer which it then deletes. The GRASS plugin worked, but it doesn't deal with overlapping polygons.
After doing some digging around (in the source of the promising-sounding ZonalStats plugin), I found that QGIS actually has zonal statistics methods built into the API, and these also have Python bindings. So as long as you are only looking for count, sum and mean statistics for your polygon features, the Python Console (Plugins > Python Console) is currently the easiest way to attach the stats to the polygon attribute table.
Select your raster layer in the TOC, and type the following into the console (it grabs the source file name of your raster layer)
>>> rasterfile = qgis.utils.iface.mapCanvas().currentLayer().source()
Select your vector layer, and execute the following command in the console (it grabs the vector layer itself)
>>> vectorlayer = qgis.utils.iface.mapCanvas().currentLayer()
Execute the following three commands in the console (they pass the vector layer and raster file to QGIS' built-in zonal stats calculator)
>>> import qgis.analysis
>>> zonalstats = qgis.analysis.QgsZonalStatistics(vectorlayer,rasterfile)
>>> zonalstats.calculateStatistics(None)
The results will be appended as extra fields in the polygon layer.
Note that if you want to append a prefix to the new fields (rather than just having them as sum, mean and count), then you would replace the first command in step 3 with the following (assuming you wanted to use 'zonal_' as the prefix:
>>> zonalstats = qgis.analysis.QgsZonalStatistics(vectorlayer,rasterfile,"zonal_")
The following code works for me (in that QGis reads the statistics correctly). Note use of ComputeStatistics instead of GetStatistics as that will force gdal to recalculate band statistics.
driver = gdal.GetDriverByName( 'GTiff' )
dst_ds = driver.Create( dst_filename, 512,512, 1, gdal.GDT_UInt16 )
raster = somenumpy.calculation(etc...)
rb=dst_ds.GetRasterBand(1)
rb.WriteArray( raster )
rb.ComputeStatistics(0)
dst_ds = None
Best Answer
One way to get statistics for each band in QGIS raster layer is:
-statsflag to a command (betweengdalinfoand path to raster)Run a command.
As a result you will see log with simple statistics for each band.