Processing/LIDARDEMs: lidar_histogram.py

File lidar_histogram.py, 4.5 KB (added by benj, 11 years ago)
Line 
1# Script to read a lidar file and output a (numeric) histogram for use in filtering
2# Author: Ben Taylor
3# Date: 15th April 2008
4# Copyright 2008 Plymouth Marine Laboratory - you may use this software as you wish,
5# but no warranty of any kind is given
6
7import re
8import math
9import sys
10
11HEIGHT_COL = 3
12NORTH_COL = 2
13EAST_COL = 1
14
15# Function to print an (approximate) ascii histogram
16# Uses log of histogram values
17def print_histogram(disp_histogram, bounds, totalpoints):
18    keylist = disp_histogram.keys()
19    keylist.sort()
20    print "Region boundaries:"
21    print "n: " + str(bounds["n"])
22    print "s: " + str(bounds["s"])
23    print "e: " + str(bounds["e"])
24    print "w: " + str(bounds["w"])
25    print "-------------------"
26    print "Log scale x5 histogram"
27    print "Total points: " + str(totalpoints)
28    for key in keylist:
29        value = disp_histogram[key]
30        logval = math.log10(value)
31        dispval = int(math.floor(5*logval)) # Scale value by 5
32        print str(key) + ": " + ("#"*dispval) + " (" + str(value) + ")"
33    # end for
34# end function
35
36# Function to work out recommended upper/lower bounds for the calculated histogram
37def print_cutoffs(disp_histogram, totalpoints):
38   
39    # -1 denotes not yet set (negative value are legal but will be multiples of 10 if valid)
40    lbound = -1
41    ubound = -1
42   
43    # Sort list of keys
44    keylist = disp_histogram.keys()
45    keylist.sort()
46   
47    # Run through list. First bracket containing more than 1/10000th of all points is minimum, less 10m
48    for key in keylist:
49        value = disp_histogram[key]
50       
51        if (lbound == -1):
52            if (value > 0.0001 * totalpoints):
53                lbound = key - 10
54                break
55            # end if           
56        # end if
57    # end for
58   
59    # Run through list again backwards. First bracket containing more than 1/10000th of all points is maximum,
60    # plus 10m (20 = width of bracket + 10m)
61    keylist.reverse()
62    for key in keylist:
63        value = disp_histogram[key]
64       
65        if (ubound == -1):
66            if (value > 0.0001 * totalpoints):
67                ubound = key + 20
68                break
69            # end if           
70        # end if
71    # end for
72   
73    print "Recommended bounds:"
74    print "Upper: " + str(ubound)
75    print "Lower: " + str(lbound)
76    print "-------------------"
77# end function
78
79if (len(sys.argv) != 2):
80    print "Usage: python lidar_histogram.py input_filename"
81    sys.exit(1)
82# end if
83
84filepath = sys.argv[1]
85
86try:
87    datafile = open(filepath)
88except:
89    print "Could not open input file: " + filepath
90    sys.exit(1)
91# end try
92
93histogram = dict() # Empty dictionary to hold output
94boundaries = dict() # Empty dictionary to hold map boundaries
95
96eof = False
97
98pointcounter = 0
99
100# For each file line, read it, split on whitespace, read the entries for height, northing and easting,
101# increment the relevant 10m height bracket in the histogram and update the map boundaries if appropriate
102while (not(eof)):
103    dataline = datafile.readline()
104    if (dataline == ""):
105        eof = True
106        break
107    # end if
108   
109    pointcounter = pointcounter + 1
110   
111    linesplit = re.split("\s+", dataline)
112    lineheight = float(linesplit[HEIGHT_COL])
113    linenorth = float(linesplit[NORTH_COL])
114    lineeast = float(linesplit[EAST_COL])
115   
116    heightindex = int(math.floor(lineheight) - (math.floor(lineheight) % 10)) # Get height bracket (round down to nearest 10m)
117   
118    if (heightindex in histogram):
119        histogram[heightindex] = histogram[heightindex] + 1
120    else:
121        histogram[heightindex] = 1
122    # end if
123   
124    if not ("n" in boundaries):
125        boundaries["n"] = int(math.ceil(linenorth))
126    # end if
127   
128    if not ("s" in boundaries):
129        boundaries["s"] = int(math.floor(linenorth))
130    # end if
131   
132    if not ("e" in boundaries):
133        boundaries["e"] = int(math.ceil(lineeast))
134    # end if
135   
136    if not ("w" in boundaries):
137        boundaries["w"] = int(math.floor(lineeast))
138    # end if
139   
140    if (linenorth > boundaries["n"]):
141        boundaries["n"] = int(math.ceil(linenorth))
142    elif (linenorth < boundaries["s"]):
143        boundaries["s"] = int(math.floor(linenorth))
144    # end if
145   
146    if (lineeast > boundaries["e"]):
147        boundaries["e"] = int(math.ceil(lineeast))
148    elif (lineeast < boundaries["w"]):
149        boundaries["w"] = int(math.floor(lineeast))
150    # end if
151       
152# end while
153
154datafile.close()
155
156print_cutoffs(histogram, pointcounter)
157print_histogram(histogram, boundaries, pointcounter)