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.
comments powered by Disqus