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 firstname.lastname@example.org 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 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