[gdal-dev] gdal python script to calculate categories statistic

[Anton Korosov]

Hi Giuseppe,

gdal_calc is built on operations with matrices using the Python module 
'numpy' http://numpy.scipy.org/

Instead of digging into it I would suggest to study numpy functionality 
a bit. The following simple script somewhat solves your problem:

import numpy as np

"""
you may data from GDAL dataset using e.g.:
ds = gdal.Open('filename_with_values_and_categories.tif')
valuesRaster = ds.GetRasterBand(1).b.ReadAsArray()
categoryRaster = ds.GetRasterBand(2).b.ReadAsArray()
but below I just suggest arrays with random values
"""
# array with FLOAT random values normally distributed
valuesRaster = np.random.randn(100, 300)
# array with INT random number uniformly distributed within [0, 9]
categoryRaster = np.uint8(np.random.random(valuesRaster.shape) * 10)
# list of all possible categories values
categories = np.unique(categoryRaster)
# loop over all categories
for category in categories:
     # take subsample from array with values based on array with categories
     subSample = valuesRaster[categoryRaster == category]
     # calculate statistics
     sumVal = np.sum(subSample)
     countVal = len(subSample)
     meanVal = np.mean(subSample)
     stdVal = np.std(subSample)
     # print to screen
     print 'Stats for category %d : %f %d %f %f' % (category, sumVal, 
countVal, meanVal, stdVal)


Anton

On 04/17/2012 09:48 PM, Giuseppe Amatulli wrote:
> Hi,
> I'm new in python and i would like to start in the hard way by modify
> an already existing script.
>
> In particular, I'm searching for a python script that is computing
> statistic of a raster base on another existing category raster.
> Something similar to r.statistics in grass.
> The out put has to be a raster with the statistical value of each
> categorise and also the summarized text file.
>
> Probably the most easy one is to modify is the gdal_calc.py but i'm
> open to other suggestions.
> Thanks in advance
>