central-coast-aerial-survey/nielsen_volumes.py

#
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