diff --git a/las_manipulation.py b/las_manipulation.py index 30537df..0865166 100644 --- a/las_manipulation.py +++ b/las_manipulation.py @@ -1,170 +1,16 @@ -# 2017088 Central Coast job - -# this is aimed at recieving a large lidar file and -#first cropping the lidar to a chosen shp (ie to a single beach) -#second try to reclassify vegetation and buildings where possible -#cropping the las file to remove the wave runup zone (and making a shapefile to crop it) --> This is a deliverable -#make a raster at a given step value and crop by the above shapefile --> this is what will be used to extract all the values and do heat maps - - -#cropping the las file to remove the wave runup zone (and making a shapefile to crop it) --> This is a deliverable - -# extract values along a given line --> This is a deliverable - - -########################### IMPORTS ########################################### +import os +import sys +import argparse import subprocess -import pandas as pd +from tqdm import tqdm import numpy as np -import shapefile from scipy import interpolate -from shapely.geometry import Point -import os -from tqdm import tqdm - -############################################################################### -########################## FIXED INPUTS ####################################### -#path to LasTools NOTE THERE CAN BE NO SPACES -path_2_lastools='C:/ProgramData/chocolatey/' -#input_file=r"J:\Project\wrl2017088 Central Coast Council Aerial Survey and Coastal Analysis\04_Working\Python\Parameter Files\Survey 2\las manipulation survey2.xlsx" -# input_file=r"J:\Project\wrl2017088 Central Coast Council Aerial Survey and Coastal Analysis\04_Working\Python\Parameter Files\Survey 2\las manipulation patonga2.xlsx" -input_file='Parameter Files/las manipulation survey2.xlsx' - -##################### UNCOMMENT THIS SECTION IF YOU WANT TO DO INDIVIDUAL BEACHES######### -##STEP ONE -#input_las1="C:\\Users\\z3331378\\Desktop\\LAStools\\XXFiles\\Wamberal2013\\wamberal2013.las" #CANNOT HAVE SPACES -#initial_crop_poly="C:\\Users\\z3331378\\Desktop\\LAStools\\XXFiles\\CC_Crop\\copacabana.shp" -#temp_las1 = "C:\\Users\\z3331378\\Desktop\\LAStools\\XXFiles\\Temp\\tmp1.las" -# -##STEP TWO -#input_las2 = temp_las1 -#temp_las2= "C:\\Users\\z3331378\\Desktop\\LAStools\\XXFiles\\Temp\\tmp2.las" -#step_veg=5 #step size to remove vege. Probably don't change -#temp_las3="C:\\Users\\z3331378\\Desktop\\LAStools\\XXFiles\\Temp\\tmp3.las" -#step_build=35 #step size to remove buildings. May need to change. bigger buildings mean bigger step etc -# -##STEP THREE -#input_las3=temp_las3 -#temp_xyz= "C:\\Users\\z3331378\\Desktop\\LAStools\\XXFiles\\Temp\\tmp_thinned.xyz" -#zone_MGA=56 #assumes it is in MGA -#check_value=1 #m AHD value above which you believe, shouldn't delete -#direct = 'west_east' #at the moment can only be north_south or west_east, wont work if it is south_north -#check_distance=1000 #maximum distance you can crop away from the desnified polygon (should probably be somewhere between 10 - 50m) -#polygon_name="copacabana_crop" -#path_2_poly="C:\\Users\\z3331378\\Desktop\\LAStools\\XXFiles\\Results\\" - -#tmp_dir="C:\\Users\\z3331378\\Desktop\\LAStools\\XXFiles\\Temp" -############################################################################### - -############################### SUB ROUTINES ################################## -def check_output(command,console): - - if console == True: - process = subprocess.Popen(command) - else: - process = subprocess.Popen(command, shell=True, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, universal_newlines=True) - output,error = process.communicate() - returncode = process.poll() - - return returncode,output - -def crop_las(las, shp, output, lastools_loc): - # output is the full path and filename (inc extension) to put in - path_2_lasclip=lastools_loc+"\\bin\\lasclip" - - command="%s -i %s -poly %s -o %s" % (path_2_lasclip, las, shp, output) - - returncode,output = check_output(command, False) - - if returncode!= 0: - print("Error. lasclip failed on %s" % shp.split('//')[-1].split('.')[0]) - else: - return None - -def colour_las(las, tif, out_las, lastools_loc): - # output is the full path and filename (inc extension) to put in - path_2_lascolor=lastools_loc+"\\bin\\lascolor" - - command="%s -i %s -image %s -o %s" % (path_2_lascolor, las, tif, out_las) - - returncode,output = check_output(command, False) - - if returncode!= 0: - print("Error. lascolor failed on") - else: - return None - -def remove_veg(las, output, lastools_loc, step=5): - # output is the full path and filename (inc extension) to put in - #step is how the features are removes - the step size should be able to distinguish between the features you are removing - #I THINK that the way lasground works is to create a tin at the step size specified and find points that don't make sense - # by comparing to the tin - #to remove veg, a step size of 3 - 5 is typically used - path_2_lasground=lastools_loc+"\\bin\\lasground_new" - - command="%s -i %s -o %s -step 5" % (path_2_lasground, las, output) - - returncode,output = check_output(command, False) - - if returncode!= 0: - print("Error. lasground failed on %s" % las.split('\\')[-1].split('.')[0]) - else: - return None - -def remove_buildings(las, output, lastools_loc, step=50): - # output is the full path and filename (inc extension) to put in - #step is how the features are removes - the step size should be able to distinguish between the features you are removing - #I THINK that the way lasground works is to create a tin at the step size specified and find points that don't make sense - # by comparing to the tin - #to remove buildings, a step size of 10 -50 is typically used - #too high a step size will miss small buildings (eg garages, SLSC huts etc), too small will miss large buildings (ie stadiums and warehouses) - # use a step size of ~25m for houses etc. - - path_2_lasground=lastools_loc+"\\bin\\lasground_new" - - command="%s -i %s -o %s -step %s" % (path_2_lasground, las, output, step) - - returncode,output = check_output(command, False) - - if returncode!= 0: - print("Error. lasground failed on %s" % las.split('\\')[-1].split('.')[0]) - else: - return None - -def make_raster(las, output, lastools_loc, keep_only_ground=False): - #not that keep ground only option only rasters points classified as "2" in the lidar (ie ground) - #this effectively creates a "bare earth dem" - #note that this should only be used after remove_veg and/or remove_buildings has been run - - path_2_blastdem=lastools_loc+"\\bin\\blast2dem" - - if keep_only_ground==False: - command="%s -i %s -o %s -step 0.2" % (path_2_blastdem, las, output) - else: - command="%s -i %s -o %s -step 0.2 -keep_class 2" % (path_2_blastdem, las, output) - - returncode,output2 = check_output(command, False) - - if returncode!= 0: - print("Error. blast2dem failed on %s" % las.split('\\')[-1].split('.')[0]) - else: - return None - -def make_xyz(las, output, lastools_loc, step=1): - #to thin the las and create regular grid to allow the computation of the wave runup polygon - #output must have location and .xyz file extension - - path_2_blastdem=lastools_loc+"\\bin\\blast2dem" - - command="%s -i %s -o %s -step %s" % (path_2_blastdem, las, output, step) - +import pandas as pd +import geopandas as gpd +from shapely.geometry import Point, Polygon - returncode,output2 = check_output(command, False) +from survey_tools import call_lastools - if returncode!= 0: - print("Error. blast2dem failed on %s" % las.split('\\')[-1].split('.')[0]) - else: - return None def remove_problems(x_list, y_list, z_list, x_now, y_now, check_value): @@ -342,8 +188,7 @@ def distance_point_to_poly(x_list, y_list, x_now, y_now): return dist -def polygon_wave_runup(xyz_1m, direction, path_2_crop_polygon, crop_poly_name, set_check_value, distance_check, zone): - +def polygon_wave_runup(xyz_1m, direction, shp_name, set_check_value, distance_check, zone): #print('starting processing of wave runup') all_data=pd.read_csv(xyz_1m, header=None, names=['X','Y','Z']) @@ -477,24 +322,11 @@ def polygon_wave_runup(xyz_1m, direction, path_2_crop_polygon, crop_poly_name, s #print('exceeded the manual distance_check %s%% of the time. manually cropping undertaken' % (round(len(exceed_list)/a,2)*100)) #making the cropping shapefile #print('making the crop polygon') - w = shapefile.Writer(shapefile.POLYGON) - w.poly(parts=[for_shape]) - fname = crop_poly_name - zone=56 - w.field('FIRST_FLD','C','40') - w.field('SECOND_FLD','C','40') - w.record('Poly','PolyTest') - w.save('%s%s' % (path_2_crop_polygon,fname)) - prjfname='%s%s.prj'%(path_2_crop_polygon,fname) - prj=open(prjfname, 'w') - if zone==56: - prj.write('PROJCS["GDA_1994_MGA_Zone_56",GEOGCS["GCS_GDA_1994",DATUM["D_GDA_1994",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",10000000.0],PARAMETER["Central_Meridian",153.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]') - elif zone==55: - prj.write('PROJCS["GDA_1994_MGA_Zone_55",GEOGCS["GCS_GDA_1994",DATUM["D_GDA_1994",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",10000000.0],PARAMETER["Central_Meridian",147.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]') - elif zone==57: - prj.write('PROJCS["GDA_1994_MGA_Zone_57",GEOGCS["GCS_GDA_1994",DATUM["D_GDA_1994",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",10000000.0],PARAMETER["Central_Meridian",159.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]') - prj.close() + # Export polygon as shapefile + df = gpd.GeoDataFrame(geometry=[Polygon(for_shape)]) + df.crs = {'init': 'epsg:283{}'.format(zone), 'no_defs': True} + df.to_file(shp_name + '.shp', driver='ESRI Shapefile') return None @@ -503,46 +335,79 @@ def remove_temp_files(directory): os.unlink(os.path.join(directory, f)) return None -############################################################################### - -############################### CODE ######################################### -# read the parameters file and scroll through it -params_file=pd.read_excel(input_file, sheet_name="PARAMS") - -for i in range(0, len(params_file)): #0, len(params_file) - print("Starting to process %s" % params_file['Beach'][i]) - input_las1=params_file['INPUT LAS1'][i] - initial_crop_poly=params_file['INITIAL CROP POLY'][i] - temp_las1=params_file['TEMP LAS1'][i] - temp_las2 = params_file['TEMP LAS2'][i] - step_veg = params_file['STEP VEG'][i] - temp_las3 = params_file['TEMP LAS3'][i] - step_build= params_file['STEP BUILD'][i] - temp_xyz=params_file['TEMP XYZ'][i] - zone_MGA=params_file['ZONE MGA'][i] - check_value=params_file['CHECK VALUE'][i] - direct=params_file['DIRECTION'][i] - check_distance=params_file['CHECK DISTANCE'][i] - polygon_name=params_file['POLYGON NAME'][i] - path_2_poly = params_file['PATH TO POLYGON'][i] - temp_las4= params_file['TEMP LAS4'][i] - picture= params_file['PICTURE'][i] - tmp_dir=params_file['TMP FOLDER'][i] - - #STEP ONE CROP TO BEACH - crop_las(input_las1,initial_crop_poly, temp_las1, path_2_lastools) - - #STEP TWO REMOVE VEG - remove_veg(temp_las1, temp_las2, path_2_lastools, step=step_veg) - remove_buildings(temp_las2, temp_las3, path_2_lastools, step=step_build) - - #STEP THREE MAKE WAVE RUNUP REMOVAL POLYGON - make_xyz(temp_las3, temp_xyz, path_2_lastools) - polygon_wave_runup(temp_xyz, direct, path_2_poly, polygon_name, check_value, check_distance, zone_MGA) - #NOTE THAT YOU NEED TO CHECK THE OUTPUT SHP FILE AND ADJUST AS REQUIRED - - #STEP FOUR COLOURISE THE LAS - #colour_las(temp_las3, picture, temp_las4, path_2_lastools) - - #delete the temp files - remove_temp_files(tmp_dir) + +def main(): + parser = argparse.ArgumentParser() + parser.add_argument( + 'input_file', + metavar='PARAMS_FILE', + help='name of parameter file', + default=None) + + # Print usage if no arguments are provided + if len(sys.argv) == 1: + parser.print_help(sys.stderr) + sys.exit(1) + + args = parser.parse_args() + + # read the parameters file and scroll through it + input_file = args.input_file + input_file = 'Parameter Files/las manipulation survey2.xlsx' + params_file=pd.read_excel(input_file, sheet_name="PARAMS") + + for i, row in params_file.iterrows(): + print("Starting to process %s" % row['Beach']) + input_las1 = row['INPUT LAS1'] + initial_crop_poly = row['INITIAL CROP POLY'] + # temp_las1 = row['TEMP LAS1'] + # temp_las2 = row['TEMP LAS2'] + step_veg = row['STEP VEG'] + # temp_las3 = row['TEMP LAS3'] + step_build = row['STEP BUILD'] + temp_xyz = row['TEMP XYZ'] + zone_MGA = row['ZONE MGA'] + check_value = row['CHECK VALUE'] + direct = row['DIRECTION'] + check_distance = row['CHECK DISTANCE'] + polygon_name = row['POLYGON NAME'] + path_2_poly = row['PATH TO POLYGON'] + # temp_las4 = row['TEMP LAS4'] + picture = row['PICTURE'] + tmp_dir = row['TMP FOLDER'] + + # Crop to beach boundary + print('Clipping...') + las_data = call_lastools('lasclip', input=input_las1, output='-stdout', + args=['-poly', initial_crop_poly], verbose=False) + + # Remove vegetation + print('Removing vegetation...') + las_data = call_lastools('lasground_new', input=las_data, output='-stdout', + args=['-step', step_veg], verbose=False) + # Remove buildings + print('Removing buildings...') + las_data = call_lastools('lasground_new', input=las_data, output='-stdout', + args=['-step', step_build], verbose=False) + + # Interpolate point cloud onto a grid + print('Interpolating to grid...') + call_lastools('las2dem', input=las_data, output=temp_xyz, + args=['-step', 1], verbose=False) + + # Make runup clipping mask + print('Calculating runup clipping mask...') + shp_name = os.path.join(path_2_poly, polygon_name) + polygon_wave_runup(temp_xyz, direct, shp_name, check_value, check_distance, zone_MGA) + + #NOTE THAT YOU NEED TO CHECK THE OUTPUT SHP FILE AND ADJUST AS REQUIRED + + #STEP FOUR COLOURISE THE LAS + #colour_las(temp_las3, picture, temp_las4, path_2_lastools) + + #delete the temp files + remove_temp_files(tmp_dir) + + +if __name__ == '__main__': + main()