Processing/LIDARDEMs: lidar_histogram.py

# Script to read a lidar file and output a (numeric) histogram for use in filtering
# Author: Ben Taylor
# Date: 15th April 2008
# Copyright 2008 Plymouth Marine Laboratory - you may use this software as you wish,
# but no warranty of any kind is given

import re
import math
import sys

HEIGHT_COL = 3
NORTH_COL = 2
EAST_COL = 1

# Function to print an (approximate) ascii histogram
# Uses log of histogram values
def print_histogram(disp_histogram, bounds, totalpoints):
    keylist = disp_histogram.keys()
    keylist.sort()
    print "Region boundaries:"
    print "n: " + str(bounds["n"])
    print "s: " + str(bounds["s"])
    print "e: " + str(bounds["e"])
    print "w: " + str(bounds["w"])
    print "-------------------"
    print "Log scale x5 histogram"
    print "Total points: " + str(totalpoints)
    for key in keylist:
        value = disp_histogram[key]
        logval = math.log10(value)
        dispval = int(math.floor(5*logval)) # Scale value by 5
        print str(key) + ": " + ("#"*dispval) + " (" + str(value) + ")"
    # end for
# end function

# Function to work out recommended upper/lower bounds for the calculated histogram
def print_cutoffs(disp_histogram, totalpoints):
    
    # -1 denotes not yet set (negative value are legal but will be multiples of 10 if valid)
    lbound = -1
    ubound = -1
    
    # Sort list of keys
    keylist = disp_histogram.keys()
    keylist.sort()
    
    # Run through list. First bracket containing more than 1/10000th of all points is minimum, less 10m
    for key in keylist:
        value = disp_histogram[key]
        
        if (lbound == -1):
            if (value > 0.0001 * totalpoints):
                lbound = key - 10
                break
            # end if            
        # end if
    # end for
    
    # Run through list again backwards. First bracket containing more than 1/10000th of all points is maximum, 
    # plus 10m (20 = width of bracket + 10m)
    keylist.reverse()
    for key in keylist:
        value = disp_histogram[key]
        
        if (ubound == -1):
            if (value > 0.0001 * totalpoints):
                ubound = key + 20
                break
            # end if            
        # end if
    # end for
    
    print "Recommended bounds:"
    print "Upper: " + str(ubound)
    print "Lower: " + str(lbound)
    print "-------------------"
# end function

if (len(sys.argv) != 2):
    print "Usage: python lidar_histogram.py input_filename"
    sys.exit(1)
# end if

filepath = sys.argv[1]

try:
    datafile = open(filepath)
except:
    print "Could not open input file: " + filepath
    sys.exit(1)
# end try

histogram = dict() # Empty dictionary to hold output
boundaries = dict() # Empty dictionary to hold map boundaries

eof = False

pointcounter = 0

# For each file line, read it, split on whitespace, read the entries for height, northing and easting,
# increment the relevant 10m height bracket in the histogram and update the map boundaries if appropriate
while (not(eof)):
    dataline = datafile.readline()
    if (dataline == ""):
        eof = True
        break
    # end if
    
    pointcounter = pointcounter + 1
    
    linesplit = re.split("\s+", dataline)
    lineheight = float(linesplit[HEIGHT_COL])
    linenorth = float(linesplit[NORTH_COL])
    lineeast = float(linesplit[EAST_COL])
    
    heightindex = int(math.floor(lineheight) - (math.floor(lineheight) % 10)) # Get height bracket (round down to nearest 10m)
    
    if (heightindex in histogram):
        histogram[heightindex] = histogram[heightindex] + 1
    else:
        histogram[heightindex] = 1
    # end if
    
    if not ("n" in boundaries):
        boundaries["n"] = int(math.ceil(linenorth))
    # end if
    
    if not ("s" in boundaries):
        boundaries["s"] = int(math.floor(linenorth))
    # end if
    
    if not ("e" in boundaries):
        boundaries["e"] = int(math.ceil(lineeast))
    # end if
    
    if not ("w" in boundaries):
        boundaries["w"] = int(math.floor(lineeast))
    # end if
    
    if (linenorth > boundaries["n"]):
        boundaries["n"] = int(math.ceil(linenorth))
    elif (linenorth < boundaries["s"]):
        boundaries["s"] = int(math.floor(linenorth))
    # end if
    
    if (lineeast > boundaries["e"]):
        boundaries["e"] = int(math.ceil(lineeast))
    elif (lineeast < boundaries["w"]):
        boundaries["w"] = int(math.floor(lineeast))
    # end if
        
# end while

datafile.close()

print_cutoffs(histogram, pointcounter)
print_histogram(histogram, boundaries, pointcounter)