# Processing/LIDARDEMs: lidar_histogram.py

File lidar_histogram.py, 4.5 KB (added by benj, 16 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 | |

7 | import re |

8 | import math |

9 | import sys |

10 | |

11 | HEIGHT_COL = 3 |

12 | NORTH_COL = 2 |

13 | EAST_COL = 1 |

14 | |

15 | # Function to print an (approximate) ascii histogram |

16 | # Uses log of histogram values |

17 | def 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 |

37 | def 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 | |

79 | if (len(sys.argv) != 2): |

80 | print "Usage: python lidar_histogram.py input_filename" |

81 | sys.exit(1) |

82 | # end if |

83 | |

84 | filepath = sys.argv[1] |

85 | |

86 | try: |

87 | datafile = open(filepath) |

88 | except: |

89 | print "Could not open input file: " + filepath |

90 | sys.exit(1) |

91 | # end try |

92 | |

93 | histogram = dict() # Empty dictionary to hold output |

94 | boundaries = dict() # Empty dictionary to hold map boundaries |

95 | |

96 | eof = False |

97 | |

98 | pointcounter = 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 |

102 | while (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 | |

154 | datafile.close() |

155 | |

156 | print_cutoffs(histogram, pointcounter) |

157 | print_histogram(histogram, boundaries, pointcounter) |