Processing/LIDARDEMs: generate_mask.py

import bz2
import re
import math
import sys

NORTH_COL = 2
EAST_COL = 1
PADDING = 2000

if (len(sys.argv) != 4):
    print "Usage: python generate_mask.py file_list mask_file point_file"
    sys.exit(1)
# end if

dirlistfilename = sys.argv[1]
maskfilename = sys.argv[2]
merged_filename = sys.argv[3]

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

try:
    maskfile = open(maskfilename, "w")
except:
    print "Could not open mask file for writing: " + maskfilename
    sys.exit(1)
# end try

try:  
    mergedfile = open(merged_filename, "w")
except:
    print "Could not open point file for writing: " + merged_filename
    sys.exit(1)
# end try

all_bounds = list()

enddir = False
while (not enddir):
    north_max = -1
    south_max = -1
    east_max = -1
    west_max = -1
    nw_max = -1
    sw_max = -1
    ne_max = -1
    se_max = -1
    linebounds = dict()
        
    datafilename = dirfile.readline()
    if (datafilename == ""):
        enddir = True
        break
    else:
        datafilename = datafilename.strip() # Strip newline character off end of line if necessary
    # end if
    
    try:
        datafile = bz2.BZ2File(datafilename)
    except:
        print "Could not open lidar input file: " + datafilename
        sys.exit(1)
    # end try
        
    for currentline in datafile:
        
        # Strip off zone number (17th-18th characters on each line)
        strippedline = ""
        for i in range(0, len(currentline)):
            if ((i != 16) and (i != 17)):
                strippedline = strippedline + currentline[i]
            # end if
        # end for
        
        strippedline = strippedline.strip()
        pointfields = re.split("\s+", strippedline) # Split line on whitespace
        point_north = float(pointfields[NORTH_COL])
        point_east = float(pointfields[EAST_COL])
        
        if ((north_max == -1) or (north_max < int(math.ceil(point_north)) + PADDING)):
            north_max = int(math.ceil(point_north)) + PADDING
            nmax_point = (int(math.floor(point_east)), north_max)
        # end if
        
        if ((south_max == -1) or (south_max > int(math.floor(point_north)) - PADDING)):
            south_max = int(math.floor(point_north)) - PADDING
            smax_point = (int(math.floor(point_east)), south_max)
        # end if
        
        if ((east_max == -1) or (east_max < int(math.ceil(point_east)) + PADDING)):
            east_max = int(math.ceil(point_east)) + PADDING
            emax_point = (east_max, int(math.floor(point_north)))
        # end if
        
        if ((west_max == -1) or (west_max > int(math.floor(point_east)) - PADDING)):
            west_max = int(math.floor(point_east)) - PADDING
            wmax_point = (west_max, int(math.floor(point_north)))
        # end if
        
        # Max north-west point defined by 2*(ycoord-xcoord) (derived from geometry against SW/NE axis)
        if ((nw_max == -1) or (nw_max < int(math.ceil(2*(point_north - point_east))))):
            nw_max = int(math.ceil(2*(point_north - point_east)))
            nwmax_point = (int(math.floor(point_east)-PADDING), int(math.ceil(point_north)+PADDING))
        # end if
        
        if ((se_max == -1) or (se_max > int(math.floor(2*(point_north - point_east))))):
            se_max = int(math.floor(2*(point_north - point_east)))
            semax_point = (int(math.floor(point_east)-PADDING), int(math.floor(point_north)-PADDING))
        # end if
        
        # Max north-east point defined by 2x(ycoord+xcoord) (derived from geometry against SE/NW axis)
        if ((ne_max == -1) or (ne_max < int(math.ceil(2*(point_north + point_east))))):
            ne_max = int(math.ceil(2*(point_north + point_east)))
            nemax_point = (int(math.ceil(point_east)+PADDING), int(math.ceil(point_north)+PADDING))
        # end if
        
        if ((sw_max == -1) or (sw_max > int(math.floor(2*(point_north + point_east))))):
            sw_max = int(math.floor(2*(point_north + point_east)))
            swmax_point = (int(math.ceil(point_east)+PADDING), int(math.floor(point_north)-PADDING))
        # end if
        
        mergedfile.write(pointfields[0] + "  " + pointfields[1] + "  " + pointfields[2] + "  " + pointfields[3] + "  " + pointfields[4] + "\n")
    # end for
    
    # Write points to mask file
    maskfile.write(str(nmax_point[0]) + "," + str(nmax_point[1]) + "\n")
    maskfile.write(str(nwmax_point[0]) + "," + str(nwmax_point[1]) + "\n")
    maskfile.write(str(wmax_point[0]) + "," + str(wmax_point[1]) + "\n")
    maskfile.write(str(swmax_point[0]) + "," + str(swmax_point[1]) + "\n")
    maskfile.write(str(smax_point[0]) + "," + str(smax_point[1]) + "\n")
    maskfile.write(str(semax_point[0]) + "," + str(semax_point[1]) + "\n")
    maskfile.write(str(emax_point[0]) + "," + str(emax_point[1]) + "\n")
    maskfile.write(str(nemax_point[0]) + "," + str(nemax_point[1]) + "\n")
# end while

dirfile.close()
mergedfile.close()
maskfile.close()