1

I'm trying to extract the count of each pixel value in a raster, like in this question: How to extract pixel value counts from a raster in QGIS? However, when I run r.report (from QGis (2.10) Processing toolbox) I get a result like this: r.report result

Even if my settings in r.report are the same as the ones in that question i linked earlier:

r.report settings

I tried to change the settings but I can't get the result I want. I don't know if I'm missing something obvious, but I've been looking around for a while and found no explanation for this.

  • could you clarify how the results differ from what you were after? Were you wanting a count for discrete raster values (20480, 20481, 20482... 611440), rather than by band? – Steven Kay Aug 13 '15 at 19:10
  • 1
    Yes, for discrete values. The raster I'm using is a quality assessment raster for Landsat data, so the pixels can only have a limited range of discrete values that correspond to cloud, cloud shadow, etc. I just want to know how many pixels belong to each value so I can decide if the image is usable for my purposes or not. I don't understand why it's giving me "from-to" results... – Simona Aug 13 '15 at 19:26
1

I don't know any way of doing this without writing custom code.. this is a use case I've come across in the past, but there doesn't seem to be an easy way to do it.

I'd be delighted if someone points out an easier/quicker solution!

There are several good stats based plugins, but they're for vector layers.

Zonal Statistics is good for finding summary stats for a raster (average, majority), but that won't help you here.

As your raster is using floats, you're likely to find a very uniform distribution unless you discretise the values, e.g. by rounding down to nearest int.

This chunk of Python code will round down the float values, then count each value. It should work if you have a recent release of QGIS, all the dependencies should be there. I've assumed it's a one-band float raster.

Be careful, this is a quick-and-dirty piece of code which loads the whole raster into memory, so it may crash if the raster is large. As with all bits of code on the internet, back up your work first!

#!/usr/bin/env python
# -*- coding: UTF-8 -*-

from osgeo import ogr
from osgeo import osr
from osgeo import gdal
from gdalconst import GA_ReadOnly

# top tip - don't include spaces in the path..it upsets gdal

ds = gdal.Open("/tmp/linear.tif",GA_ReadOnly)
band = ds.GetRasterBand(1) # first band
cols = ds.RasterXSize               
rows = ds.RasterYSize

# reads in whole raster at once. May run out of memory on big rasters!

data = band.ReadAsArray(0, 0, cols, rows)
counts = {}
rownum = 0

# tally up pixel values

for scanrow in data:
    print "Processing scan line %d" % rownum
    for pixelval in scanrow:
        pixelval = int(pixelval)
        if not pixelval in counts:
            counts[pixelval] = 1
        else:
            counts[pixelval] += 1
    rownum += 1

# now dump the values to console in CSV format 

print "Value,Frequency"
for pixelval in sorted(counts):
    print "%d,%d" % (pixelval,counts[pixelval])

Here's some example output from a DEM...

Value,Frequency
-2,251658
-1,163897
0,810264
1,4993
2,6912
3,5666
4,5984
5,7326
6,7861
7,9599
8,11716
9,14051
10,9079
11,7636
....
  • Thank you very much for your help! It looks like it should work, but unfortunately I have never used python and I don't know how to run codes... if you know of some very basic tutorial to get started with it I'd be grateful because I tried a few but I got stuck immediately (also, I don't know how up-to-date they are). – Simona Aug 14 '15 at 7:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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