From bf9fd6545d40bd2ccb8ca9d88928aec77236b2c7 Mon Sep 17 00:00:00 2001
From: Dan Howe <d.howe@wrl.unsw.edu.au>
Date: Tue, 27 Aug 2019 12:58:02 +1000
Subject: [PATCH] Improve robustness of shoreline detection

---
 extract_contours.py | 52 ++++++++++++++++++++++++++++-----------------
 1 file changed, 32 insertions(+), 20 deletions(-)

diff --git a/extract_contours.py b/extract_contours.py
index 5a271f4..8093601 100644
--- a/extract_contours.py
+++ b/extract_contours.py
@@ -5,7 +5,7 @@ import numpy as np
 import pandas as pd
 import geopandas as gpd
 import matplotlib.pyplot as plt
-from shapely.geometry import LineString
+from shapely.geometry import Point, LineString
 
 # Set contour elevation
 contour_z = 0.7
@@ -32,32 +32,44 @@ contours = []
 for beach in beaches:
     for block in blocks:
         for date in dates:
-            x = []
-            y = []
+            points = []
             for profile in profiles:
                 try:
-                    # Get data for current profile
+                    # Get data for current block/profile (if it exists)
                     survey = df.loc[(beach, block, profile, date), :]
-                    survey = survey.set_index('Chainage')
+                except KeyError:
+                    continue
 
-                    # Reverse survey chainage so elevation is increasing
-                    survey = survey[::-1]
+                # Reverse profile direction so elevation is increasing
+                survey = survey[::-1]
+                survey = survey.set_index('Chainage')
 
-                    # Find largest chainage in profile above contour elevation
-                    idx = survey.where(
-                        survey['Elevation'] > contour_z).first_valid_index()
-                    elevation = survey.loc[:idx, 'Elevation']
-                    eastings = survey.loc[:idx, 'Easting']
-                    northings = survey.loc[:idx, 'Northing']
-                    x = np.interp(contour_z, elevation, eastings)
-                    y = np.interp(contour_z, elevation, northings)
+                # Ignore profile information past high point
+                high_idx = survey['Elevation'].idxmax()
+                survey = survey[:high_idx]
 
-                except KeyError:
-                    pass
+                # Extract coordinates for current profile
+                elevation = survey['Elevation'].values
+                eastings = survey['Easting'].values
+                northings = survey['Northing'].values
+
+                # Check if entire profile is above reference elevation
+                if np.all(elevation > contour_z):
+                    continue
+
+                # Find first points on either side of reference elevation
+                crossing_idx = np.where(elevation < contour_z)[0][-1]
+                idx = slice(crossing_idx, crossing_idx + 2)
+
+                # Calculate shoreline intersection coordinates
+                x = np.interp(contour_z, elevation[idx], eastings[idx])
+                y = np.interp(contour_z, elevation[idx], northings[idx])
+                if x and y:
+                    points.append(Point(x, y))
 
-            if x and y:
-                # Join points from same date into a line
-                line = LineString([(x1, y1) for x1, y1 in zip(x, y)])
+            # Join points from same date into a line
+            if len(points) > 1:
+                line = LineString(points)
                 contours.append({
                     'date': date,
                     'beach': beach,