Playground notebook for volume calculations


In [29]:
# imports
import math as m
import numpy as np
import pandas as pd
import pickle
import os
from scipy import integrate

In [2]:
# creating dataframe from csv file
print('Loding survey data')
dfCoordinates = pd.read_csv('gpsDataUTM.csv') # reads data into dataframe


Loding survey data

In [3]:
dfCoordinates.head()


Out[3]:
x y z zoneLetter zoneNumber x_rel y_rel z_rel
0 460564.162280 7.242060e+06 883.16 J 22.0 213.860510 127.527218 14.12
1 460564.002393 7.242060e+06 883.16 J 22.0 213.700623 127.211200 14.12
2 460564.321028 7.242032e+06 883.18 J 22.0 214.019258 98.792774 14.14
3 460564.563586 7.242003e+06 880.68 J 22.0 214.261816 70.346291 11.64
4 460567.438707 7.241978e+06 878.91 J 22.0 217.136937 44.608228 9.87

In [4]:
# Adding new info to dataframe: distance from origin (common reference to all points)
def distanceFromOrigin(x,y,z):
    # compute distance from origin <0,0,0>
    point = np.array([x,y,z]) # place columns into an array <x,y,z>
    origin = np.zeros(3) # create origin <0,0,0>
    return np.linalg.norm(point-origin) # compute norm

dfCoordinates['oDist'] = np.vectorize(distanceFromOrigin)(dfCoordinates['x_rel'],
                                                          dfCoordinates['y_rel'],
                                                          dfCoordinates['z_rel'])

In [5]:
dfCoordinates.head()


Out[5]:
x y z zoneLetter zoneNumber x_rel y_rel z_rel oDist
0 460564.162280 7.242060e+06 883.16 J 22.0 213.860510 127.527218 14.12 249.397040
1 460564.002393 7.242060e+06 883.16 J 22.0 213.700623 127.211200 14.12 249.098414
2 460564.321028 7.242032e+06 883.18 J 22.0 214.019258 98.792774 14.14 236.144436
3 460564.563586 7.242003e+06 880.68 J 22.0 214.261816 70.346291 11.64 225.814562
4 460567.438707 7.241978e+06 878.91 J 22.0 217.136937 44.608228 9.87 221.891326

In [6]:
# sorts accordint to distance from the origin.
dfCoordinates = dfCoordinates.sort_values(by=['oDist'],
                                          ascending=True)

In [11]:
dfCoordinates.head()


Out[11]:
x y z zoneLetter zoneNumber x_rel y_rel z_rel oDist
19 460351.568925 7.241944e+06 877.68 J 22.0 1.267155 10.898245 8.64 13.965207
20 460350.322685 7.241947e+06 878.51 J 22.0 0.020915 14.597884 9.47 17.400562
18 460368.718418 7.241934e+06 876.81 J 22.0 18.416649 0.664109 7.77 19.999672
17 460370.150521 7.241933e+06 876.80 J 22.0 19.848751 0.000000 7.76 21.311746
72 460377.487691 7.241936e+06 880.45 J 22.0 27.185922 2.907675 11.41 29.626289

In [30]:
def computeLineEquation(A,B):
    # computes line equation given 2 points: A and B
    # Both A and B need to be numpy arrays
    
    if (B[0]-A[0]) == 0: # Check if line is parallel to y axis (zero divide)
        #print("Zero division. Line is parallel to y axis. x={x}".format(x=A[0]))
        return np.array([None,None])
    else: # detect slope and intersection
        m = (B[1]-A[1])/(B[0]-A[0]) # slope
        xo = A[0]
        yo = A[1]
        c = -m*xo+yo    
        #print("The line equation is y = {m}x + {c}".format(m=m,c=c))
        return np.array([m,c]) # returns slope and constant in numpy array format

In [150]:
def computeVolume(A,B,C):
    # computes volume given 3 points using double integral, line and plane equations
    # A, B and C are numpy arrays
    
    # Setting up plane equation
    AB = B-A # vector within plane
    BC = C-B # vector within plane
    NV = np.cross(AB,BC) # vector that is normal to plane
    
    # Plane values; plane equation is a*(x-xo)+b*(y-yo)+c*(z-zo)=0
    a = NV[0]
    b = NV[1]
    c = NV[2]
    xo = A[0]
    yo = A[1]
    zo = A[2]
    
    # Compute line equations:
    slopeAB, constantAB = computeLineEquation(A,B)
    slopeBC, constantBC = computeLineEquation(B,C)
    slopeAC, constantAC = computeLineEquation(A,C)
    
    # Compute double integral
    if ((slopeAC is not None) and (slopeAB is not None)): # if slopes are defined
        volume1, error1 = integrate.dblquad(lambda y, x: ((-a*(x-xo)-b*(y-yo))/c)+zo, 
                                            A[0],
                                            B[0], 
                                            lambda x: slopeAB*x+constantAB,
                                            lambda x: slopeAC*x+constantAC)
    else:
        volume1 = 0.0
    if ((slopeAC is not None) and (slopeBC is not None)): # if slopes are defined
        volume2, error2 = integrate.dblquad(lambda y, x: ((-a*(x-xo)-b*(y-yo))/c)+zo, 
                                            B[0],
                                            C[0], 
                                            lambda x: slopeAC*x+constantAC,
                                            lambda x: slopeBC*x+constantBC)
    else:
        volume2 = 0.0
    return abs(volume1)+abs(volume2) # return total volume

In [156]:
# Required number of loops to compute total volume = number of points - 1

totalVolume = 0.0 # initialize volume summation holder
for i in range(dfCoordinates['x'].count()-2): # iterate thru all groups of 3 points in ascending order from origin.
    s1 = dfCoordinates.iloc[i] # get series (row entry)
    s2 = dfCoordinates.iloc[i+1]
    s3 = dfCoordinates.iloc[i+2]
    
    dfAnalysis = pd.DataFrame([]) # build dataframe composed of these 3 entry set
    dfAnalysis = dfAnalysis.append(s1,ignore_index=True)
    dfAnalysis = dfAnalysis.append(s2,ignore_index=True)
    dfAnalysis = dfAnalysis.append(s3,ignore_index=True)
    dfAnalysis = dfAnalysis.sort_values(by=['x_rel'], ascending=True)# order by lowest x_rel >> A,B,C
    
    A = np.array([dfAnalysis.iloc[0]['x_rel'],
                  dfAnalysis.iloc[0]['y_rel'],
                  dfAnalysis.iloc[0]['z_rel']])
    B = np.array([dfAnalysis.iloc[1]['x_rel'],
                  dfAnalysis.iloc[1]['y_rel'],
                  dfAnalysis.iloc[1]['z_rel']])
    C = np.array([dfAnalysis.iloc[2]['x_rel'],
                  dfAnalysis.iloc[2]['y_rel'],
                  dfAnalysis.iloc[2]['z_rel']])
    
    # Compute line equations, plane equations and volume
    volume = computeVolume(A,B,C)
    totalVolume += volume # add to total volume
print('Total volume is: {v} cubic meters.'.format(v=totalVolume))


Total volume is: 768401.7102283855 cubic meters.

In [160]:
18.5*194*220


Out[160]:
789580.0