diff --git a/extract_contours.py b/extract_contours.py
new file mode 100644
index 0000000..d69128c
--- /dev/null
+++ b/extract_contours.py
@@ -0,0 +1,56 @@
+import os
+import re
+import fiona
+from tqdm import tqdm
+import numpy as np
+import pandas as pd
+import geopandas as gpd
+from shapely.geometry import MultiLineString
+from survey_tools import call_lastools
+
+mga56 = fiona.crs.from_epsg(28356)
+
+input_dir = 'output/01-las-classified'
+output_dir = 'output/11-shp-shorelines'
+
+contour = 0.7
+smoothing = 2
+clean = 5  # any lines shorted that this are automatically removed
+
+las_names = os.listdir(input_dir)
+for las_name in tqdm(las_names):
+    shp_name = las_name.replace('.las', '_0_7_contour.shp')
+    call_lastools(
+        'blast2iso',
+        input=os.path.join(input_dir, las_name),
+        output=os.path.join(output_dir, shp_name),
+        args=['-iso', contour, '-smooth', smoothing, '-clean', clean])
+
+# Consolidate all surveys for each beach
+shp_names = [f for f in os.listdir(output_dir) if f.endswith('contour.shp')]
+
+master = gpd.GeoDataFrame({'shp_name': shp_names})
+master.geometry = [[] for s in shp_names]
+
+# Get date information
+master['beach'] = master['shp_name'].str.extract('(.*)_\d{8}')
+yyyy = master['shp_name'].str.extract('(\d{4})')
+mm = master['shp_name'].str.extract('\d{4}(\d{2})')
+dd = master['shp_name'].str.extract('\d{6}(\d{2})')
+
+# Use US date format (required by DSAS)
+master['date'] = mm + '/' + dd + '/' + yyyy
+
+# Read shapefiles
+for i, row in master.iterrows():
+    geom = gpd.read_file(os.path.join(output_dir, row['shp_name'])).geometry
+    # Combine all lines into mulit line string
+    row.geometry = MultiLineString(list(geom))
+
+    # Update fields
+    master.loc[i, 'uncertainty'] = 0
+
+# Export summary
+master.crs = {'init': 'epsg:28356'}
+master.to_file(os.path.join(output_dir, 'all_beaches.shp'),
+               driver='ESRI Shapefile')