From fd5ddea7f22e28c95679bf4cdd9eacb1c44454b1 Mon Sep 17 00:00:00 2001
From: Dan Howe <d.howe@wrl.unsw.edu.au>
Date: Mon, 2 Jul 2018 15:45:43 +1000
Subject: [PATCH] Add las_manipulation.py

---
 las_manipulation.py | 548 ++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 548 insertions(+)
 create mode 100644 las_manipulation.py

diff --git a/las_manipulation.py b/las_manipulation.py
new file mode 100644
index 0000000..30537df
--- /dev/null
+++ b/las_manipulation.py
@@ -0,0 +1,548 @@
+# 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 subprocess
+import pandas as pd
+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)
+
+
+    returncode,output2 = check_output(command, False)
+
+    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):
+
+    z_ave=nine_pt_moving_average(z_list)
+    deriv_ave, chainage=forward_derivative(x_list, y_list, z_ave)
+    deriv_real, chainage=two_point_derivative(x_list, y_list, z_list)
+
+    #first find the reference contour on the beach
+    #index_contour, x_now, y_now, distance=find_beach_reference_contour_choose_closest(chainage, z_ave, x_list, y_list, x_now, y_now,deriv_ave, check_value)
+    index_contour, x_now, y_now, distance=find_beach_reference_contour(chainage, z_ave, x_list, y_list, x_now, y_now,deriv_ave,deriv_real,check_value)
+    if index_contour<len(chainage): #other wise keep everthing
+        #find the beach slope, get the interpolated line (beach line) and the index of the reference contour +1
+        beach_slope, beach_line, index_high=find_beach_slope(chainage, z_ave,index_contour, check_value)
+        #find the natural deviation of the lower beach
+        nat_dev=get_natural_deviation(chainage, z_list, index_contour, index_high, beach_line)
+
+        for i in range(index_contour, len(z_list)):
+            if abs(z_list[i]-float(beach_line(chainage[i])))>nat_dev:
+                break
+    else:
+        i=index_contour
+
+    z_return=z_list[0:i]
+    chainage_return=chainage[0:i]
+    return z_return, chainage_return, x_now, y_now, distance
+
+def two_point_derivative(x_list, y_list, z_list):
+    chain=[((x_list[0]-x_list[i])**2+(y_list[0]-y_list[i])**2)**0.5 for i in range(0,len(x_list))]
+    deriv=[(z_list[i+1]-z_list[i-1])/(chain[i+1]-chain[i-1]) for i in range(1, len(z_list)-1)]
+
+    deriv.insert(0,0)
+    deriv.append(0)
+
+    return deriv, chain
+
+def forward_derivative(x_list, y_list, z_list):
+    chain=[((x_list[0]-x_list[i])**2+(y_list[0]-y_list[i])**2)**0.5 for i in range(0,len(x_list))]
+    deriv=[(z_list[i]-z_list[i-1])/(chain[i]-chain[i-1]) for i in range(0, len(z_list)-1)]
+
+    deriv.insert(0,0)
+
+    return deriv, chain
+
+def find_first_over_reference(z_list, value):
+    i=len(z_list)-1
+
+    while i>0 and z_list[i]<value:
+        i=i-1
+
+    return i
+
+def nine_pt_moving_average(z_list):
+    i=0
+    move_ave=[]
+    while i<len(z_list):
+        if i<5:
+            ave=np.mean([z_list[j] for j in range(0,i+5)])
+        elif i>len(z_list)-5:
+            ave=np.mean([z_list[j] for j in range(i-4,len(z_list))])
+        else:
+            ave=np.mean([z_list[j] for j in range(i-4,i+5)])
+
+        move_ave.append(ave)
+        i=i+1
+    return move_ave
+
+def find_neg_derivative(z_list, deriv_list):
+    i=len(z_list)-5
+    while z_list[i]>=0 and z_list[i+1]>=0 and z_list[i+2]>=0 and z_list[i+3]>=0 and z_list[i+4]>=0:
+        i=i-1
+
+    return i
+
+def find_beach_reference_contour_choose_closest(chain_list, z_ave_list, x_list, y_list, x_last, y_last, deriv_ave_list, check_value):
+    #note that z_list should be the 9 point moving average
+    #assumes that beaches are shallow (|derivative|<0.3), sloping and between 0 - 4 m AHD
+    i=0
+    choice_list=[]
+    distance_list=[]
+    if z_ave_list[i]>check_value:
+        state_now='over'
+    else:
+        state_now='under'
+
+    while i<len(z_ave_list):
+        if state_now=='under' and z_ave_list[i]>check_value: #only keep if it is downward sloping
+            state_now='over'
+        elif state_now=='over' and z_ave_list[i]<check_value:
+            choice_list.append(i)
+            state_now='under'
+            if x_last!=None:
+                distance_list.append(((x_last - x_list[i])**2+(y_last - y_list[i])**2)**0.5)
+        i=i+1
+
+    if len(choice_list)>0 and x_last==None: #choose the first time for the first point
+        i=choice_list[0]
+        distance=0
+    elif len(choice_list)>0 and x_last!=None:
+        assert(len(choice_list)==len(distance_list))
+        i=choice_list[distance_list.index(min(distance_list))]
+        distance=min(distance_list)
+
+
+    if i>=len(x_list):
+        i=len(x_list)-1
+        if x_last!=None:
+            distance=((x_last - x_list[i])**2+(y_last - y_list[i])**2)**0.5
+        else:
+            distance=0
+
+    x=x_list[i]
+    y=y_list[i]
+
+    return i, x, y, distance
+
+def find_beach_reference_contour(chain_list, z_ave_list, x_list, y_list, x_last, y_last, deriv_ave_list,deriv_real_list, check_value):
+    #note that z_list should be the 9 point moving average
+    #assumes that beaches are shallow (|derivative|<0.3), sloping and between 0 - 4 m AHD
+
+    i=len(z_ave_list)-1
+    while i>=0 and (z_ave_list[i]>check_value+2 or z_ave_list[i]<check_value-2  or deriv_ave_list[i]>0 or max([abs(i) for i in deriv_real_list[max(0,i-7):i]]+[0])>0.3):#beaches are shallow sloping, low
+        i=i-1
+
+    #find the first time it gets to check_value after this
+    while i>=0 and z_ave_list[i]<check_value:
+        i=i-1
+
+    if i==0:
+        i=len(z_ave_list)-1 # the whole this is above the beach
+
+    if x_last!=None:
+        distance=((x_last - x_list[i])**2+(y_last - y_list[i])**2)**0.5
+    else:
+        distance=0
+    x=x_list[i]
+    y=y_list[i]
+
+    return i, x, y, distance
+
+
+def find_beach_slope(chain_list, z_ave_list, ref_index, check_value):
+    #ref index is the index of the check value
+    #find the beach slope between this point and 1 m above this point
+
+    i=ref_index
+    while i>0 and z_ave_list[i]<check_value+1:
+        i=i-1
+
+    slope=(z_ave_list[i]-z_ave_list[ref_index])/(chain_list[i]-chain_list[ref_index])
+
+    beach_ave=interpolate.interp1d([min(chain_list),max(chain_list)], [(min(chain_list)-chain_list[ref_index])*slope+z_ave_list[ref_index], (z_ave_list[ref_index]-(chain_list[ref_index]-max(chain_list))*slope)])
+
+    return slope, beach_ave, i
+
+def get_natural_deviation(chain_list, z_list, ref_index, ref_high, beach_ave):
+    #for the points considered to be on the beach (reference contour to reference contour +1), find the average natural deviation
+    deviation=[]
+    for i in range(ref_high, ref_index+1):
+        dev_tmp=abs(z_list[i] - float(beach_ave(chain_list[i])))
+        deviation.append(dev_tmp)
+
+    natural_deviation=min(np.max(deviation),0.4) #THIS MAY BE TOO CONSERVATIVE
+
+    return natural_deviation
+
+def distance_point_to_poly(x_list, y_list, x_now, y_now):
+    #make a line from the mid of x_list, y_list
+    end=Point(x_list[-1], y_list[-1])
+
+    point=Point(x_now, y_now)
+
+
+    dist=point.distance(end)
+
+
+    return dist
+
+def polygon_wave_runup(xyz_1m, direction, path_2_crop_polygon, crop_poly_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'])
+
+    if direction=='north_south':
+        all_data_sorted=all_data.sort_values(by=['X', 'Y'], ascending=[1,0])
+    elif direction=='west_east':
+        all_data_sorted=all_data.sort_values(by=['Y', 'X'], ascending=[0,1])
+
+    fixed_now=0
+    a=0
+    X_tmp=[]
+    processed_data = pd.DataFrame(columns=['X','Y','Z'])
+    list_to_print=[10,20,30,40,50,60,70,80,90]
+    crop_line=[]
+    top_line=[]
+    tmp_x_last=None
+    tmp_y_last=None
+    exceed_list=[]
+
+    # Create progress bar
+    pbar = tqdm(all_data_sorted.iterrows(), total=all_data_sorted.shape[0])
+    for index, line in pbar:
+        a=a+1
+        percent_done=round(a/len(all_data_sorted)*100,1)
+        if percent_done in list_to_print:
+            #print("Finished %s%% of the processing" % percent_done)
+            list_to_print=list_to_print[1:len(list_to_print)]
+        if direction=='north_south':
+            check_this=line['X']
+        elif direction=='west_east':
+            check_this=line['Y']
+        if check_this==fixed_now:
+            X_tmp.append(line['X'])
+            Y_tmp.append(line['Y'])
+            Z_tmp.append(line['Z'])
+        else:
+            if len(X_tmp)!=0:
+                #try: ########may need to change!~!
+                if len(X_tmp)>10:
+                    Z_set, chainage_tmp, temp_x, temp_y, distance=remove_problems(X_tmp, Y_tmp, Z_tmp,tmp_x_last, tmp_y_last, set_check_value)
+                #except:
+                else:
+                    Z_set=Z_tmp
+                    temp_x=X_tmp[len(Z_set)-1]
+                    temp_y=Y_tmp[len(Z_set)-1]
+                    distance=0
+
+                distance_2_old=distance_point_to_poly(X_tmp, Y_tmp, temp_x, temp_y)
+                if distance_2_old<distance_check: # find a way to change so it is checking the distance from the first crop polyogn, concave_now.buffer(buffer)
+                    tmp_x_last=temp_x
+                    tmp_y_last=temp_y
+                    crop_line.append([X_tmp[len(Z_set)-1], Y_tmp[len(Z_set)-1]])
+                    top_line.append([X_tmp[0], Y_tmp[0]])
+
+                    #otherwise crop by the distance_check
+                else:
+                    exceed_list.append(1)
+                    try:
+                        tmp_x_last=X_tmp[len(X_tmp)-distance_check] #beacuse this is a 1m DSM, this works
+                        tmp_y_last=Y_tmp[len(Y_tmp)-distance_check]
+
+                        crop_line.append([tmp_x_last, tmp_y_last])
+                        top_line.append([X_tmp[0], Y_tmp[0]])
+                    except:
+                        print('problem with the last crop point, keeping whole line')
+                        crop_line.append([X_tmp[-1], Y_tmp[-1]])
+                        top_line.append([X_tmp[0], Y_tmp[0]])
+
+                if direction=='north_south':
+                    fixed_now=line['X']
+                elif direction=='west_east':
+                    fixed_now=line['Y']
+                X_tmp=[line['X']]
+                Y_tmp=[line['Y']]
+                Z_tmp=[line['Z']]
+
+            else:
+                if direction=='north_south':
+                    fixed_now=line['X']
+                elif direction=='west_east':
+                    fixed_now=line['Y']
+                X_tmp=[line['X']]
+                Y_tmp=[line['Y']]
+                Z_tmp=[line['Z']]
+
+
+    #for the last line
+    derivative, chainage=forward_derivative(X_tmp, Y_tmp, Z_tmp)
+    if len(X_tmp)>10:
+        Z_set, chainage_tmp, temp_x, temp_y, distance=remove_problems(X_tmp, Y_tmp, Z_tmp,tmp_x_last, tmp_y_last, set_check_value)
+    #except:
+    else:
+        Z_set=Z_tmp
+        temp_x=X_tmp[len(Z_set)-1]
+        temp_y=Y_tmp[len(Z_set)-1]
+        distance=0
+    X_set=X_tmp[0:len(Z_set)]
+    Y_set=Y_tmp[0:len(Z_set)]
+
+
+    #write to new data frame
+    #if len(Z_set)>0:
+    #    for i in range(0, len(Z_set)):
+    #        processed_data =processed_data.append({'X':X_set[i],'Y':Y_set[i],'Z':Z_set[i],'r':r_set[i],'g':g_set[i],'b':b_set[i]}, ignore_index=True)
+    #add to crop line
+    distance_2_old=distance_point_to_poly(X_tmp, Y_tmp, temp_x, temp_y)
+    if distance_2_old<distance_check: # find a way to change so it is checking the distance from the first crop polyogn, concave_now.buffer(buffer)
+        tmp_x_last=temp_x
+        tmp_y_last=temp_y
+        crop_line.append([X_tmp[len(Z_set)-1], Y_tmp[len(Z_set)-1]])
+        top_line.append([X_tmp[0], Y_tmp[0]])
+        #otherwise crop by the distance_check
+    else:
+        exceed_list.append(1)
+        tmp_x_last=X_tmp[len(X_tmp)-distance_check]
+        tmp_y_last=Y_tmp[len(Y_tmp)-distance_check]
+        crop_line.append(tmp_x_last, tmp_y_last)
+        top_line.append([X_tmp[0], Y_tmp[0]])
+
+        #otherwise dont add.  straight line is better
+    if direction=='north_south':
+        y_filtered=nine_pt_moving_average([i[1] for i in crop_line])
+        crop_new=[[crop_line[i][0],y_filtered[i]] for i in range(0, len(crop_line))]
+    elif direction=='west_east':
+        x_filtered=nine_pt_moving_average([i[0] for i in crop_line])
+        crop_new=[[x_filtered[i],crop_line[i][1]] for i in range(0, len(crop_line))]
+
+    for_shape=crop_new+top_line[::-1]
+    for_shape.append(crop_new[0])
+    #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()
+
+
+    return None
+
+def remove_temp_files(directory):
+    for f in os.listdir(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)