|
|
|
@ -15,11 +15,13 @@ import io
|
|
|
|
|
import re
|
|
|
|
|
import sys
|
|
|
|
|
import math
|
|
|
|
|
import yaml
|
|
|
|
|
import argparse
|
|
|
|
|
import datetime
|
|
|
|
|
import subprocess
|
|
|
|
|
import numpy as np
|
|
|
|
|
import pandas as pd
|
|
|
|
|
from glob import glob
|
|
|
|
|
from cycler import cycler
|
|
|
|
|
import matplotlib.pyplot as plt
|
|
|
|
|
from matplotlib.ticker import MultipleLocator
|
|
|
|
@ -166,102 +168,119 @@ def calculate_volumes(profile_name, survey_date, csv_output_dir, ch_limits, volu
|
|
|
|
|
return delta_vol
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
def process(yaml_file):
|
|
|
|
|
with open(yaml_file, 'r') as f:
|
|
|
|
|
params = yaml.safe_load(f)
|
|
|
|
|
|
|
|
|
|
print("Starting to process %s" % params['BEACH'])
|
|
|
|
|
beach = params['BEACH']
|
|
|
|
|
survey_date = str(params['SURVEY DATE'])
|
|
|
|
|
original_las = params['INPUT LAS']
|
|
|
|
|
classified_las_dir = params['LAS CLASSIFIED FOLDER']
|
|
|
|
|
shp_swash_dir = params['SHP SWASH FOLDER']
|
|
|
|
|
crop_heatmap_poly = params['HEATMAP CROP POLY']
|
|
|
|
|
output_las_dir = params['LAS OUTPUT FOLDER']
|
|
|
|
|
zone_MGA = params['ZONE MGA']
|
|
|
|
|
output_poly_dir = params['SHP RASTER FOLDER']
|
|
|
|
|
output_tif_dir = params['TIF OUTPUT FOLDER']
|
|
|
|
|
cp_csv = params['INPUT CSV']
|
|
|
|
|
profile_limit_file = params['PROFILE LIMIT FILE']
|
|
|
|
|
csv_output_dir = params['CSV OUTPUT FOLDER']
|
|
|
|
|
graph_loc = params['PNG OUTPUT FOLDER']
|
|
|
|
|
volume_output_dir = params['CSV VOLUMES FOLDER']
|
|
|
|
|
tmp_dir = params['TMP FOLDER']
|
|
|
|
|
|
|
|
|
|
# Get base name of input las
|
|
|
|
|
las_basename = os.path.splitext(os.path.basename(original_las))[0]
|
|
|
|
|
|
|
|
|
|
# Get name of input point cloud
|
|
|
|
|
input_las = os.path.join(classified_las_dir, las_basename + '.las')
|
|
|
|
|
|
|
|
|
|
# Get name of swash cropping polygon
|
|
|
|
|
crop_swash_poly = os.path.join(shp_swash_dir, las_basename + '.shp')
|
|
|
|
|
|
|
|
|
|
# Crop point cloud to swash boundary
|
|
|
|
|
print('Cropping swash...')
|
|
|
|
|
las_data = call_lastools('lasclip', input=input_las, output='-stdout',
|
|
|
|
|
args=['-poly', crop_swash_poly], verbose=False)
|
|
|
|
|
|
|
|
|
|
# Apply sea-side clipping polygon
|
|
|
|
|
print('Cropping back of beach...')
|
|
|
|
|
las_data = call_lastools('lasclip', input=las_data, output='-stdout',
|
|
|
|
|
args=['-poly', crop_heatmap_poly], verbose=False)
|
|
|
|
|
|
|
|
|
|
# Create clipping polygon for heatmap raster
|
|
|
|
|
print('Creating heat map cropping polygon...')
|
|
|
|
|
shp_name = os.path.join(output_poly_dir, las_basename + '.shp')
|
|
|
|
|
call_lastools('lasboundary', input=las_data, output=shp_name, verbose=False)
|
|
|
|
|
|
|
|
|
|
# Make a raster from point cloud
|
|
|
|
|
print('Creating heat map raster...')
|
|
|
|
|
tif_name = os.path.join(output_tif_dir, las_basename + '.tif')
|
|
|
|
|
call_lastools('las2dem', input=las_data, output=tif_name,
|
|
|
|
|
args=['-step', 0.2, '-keep_class', 2], verbose=False)
|
|
|
|
|
|
|
|
|
|
# Extract elevations along profiles from triangulated surface
|
|
|
|
|
print('Extracting profile elevations...')
|
|
|
|
|
df = extract_pts(
|
|
|
|
|
las_data,
|
|
|
|
|
cp_csv,
|
|
|
|
|
survey_date,
|
|
|
|
|
beach,
|
|
|
|
|
args=['-parse', 'sxyz', '-keep_class', '2'],
|
|
|
|
|
verbose=False)
|
|
|
|
|
|
|
|
|
|
# Update survey profiles
|
|
|
|
|
update_survey_output(df, csv_output_dir)
|
|
|
|
|
|
|
|
|
|
# Get landward limit of surveys
|
|
|
|
|
ch_limits = pd.read_excel(profile_limit_file, index_col='Profile')
|
|
|
|
|
|
|
|
|
|
# Plot profiles, and save sand volumes for current beach
|
|
|
|
|
print('Updating figures...')
|
|
|
|
|
profile_names = df['Profile'].unique()
|
|
|
|
|
for profile_name in profile_names:
|
|
|
|
|
delta_vol = calculate_volumes(profile_name, survey_date,
|
|
|
|
|
csv_output_dir, ch_limits,
|
|
|
|
|
volume_output_dir)
|
|
|
|
|
plot_profiles(profile_name, csv_output_dir, graph_loc, ch_limits,
|
|
|
|
|
delta_vol, survey_date)
|
|
|
|
|
|
|
|
|
|
# Remove temprary 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)
|
|
|
|
|
example_text = """examples:
|
|
|
|
|
|
|
|
|
|
# Process single survey at specific beach
|
|
|
|
|
python las_outputs.py survey-1-avoca.yaml
|
|
|
|
|
|
|
|
|
|
# Process all surveys at specific beach
|
|
|
|
|
python las_outputs.py *avoca.yaml
|
|
|
|
|
|
|
|
|
|
# Process all beaches for specific survey date
|
|
|
|
|
python las_outputs.py survey-1*.yaml
|
|
|
|
|
"""
|
|
|
|
|
|
|
|
|
|
# Set up command line arguments
|
|
|
|
|
parser = argparse.ArgumentParser(
|
|
|
|
|
epilog=example_text,
|
|
|
|
|
formatter_class=argparse.RawDescriptionHelpFormatter)
|
|
|
|
|
parser.add_argument('input', help='path to yaml file(s)', nargs='*')
|
|
|
|
|
|
|
|
|
|
# Print usage if no arguments are provided
|
|
|
|
|
if len(sys.argv) == 1:
|
|
|
|
|
parser.print_help(sys.stderr)
|
|
|
|
|
sys.exit(1)
|
|
|
|
|
|
|
|
|
|
# Parse arguments
|
|
|
|
|
args = parser.parse_args()
|
|
|
|
|
yaml_files = []
|
|
|
|
|
[yaml_files.extend(glob(f)) for f in args.input]
|
|
|
|
|
|
|
|
|
|
# read the parameters file and scroll through it
|
|
|
|
|
input_file = args.input_file
|
|
|
|
|
params_file=pd.read_excel(input_file, sheet_name="PARAMS")
|
|
|
|
|
|
|
|
|
|
for i, row in params_file.iterrows():
|
|
|
|
|
print("Starting to process %s" % row['BEACH'])
|
|
|
|
|
beach=row['BEACH']
|
|
|
|
|
survey_date = str(row['SURVEY DATE'])
|
|
|
|
|
original_las = row['INPUT LAS']
|
|
|
|
|
classified_las_dir = row['LAS CLASSIFIED FOLDER']
|
|
|
|
|
shp_swash_dir = row['SHP SWASH FOLDER']
|
|
|
|
|
crop_heatmap_poly = row['HEATMAP CROP POLY']
|
|
|
|
|
output_las_dir = row['LAS OUTPUT FOLDER']
|
|
|
|
|
zone_MGA = row['ZONE MGA']
|
|
|
|
|
output_poly_dir = row['SHP RASTER FOLDER']
|
|
|
|
|
output_tif_dir = row['TIF OUTPUT FOLDER']
|
|
|
|
|
cp_csv = row['INPUT CSV']
|
|
|
|
|
profile_limit_file = row['PROFILE LIMIT FILE']
|
|
|
|
|
csv_output_dir = row['CSV OUTPUT FOLDER']
|
|
|
|
|
graph_loc = row['PNG OUTPUT FOLDER']
|
|
|
|
|
volume_output_dir = row['CSV VOLUMES FOLDER']
|
|
|
|
|
tmp_dir = row['TMP FOLDER']
|
|
|
|
|
|
|
|
|
|
# Get base name of input las
|
|
|
|
|
las_basename = os.path.splitext(os.path.basename(original_las))[0]
|
|
|
|
|
|
|
|
|
|
# Get name of input point cloud
|
|
|
|
|
input_las = os.path.join(classified_las_dir, las_basename + '.las')
|
|
|
|
|
|
|
|
|
|
# Get name of swash cropping polygon
|
|
|
|
|
crop_swash_poly = os.path.join(shp_swash_dir, las_basename + '.shp')
|
|
|
|
|
|
|
|
|
|
# Crop point cloud to swash boundary
|
|
|
|
|
print('Cropping swash...')
|
|
|
|
|
las_data = call_lastools('lasclip', input=input_las, output='-stdout',
|
|
|
|
|
args=['-poly', crop_swash_poly], verbose=False)
|
|
|
|
|
|
|
|
|
|
# Apply sea-side clipping polygon
|
|
|
|
|
print('Cropping back of beach...')
|
|
|
|
|
las_data = call_lastools('lasclip', input=las_data, output='-stdout',
|
|
|
|
|
args=['-poly', crop_heatmap_poly], verbose=False)
|
|
|
|
|
|
|
|
|
|
# Create clipping polygon for heatmap raster
|
|
|
|
|
print('Creating heat map cropping polygon...')
|
|
|
|
|
shp_name = os.path.join(output_poly_dir, las_basename + '.shp')
|
|
|
|
|
call_lastools('lasboundary', input=las_data, output=shp_name, verbose=False)
|
|
|
|
|
|
|
|
|
|
# Make a raster from point cloud
|
|
|
|
|
print('Creating heat map raster...')
|
|
|
|
|
tif_name = os.path.join(output_tif_dir, las_basename + '.tif')
|
|
|
|
|
call_lastools('las2dem', input=las_data, output=tif_name,
|
|
|
|
|
args=['-step', 0.2, '-keep_class', 2], verbose=False)
|
|
|
|
|
|
|
|
|
|
# Extract elevations along profiles from triangulated surface
|
|
|
|
|
print('Extracting profile elevations...')
|
|
|
|
|
df = extract_pts(
|
|
|
|
|
las_data,
|
|
|
|
|
cp_csv,
|
|
|
|
|
survey_date,
|
|
|
|
|
beach,
|
|
|
|
|
args=['-parse', 'sxyz', '-keep_class', '2'],
|
|
|
|
|
verbose=False)
|
|
|
|
|
|
|
|
|
|
# Update survey profiles
|
|
|
|
|
update_survey_output(df, csv_output_dir)
|
|
|
|
|
|
|
|
|
|
# Get landward limit of surveys
|
|
|
|
|
ch_limits = pd.read_excel(profile_limit_file, index_col='Profile')
|
|
|
|
|
|
|
|
|
|
# Plot profiles, and save sand volumes for current beach
|
|
|
|
|
print('Updating figures...')
|
|
|
|
|
profile_names = df['Profile'].unique()
|
|
|
|
|
for profile_name in profile_names:
|
|
|
|
|
delta_vol = calculate_volumes(profile_name, survey_date,
|
|
|
|
|
csv_output_dir, ch_limits,
|
|
|
|
|
volume_output_dir)
|
|
|
|
|
plot_profiles(profile_name, csv_output_dir, graph_loc, ch_limits,
|
|
|
|
|
delta_vol, survey_date)
|
|
|
|
|
|
|
|
|
|
# Remove temprary files
|
|
|
|
|
remove_temp_files(tmp_dir)
|
|
|
|
|
for yaml_file in yaml_files:
|
|
|
|
|
process(yaml_file)
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
if __name__ == '__main__':
|
|
|
|
|