Processing/LIDARDEMs: generate_mask.py

File generate_mask.py, 5.0 KB (added by benj, 11 years ago)
Line 
1import bz2
2import re
3import math
4import sys
5
6NORTH_COL = 2
7EAST_COL = 1
8PADDING = 2000
9
10if (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
15dirlistfilename = sys.argv[1]
16maskfilename = sys.argv[2]
17merged_filename = sys.argv[3]
18
19try:
20    dirfile = open(dirlistfilename)
21except:
22    print "Could not open input file list: " + dirlistfilename
23    sys.exit(1)
24# end try
25
26try:
27    maskfile = open(maskfilename, "w")
28except:
29    print "Could not open mask file for writing: " + maskfilename
30    sys.exit(1)
31# end try
32
33try: 
34    mergedfile = open(merged_filename, "w")
35except:
36    print "Could not open point file for writing: " + merged_filename
37    sys.exit(1)
38# end try
39
40all_bounds = list()
41
42enddir = False
43while (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
140dirfile.close()
141mergedfile.close()
142maskfile.close()