From 8ca180f93f073a8920d28e990f47db73495258c2 Mon Sep 17 00:00:00 2001
From: Dan Howe <d.howe@wrl.unsw.edu.au>
Date: Thu, 19 Jul 2018 08:56:48 +1000
Subject: [PATCH] Add nielsen_volumes.py

---
 las_outputs.py     |   2 +-
 nielsen_volumes.py | 188 +++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 189 insertions(+), 1 deletion(-)
 create mode 100644 nielsen_volumes.py

diff --git a/las_outputs.py b/las_outputs.py
index 04dfce1..ec6adb2 100644
--- a/las_outputs.py
+++ b/las_outputs.py
@@ -22,7 +22,7 @@ from cycler import cycler
 import matplotlib.pyplot as plt
 from matplotlib.ticker import MultipleLocator
 
-import neilson_volumes
+import nielsen_volumes
 from survey_tools import call_lastools, extract_pts, update_survey_output
 
 
diff --git a/nielsen_volumes.py b/nielsen_volumes.py
new file mode 100644
index 0000000..81a9389
--- /dev/null
+++ b/nielsen_volumes.py
@@ -0,0 +1,188 @@
+#
+import scipy.interpolate
+import numpy as np
+import copy
+import math
+
+
+repose_angle=34 #angle of repose of sand
+repose_slope=1/math.tan(math.radians(repose_angle))
+
+def insertZerochainage(chainage, elevation):
+    #this inserts a 0 every time the profile crosses the axis
+    #need to update the chainage elevation for every point that drops below the required chainage
+    ii=0
+    while ii<len(elevation):
+        if ii==0 and elevation[ii]<0:
+            while elevation[ii]<0:
+                ii=ii+1
+                
+            try:
+                beachslope = 1/((elevation[ii-1] - elevation[ii])/(chainage[ii] - chainage[ii-1]))
+            except:
+                beachslope= 999
+                
+            chainageX= (elevation[ii-1]-0)*beachslope+chainage[ii-1]
+            elevation.insert(ii,0)
+            chainage.insert(ii,chainageX)
+            ii=ii+1
+        if elevation[ii]<0:
+            try:
+                beachslope = 1/((elevation[ii-1] - elevation[ii])/(chainage[ii] - chainage[ii-1]))
+            except ZeroDivisionError:
+                beachslope= 999
+            chainageX= (elevation[ii-1]-0)*beachslope+chainage[ii-1]
+            elevation.insert(ii,0)
+            chainage.insert(ii,chainageX)
+            ii=ii+1
+            
+            try:
+                #finding where it crosses back over zero.  If you reach the end of the list, return the list
+                while elevation[ii]<0:
+                    ii=ii+1
+                try:
+                    beachslope = 1/((elevation[ii-1] - elevation[ii])/(chainage[ii] - chainage[ii-1]))
+                except ZeroDivisionError:
+                    beachslope= 999
+                    
+                chainageX= (elevation[ii-1]-0)*beachslope+chainage[ii-1]
+                elevation.insert(ii,0)
+                chainage.insert(ii,chainageX)
+                ii=ii+1
+            except IndexError:
+                return chainage, elevation
+        else:
+            ii=ii+1
+            
+    return chainage, elevation
+    
+def custom_trap_rule(elevation,chainage):
+    #so that it doesn't count any volume below zero as a negative
+    
+    #add a zero point every time the profile crosses the axis (note that the seaward zero has already been inserted)
+    x,y=insertZerochainage(copy.copy(chainage), copy.copy(elevation))
+    # add the final 0
+    
+    x,y = findZerochainage(x,y)
+    
+    i=1
+    area=0
+    while i<len(y):
+        if y[i]>=0 and y[i-1]>=0:
+            area=area+0.5*(x[i]-x[i-1])*(y[i]+y[i-1])
+        i=i+1
+        
+    return area
+    
+def findIntersection(x11,y11,x12,y12,x21,y21,x22,y22):
+    #find the intersection of a two lines defined by four points
+    try:
+        m1=(y12-y11)/(x12-x11)
+        b1=y11-m1*x11
+    except (ZeroDivisionError,FloatingPointError):
+        x_inter=x11
+        m2=(y22-y21)/(x22-x21)
+        b2=y21-m2*x21
+        y_inter=m2*x_inter+b2
+        
+        return x_inter, y_inter
+    
+    try:
+        m2=(y22-y21)/(x22-x21)
+        b2=y21-m2*x21
+    except (ZeroDivisionError, FloatingPointError):
+        x_inter=x21
+        y_inter=m1*x_inter +b1
+        
+        return x_inter, y_inter
+    
+    x_inter=(b2-b1)/(m1-m2)
+    y_inter=m1*x_inter+b1
+    
+    return x_inter, y_inter
+    
+def findZerochainage(chainage0, elevation0, max_beach_slope=5, min_beach_slope=50):
+	#often the photogrammetry has a repeated point at the end.  Check for this and remove extra if required
+    while chainage0[-1]==chainage0[-2] and elevation0[-1]==elevation0[-2]:
+        chainage0 = chainage0[:-1]
+        elevation0 = elevation0[:-1]
+    
+    #check if the profile ever goes below zero
+    ii=len(elevation0)-1
+    below_0=False
+    while ii>len(elevation0)/2:
+        if elevation0[ii]<0:
+            below_0=True
+            while elevation0[ii]<0:
+                ii=ii-1
+            j=ii+1 #so you don't have to search for this again
+            ii=0
+        else:
+            ii=ii-1
+    
+    # if Below_0 = True, we can use the actual zero value
+    if below_0==True:
+        chainage0=chainage0[0:j+1]
+        elevation0=elevation0[0:j+1]
+        elevation_inter=scipy.interpolate.interp1d(elevation0[-2:],chainage0[-2:])
+        chainage_0=float(elevation_inter(0))
+        chainage0[j]=chainage_0
+        elevation0[j]=0
+    #otherwise need to extrapolate
+    else:
+        #find the slope between the last two points
+        try:
+            beachslope = 1/((elevation0[-2] - elevation0[-1])/(chainage0[-1] - chainage0[-2]))
+        except ZeroDivisionError:
+            beachslope= 999
+        if beachslope>min_beach_slope:
+            beachslope = min_beach_slope
+        elif beachslope<max_beach_slope:
+            beachslope = max_beach_slope
+            
+        chainageX=(elevation0[-1])*beachslope+chainage0[-1]
+        if chainageX!=chainage0[-1]:
+            chainage0.append(chainageX)
+            elevation0.append(0) 
+        
+    return chainage0,elevation0
+    
+
+def insert_LL(chainage, elevation, LL):
+    # insert an additional element to make LL an actual point
+    elev_intr=scipy.interpolate.interp1d(chainage, elevation)
+    elev_LL=float(elev_intr(LL))
+    data=[[chainage[i], elevation[i]] for i in range(0, len(chainage))]
+    
+    data.append([LL,elev_LL])
+    data.sort()
+        
+    chainage=[data[i][0] for i in range(0, len(data))]
+    elevation=[data[i][1] for i in range(0, len(data))]
+    
+    return chainage, elevation, elev_LL
+    
+def volume_available (chainage, elevation, LL):
+    chainage, elevation, LL_z = insert_LL(chainage, elevation, LL)
+    
+    index_LL=0
+    while chainage[index_LL]!=LL:
+        index_LL=index_LL+1
+        
+    bottom_ZSA_x = LL + LL_z/repose_slope
+    bottom_ZSA_y=0
+    
+    mid_ZWI_x= LL+(LL_z - 2)/(2*math.tan(math.radians(repose_angle)))
+    mid_ZWI_y = 2
+    
+    bottom_ZWI_x=mid_ZWI_x+ mid_ZWI_y*10
+    bottom_ZWI_y = 0
+    
+    x,y=findIntersection(LL,LL_z,bottom_ZSA_x,bottom_ZSA_y, mid_ZWI_x, mid_ZWI_y, bottom_ZWI_x, bottom_ZWI_y )
+    
+    V1=custom_trap_rule(elevation[index_LL:], chainage[index_LL:])
+    V2=custom_trap_rule([ LL_z,mid_ZWI_y, bottom_ZWI_y],[LL,mid_ZWI_x,bottom_ZWI_x])
+    
+    volume=V1-V2
+
+    return volume