Extrapolate towards land when points are missing

nielsen_volumes.py

@@ -16,12 +16,12 @@ def insertZerochainage(chainage, 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)
@@ -35,7 +35,7 @@ def insertZerochainage(chainage, elevation):
             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:
@@ -44,7 +44,7 @@ def insertZerochainage(chainage, elevation):
                     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)
@@ -53,27 +53,27 @@ def insertZerochainage(chainage, elevation):
                 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:
@@ -84,29 +84,29 @@ def findIntersection(x11,y11,x12,y12,x21,y21,x22,y22):
         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
@@ -119,7 +119,7 @@ def findZerochainage(chainage0, elevation0, max_beach_slope=5, min_beach_slope=5
             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]
@@ -139,50 +139,50 @@ def findZerochainage(chainage0, elevation0, max_beach_slope=5, min_beach_slope=5
             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) 
-        
+            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_intr=scipy.interpolate.interp1d(chainage, elevation, fill_value='extrapolate')
     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