@ -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,42 +335,71 @@ 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 )
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
@ -546,3 +407,7 @@ for i in range(0, len(params_file)): #0, len(params_file)
#delete the temp files
remove_temp_files ( tmp_dir )
if __name__ == ' __main__ ' :
main ( )
