# Processing/LIDARDEMs: generate_mask.py

File generate_mask.py, 5.0 KB (added by benj, 13 years ago) |
---|

Line | |
---|---|

1 | import bz2 |

2 | import re |

3 | import math |

4 | import sys |

5 | |

6 | NORTH_COL = 2 |

7 | EAST_COL = 1 |

8 | PADDING = 2000 |

9 | |

10 | if (len(sys.argv) != 4): |

11 | print "Usage: python generate_mask.py file_list mask_file point_file" |

12 | sys.exit(1) |

13 | # end if |

14 | |

15 | dirlistfilename = sys.argv[1] |

16 | maskfilename = sys.argv[2] |

17 | merged_filename = sys.argv[3] |

18 | |

19 | try: |

20 | dirfile = open(dirlistfilename) |

21 | except: |

22 | print "Could not open input file list: " + dirlistfilename |

23 | sys.exit(1) |

24 | # end try |

25 | |

26 | try: |

27 | maskfile = open(maskfilename, "w") |

28 | except: |

29 | print "Could not open mask file for writing: " + maskfilename |

30 | sys.exit(1) |

31 | # end try |

32 | |

33 | try: |

34 | mergedfile = open(merged_filename, "w") |

35 | except: |

36 | print "Could not open point file for writing: " + merged_filename |

37 | sys.exit(1) |

38 | # end try |

39 | |

40 | all_bounds = list() |

41 | |

42 | enddir = False |

43 | while (not enddir): |

44 | north_max = -1 |

45 | south_max = -1 |

46 | east_max = -1 |

47 | west_max = -1 |

48 | nw_max = -1 |

49 | sw_max = -1 |

50 | ne_max = -1 |

51 | se_max = -1 |

52 | linebounds = dict() |

53 | |

54 | datafilename = dirfile.readline() |

55 | if (datafilename == ""): |

56 | enddir = True |

57 | break |

58 | else: |

59 | datafilename = datafilename.strip() # Strip newline character off end of line if necessary |

60 | # end if |

61 | |

62 | try: |

63 | datafile = bz2.BZ2File(datafilename) |

64 | except: |

65 | print "Could not open lidar input file: " + datafilename |

66 | sys.exit(1) |

67 | # end try |

68 | |

69 | for currentline in datafile: |

70 | |

71 | # Strip off zone number (17th-18th characters on each line) |

72 | strippedline = "" |

73 | for i in range(0, len(currentline)): |

74 | if ((i != 16) and (i != 17)): |

75 | strippedline = strippedline + currentline[i] |

76 | # end if |

77 | # end for |

78 | |

79 | strippedline = strippedline.strip() |

80 | pointfields = re.split("\s+", strippedline) # Split line on whitespace |

81 | point_north = float(pointfields[NORTH_COL]) |

82 | point_east = float(pointfields[EAST_COL]) |

83 | |

84 | if ((north_max == -1) or (north_max < int(math.ceil(point_north)) + PADDING)): |

85 | north_max = int(math.ceil(point_north)) + PADDING |

86 | nmax_point = (int(math.floor(point_east)), north_max) |

87 | # end if |

88 | |

89 | if ((south_max == -1) or (south_max > int(math.floor(point_north)) - PADDING)): |

90 | south_max = int(math.floor(point_north)) - PADDING |

91 | smax_point = (int(math.floor(point_east)), south_max) |

92 | # end if |

93 | |

94 | if ((east_max == -1) or (east_max < int(math.ceil(point_east)) + PADDING)): |

95 | east_max = int(math.ceil(point_east)) + PADDING |

96 | emax_point = (east_max, int(math.floor(point_north))) |

97 | # end if |

98 | |

99 | if ((west_max == -1) or (west_max > int(math.floor(point_east)) - PADDING)): |

100 | west_max = int(math.floor(point_east)) - PADDING |

101 | wmax_point = (west_max, int(math.floor(point_north))) |

102 | # end if |

103 | |

104 | # Max north-west point defined by 2*(ycoord-xcoord) (derived from geometry against SW/NE axis) |

105 | if ((nw_max == -1) or (nw_max < int(math.ceil(2*(point_north - point_east))))): |

106 | nw_max = int(math.ceil(2*(point_north - point_east))) |

107 | nwmax_point = (int(math.floor(point_east)-PADDING), int(math.ceil(point_north)+PADDING)) |

108 | # end if |

109 | |

110 | if ((se_max == -1) or (se_max > int(math.floor(2*(point_north - point_east))))): |

111 | se_max = int(math.floor(2*(point_north - point_east))) |

112 | semax_point = (int(math.floor(point_east)-PADDING), int(math.floor(point_north)-PADDING)) |

113 | # end if |

114 | |

115 | # Max north-east point defined by 2x(ycoord+xcoord) (derived from geometry against SE/NW axis) |

116 | if ((ne_max == -1) or (ne_max < int(math.ceil(2*(point_north + point_east))))): |

117 | ne_max = int(math.ceil(2*(point_north + point_east))) |

118 | nemax_point = (int(math.ceil(point_east)+PADDING), int(math.ceil(point_north)+PADDING)) |

119 | # end if |

120 | |

121 | if ((sw_max == -1) or (sw_max > int(math.floor(2*(point_north + point_east))))): |

122 | sw_max = int(math.floor(2*(point_north + point_east))) |

123 | swmax_point = (int(math.ceil(point_east)+PADDING), int(math.floor(point_north)-PADDING)) |

124 | # end if |

125 | |

126 | mergedfile.write(pointfields[0] + " " + pointfields[1] + " " + pointfields[2] + " " + pointfields[3] + " " + pointfields[4] + "\n") |

127 | # end for |

128 | |

129 | # Write points to mask file |

130 | maskfile.write(str(nmax_point[0]) + "," + str(nmax_point[1]) + "\n") |

131 | maskfile.write(str(nwmax_point[0]) + "," + str(nwmax_point[1]) + "\n") |

132 | maskfile.write(str(wmax_point[0]) + "," + str(wmax_point[1]) + "\n") |

133 | maskfile.write(str(swmax_point[0]) + "," + str(swmax_point[1]) + "\n") |

134 | maskfile.write(str(smax_point[0]) + "," + str(smax_point[1]) + "\n") |

135 | maskfile.write(str(semax_point[0]) + "," + str(semax_point[1]) + "\n") |

136 | maskfile.write(str(emax_point[0]) + "," + str(emax_point[1]) + "\n") |

137 | maskfile.write(str(nemax_point[0]) + "," + str(nemax_point[1]) + "\n") |

138 | # end while |

139 | |

140 | dirfile.close() |

141 | mergedfile.close() |

142 | maskfile.close() |