diff --git a/nielsen_volumes.py b/nielsen_volumes.py index 81a9389..858ac58 100644 --- a/nielsen_volumes.py +++ b/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=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