Histogram of an ascii raster

March 6, 2007

Warning, esoteric GIS related post follows.

I often like to be able to generate a quick histogram of an ASCII raster file that I use for some of my models. I can do this inside of ArcGIS on a GRID file, but I like to be able to do this directly on an ASCII sometimes. Python to the rescue…

In a few lines of python I whipped this up:

#!/usr/bin/python
# By Jonah Duckles jonah@purdue.edu 7/17/2006
# Histogram Ascii will take a standard ESRI ascii file and output a histogram
#   of the values found in the file.

import sys, string, os
from operator import *

file = sys.argv[1]
histo = dict()
fh = open(file)
currow = 0
for line in fh.xreadlines():
    if currow > 6:
        line = string.rstrip(line, '/n')
        row = line.split(" ")
        row = row[:len(row)-1]
        for cell in row:
            if cell in histo:
                histo[cell] = histo[cell] + 1
            else:
                histo[cell] = 1
    currow += 1

k = histo.keys()
k.sort()

print "#Histogram for " + file
for element in k:
    print str(element) + "\t" + str(histo[element])

Now what if I want to be able to do this directly on a GRID file on Linux with no ArcGIS. I admit this is a complete hack, but here it is. A script that calls GDAL bin tools to convert the GRID to ASCII (it does it pretty fast actually), then calls hasc.py above and saves the histogram to a file.

#!/bin/bash
tmp=00-$1.asc
out=hist-$1.csv
gdal_translate -of AAIGrid $1 $tmp
echo "Computing histogram"
hasc $tmp > $out
echo "Histogram output to $out"
echo "Cleaning up"
rm $tmp

EDIT: Fixed a minor bug per pyther’s comment.