Update 'las_outputs.py'

Move legend in profile plots to RHS

README.md

@@ -63,7 +63,7 @@ python las_outputs.py param-files/survey-1*.yaml param-files/survey-2*.yaml
 Example usage:
 
 ```
-C:/Python27/ArcGIS10.4/python generate_heatmaps.py param-files/survey-2*.yaml
+C:\Python27\ArcGIS10.5\python generate_heatmaps.py param-files\survey-2*.yaml
 
 ```
 
@@ -76,5 +76,29 @@ C:/Python27/ArcGIS10.4/python generate_heatmaps.py param-files/survey-2*.yaml
 Example usage:
 
 ```
-C:/Python27/ArcGIS10.4/python plot_heatmaps.py param-files/survey-2-avoca.yaml
+C:\Python27\ArcGIS10.5\python plot_heatmaps.py param-files\survey-2-avoca.yaml
+```
+
+### 7. Export 0.7 m contour for DSAS trend analysis
+`extract_contours.py` extracts a predefined contour (say 0.7 m AHD) and outputs it in a format suitable for DSAS shoreline trend analysis.
+
+Example usage:
+```
+python extract_contours.py
+```
+
+## polyline_to_points.py
+This script interpolates points along transects in a shapefile.
+
+Example usage:
+
+```
+# Extract points at default spacing (1m)
+$ python polyline_to_points.py path/to/shp
+
+# Extract points at 5m increments
+$ python polyline_to_points.py -s 5 path/to/shp
+
+# Use profile names from field "ProfileID" in the attribute table at 1 m spacing
+$ python polyline_to_points.py -f ProfileID -s 1 path/to/shp
 ```

generate_heatmaps.py

@@ -72,10 +72,10 @@ def process(yaml_name):
         raise ValueError('No previous survey date provided')
 
     # Set paths for raster files
-    current_raster = os.path.join(input_tif_dir, base_name + '.tif')
+    current_raster = os.path.join(input_tif_dir, base_name + '_DEM.tif')
     previous_base_name = re.sub('\d+', previous_date, base_name)
-    previous_raster = os.path.join(input_tif_dir, previous_base_name + '.tif')
+    previous_raster = os.path.join(input_tif_dir, previous_base_name + '_DEM.tif')
-    heatmap_raster = os.path.join(output_tif_dir, base_name + '.tif')
+    heatmap_raster = os.path.join(output_tif_dir, base_name + '_heatmap.tif')
 
     print('processing {}'.format(beach))
 

las_manipulation.py

@@ -406,7 +406,7 @@ def process(yaml_file):
     # Interpolate point cloud onto a grid
     print('Interpolating to grid...')
     xyz_name = os.path.join(tmp_dir, las_basename + '.xyz')
-    call_lastools('las2dem', input=las_classified_name, output=xyz_name,
+    call_lastools('blast2dem', input=las_classified_name, output=xyz_name,
                   args=['-step', 1], verbose=False)
 
     # Make runup clipping mask from gridded point cloud

las_outputs.py

@@ -87,14 +87,14 @@ def plot_profiles(profile_name, csv_output_dir, graph_loc, ch_limits, delta_vol,
         fig_h = profiles.dropna().values.max() / m_per_inch * vertical_exag
         fig_w = (profiles.index.max() - ch_min) / m_per_inch
     except ValueError:
-        fig_h = 5
+        fig_h = 2.3
         fig_w = 10
 
     if scale_figures:
         fig, ax = plt.subplots(figsize=(fig_w, fig_h))
         ax.set_aspect(vertical_exag)
     else:
-        fig, ax = plt.subplots(figsize=(10, 5))
+        fig, ax = plt.subplots(figsize=(10, 2.3))
 
 
     for col in profiles.columns:
@@ -129,7 +129,7 @@ def plot_profiles(profile_name, csv_output_dir, graph_loc, ch_limits, delta_vol,
                     backgroundcolor='#ffffff', linespacing=1.5)
 
     ax.legend(edgecolor='none', facecolor='#ffffff', fontsize=9, labelspacing=1)
-
+              loc='lower left', bbox_to_anchor=(1.02,0))
     ax.xaxis.set_minor_locator(MultipleLocator(5))
     ax.yaxis.set_minor_locator(MultipleLocator(0.5))
     ax.xaxis.grid(True, which='minor', color='k', linewidth=0.5, alpha=0.3)
@@ -262,7 +262,7 @@ def process(yaml_file):
 
     # Make a raster from point cloud
     print('Creating heat map raster...')
-    tif_name = os.path.join(output_tif_dir, las_basename + '.tif')
+    tif_name = os.path.join(output_tif_dir, las_basename + '_DEM.tif')
     call_lastools('las2dem', input=las_data, output=tif_name,
                   args=['-step', 1, '-keep_class', 2], verbose=False)
     # IF THIS STEP ISN'T WORKING:

plot_heatmaps.py

@@ -70,7 +70,7 @@ def process(yaml_name):
         raise ValueError('No previous survey date provided')
 
     # Set paths for raster files
-    heatmap_raster = os.path.join(input_tif_dir, base_name + '.tif')
+    heatmap_raster = os.path.join(input_tif_dir, base_name + '_heatmap.tif')
 
     print('processing {}'.format(beach))
     mxd = arcpy.mapping.MapDocument(mxd_name)

polyline_to_points.py

@@ -0,0 +1,92 @@
+"""polyline_to_points.py
+Extract interpolated points along transects in a shapefile.
+
+D. Howe
+d.howe@wrl.unsw.edu.au
+2020-02-19
+"""
+import sys
+import argparse
+import numpy as np
+import pandas as pd
+import geopandas as gpd
+from shapely.geometry import LineString
+
+
+def extract(shp_path, spacing=1, field=None):
+
+    rows = []
+    shp = gpd.read_file(shp_path)
+
+    if field is None:
+        # Assume profile name is second field in shapefile
+        field = shp.columns[1]
+
+    for i, line in shp.iterrows():
+        g = line.geometry
+        chainages = np.arange(0, g.length, step=spacing)
+        for chainage in chainages:
+            easting, northing = g.interpolate(chainage).xy
+
+            row = {
+                'ProfileNum': line[field],
+                'Easting': easting[0],
+                'Northing': northing[0],
+                'Chainage': chainage,
+            }
+            rows.append(row)
+
+    # Create output dataframe
+    df = pd.DataFrame(rows)
+
+    # Re-order columns
+    df = df[['ProfileNum', 'Easting', 'Northing', 'Chainage']]
+
+    # Export to csv
+    csv_path = shp_path.replace('.shp', '.csv')
+    df.to_csv(csv_path, index=False)
+
+
+def main():
+    example_text = """examples:
+    # Extract points at default spacing (1m)
+    $ python polyline_to_points.py path/to/shp
+
+    # Extract points at 5m increments
+    $ python polyline_to_points.py -s 5 path/to/shp
+
+    # Use profile names from field "ProfileID" in the attribute table
+    $ python polyline_to_points.py -f ProfileID path/to/shp
+    """
+    # Set up command line arguments
+    parser = argparse.ArgumentParser(
+        epilog=example_text,
+        formatter_class=argparse.RawDescriptionHelpFormatter)
+    parser.add_argument('shp_path',
+                        metavar='SHAPEFILE',
+                        help='path to input shapefile')
+    parser.add_argument('-s',
+                        '--spacing',
+                        metavar='SPACING',
+                        default=1,
+                        type=int,
+                        help='space between points (default=1)')
+    parser.add_argument('-f',
+                        '--field',
+                        metavar='FIELDNAME',
+                        help='profile field name in attribute table')
+
+    # 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()
+
+    # Combine images
+    extract(**vars(args))
+
+
+if __name__ == '__main__':
+    main()