From 061805791aaa2f882e988457ab6d3f7b42871fa6 Mon Sep 17 00:00:00 2001
From: Dan Howe <d.howe@wrl.unsw.edu.au>
Date: Mon, 26 Aug 2019 15:52:31 +1000
Subject: [PATCH] Initial commit

---
 extract_contours.py | 73 +++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 73 insertions(+)
 create mode 100644 extract_contours.py

diff --git a/extract_contours.py b/extract_contours.py
new file mode 100644
index 0000000..0602d6b
--- /dev/null
+++ b/extract_contours.py
@@ -0,0 +1,73 @@
+import os
+import re
+import fiona
+import numpy as np
+import pandas as pd
+import geopandas as gpd
+import matplotlib.pyplot as plt
+from shapely.geometry import LineString
+
+# Set contour elevation
+contour_z = 0.7
+
+# Set coordinate reference system
+epsg = 28356
+
+csv_name = 'photogrammetry_Xsections_Wamberal and Terrigal.csv'
+
+# Load csv
+df = pd.read_csv(csv_name, header=5)
+df = df.set_index(['Beach', 'Block', 'Profile', 'Year/Date'])
+
+# Get unique survey parameters
+beaches = df.index.get_level_values('Beach').unique()
+blocks = df.index.get_level_values('Block').unique()
+profiles = df.index.get_level_values('Profile').unique()
+dates = df.index.get_level_values('Year/Date').unique()
+
+# Sort dataframe
+df = df.sort_index()
+
+contours = []
+for beach in beaches:
+    for block in blocks:
+        for date in dates:
+            x = []
+            y = []
+            for profile in profiles:
+                try:
+                    # Get data for current profile
+                    survey = df.loc[(beach, block, profile, date), :]
+                    survey = survey.set_index('Chainage')
+
+                    # Find last point in profile above ref contour elevation
+                    last_idx = survey.where(
+                        survey['Elevation'] > contour_z).last_valid_index()
+                    elevation = survey.loc[last_idx:, 'Elevation']
+                    eastings = survey.loc[last_idx:, 'Easting']
+                    northings = survey.loc[last_idx:, 'Northing']
+                    x.append(np.interp(contour_z, elevation, eastings))
+                    y.append(np.interp(contour_z, elevation, northings))
+
+                except KeyError:
+                    pass
+
+            if x and y:
+                # Join points from same date into a line
+                line = LineString([(x1, y1) for x1, y1 in zip(x, y)])
+                contours.append({
+                    'date': date,
+                    'beach': beach,
+                    'block': block,
+                    'geometry': line
+                })
+
+# Make geodataframe
+gdf = gpd.GeoDataFrame(contours)
+
+# Set coordinate reference system
+gdf.crs = fiona.crs.from_epsg(epsg)
+
+beach_name = csv_name.replace('photogrammetry_Xsections_', '')
+shp_name = beach_name.replace('.csv', f' {contour_z}m contour.shp')
+gdf.to_file(shp_name, driver='ESRI Shapefile')