Playground notebook for volume calculations


In [1]:
# imports
import math as m
import numpy as np
import pandas as pd
import pickle
import os

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]:
# Example: computing distance between points
A = dfCoordinates[['x_rel','y_rel','z_rel']].values[0]
B = dfCoordinates[['x_rel','y_rel','z_rel']].values[50]
O = np.array([0.0,0.0,0.0]) # origin

In [5]:
distanceAB = np.linalg.norm(B-A)
distanceOA = np.linalg.norm(A-O)
distanceOB = np.linalg.norm(B-O)
print('Distance AB: %s' % distanceAB)
print('Distance OA: %s' % distanceOA)
print('Distance OB: %s' % distanceOB)


Distance AB: 62.0909745878
Distance OA: 249.397040219
Distance OB: 191.342243601

In [6]:
# Adding new info to dataframe
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 [7]:
dfCoordinates.head()


Out[7]:
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 [8]:
# sorts accordint to distance from the origin.
dfCoordinates = dfCoordinates.sort_values(by=['oDist'],
                                          ascending=True)

In [9]:
pd.__version__


Out[9]:
'0.20.1'

In [10]:
dfCoordinates.head()


Out[10]:
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 [11]:
m.sqrt(m.pow(19.848751,2)+m.pow(7.76,2)) # check!


Out[11]:
21.311745969300613

Learning to compute line equations


In [15]:
# Useful computing. Compute line equation given 2 points A, B.
A = np.array([1,2])
B = np.array([4,5])
# line equation: (y-yo) = m*(x-xo)
m = (A[1]-B[1])/(A[0]/B[0]) # slope
xo = A[0]
yo = A[1]
c = -m*xo+yo
print("Line Solution is y = {m}x + {c}".format(m=m,c=c))


Line Solution is y = -12.0x + 14.0

In [9]:
# Defining a generic function to comput the above
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 [11]:
slope, constant = computeLineEquation(np.array([1,2]), np.array([1,15]))


Zero division. Line is parallel to y axis. x=1

In [5]:
# the same but different...
slope, constant = computeLineEquation(np.array([1,2]), np.array([4,5]))


The line equation is y = 1.0x + 1.0

In [6]:
slope, constant = computeLineEquation(np.array([4,5]), np.array([10,1]))


The line equation is y = -0.6666666666666666x + 7.666666666666666

In [7]:
slope, constant = computeLineEquation(np.array([1,2]), np.array([10,1]))


The line equation is y = -0.1111111111111111x + 2.111111111111111

Computing plane equations from 3 points


In [46]:
# Example
x = np.linspace(0.0,100.0,num=500) # generate linear space for x values
y = np.linspace(0.0,100.0,num=500)
xGrid, yGrid = np.meshgrid(x,y) # generate mesh grid for plotting sample data

point1 = np.array([[1.0,
                    2.0,
                    1.0]])
point2 = np.array([[4.0,
                    5.0,
                    0.0]])
point3 = np.array([[10.0,
                    1.0,
                    3.0]])

planeVectorOne = point2 - point1 # a vector within the plane
planeVectorTwo = point3 - point2 # a vector within the plane
perpendicularVector = np.cross(planeVectorOne, planeVectorTwo) # plane equation coeficients obtained <a,b,c>

# print(perpendicularVector)
# print(type(perpendicularVector))
# print(perpendicularVector.shape)
# print(point1.shape)
# print(perpendicularVector[0][2])
# saving constants to literals
a = perpendicularVector[0][0]
b = perpendicularVector[0][1]
c = perpendicularVector[0][2]
xo = point1[0][0]
yo = point1[0][1]
zo = point1[0][2]

# Plane equation a(x-xo)+b(y-yo)+c(z-zo)=0 with isolated z:
z = ((-a*(xGrid-xo)-b*(yGrid-yo))/c)+zo

Computing volume given 3 sets of points (x,y,z)


In [47]:
# Steps:
# 1. Compute plane equation from 3 points. 1 plane.
# 2. Compute line equations. 3 equations.
# 3. Sort them according to distance. A >> B >> C
# 4. Volume1 = Double integral xA-xB lineEqAC lineEqAB
# 5. Volume2 = Double integral xB-xC lineEqAC lineEqBC
# 6. VolumeTotal = Volume1 + Volume2

In [49]:
# Setting up sample space
x = np.linspace(0.0,11.0,num=50) # generate linear space for x values
y = np.linspace(0.0,11.0,num=50)
xGrid, yGrid = np.meshgrid(x,y) # generate mesh grid for plotting sample data

In [62]:
# Setting up points in 3D space
A = np.array([1.0,
              2.0,
              1.0])
B = np.array([4.0,
              5.0,
              0.0])
C = np.array([10.0,
              1.0,
              3.0])

In [63]:
# 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

In [71]:
# Plane values
a = NV[0]
b = NV[1]
c = NV[2]
xo = A[0]
yo = A[1]
zo = A[2]
print("The plane equation is {a}*(x-{xo})+{b}*(y-{yo})+{c}*(z-{zo})=0".format(a=a,b=b,c=c,xo=xo,yo=yo,zo=zo))
# With these values, you can fully define the plane equation


The plane equation is 5.0*(x-1.0)+-15.0*(y-2.0)+-30.0*(z-1.0)=0

In [68]:
# computing plane equations
slopeAB, constantAB = computeLineEquation(A,B)
slopeBC, constantBC = computeLineEquation(B,C)
slopeAC, constantAC = computeLineEquation(A,C)


The line equation is y = 1.0x + 1.0
The line equation is y = -0.6666666666666666x + 7.666666666666666
The line equation is y = -0.1111111111111111x + 2.111111111111111

In [78]:
# Compute volumes
from scipy import integrate
volume1, error1 = integrate.dblquad(lambda y, x: ((-a*(x-xo)-b*(y-yo))/c)+zo, 
                                    A[0],
                                    B[0], 
                                    lambda x: slopeAC*x+constantAC,
                                    lambda x: slopeAB*x+constantAB)
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)
volumeTotal = volume1+volume2
print("Total computed volume is {v}".format(v=volumeTotal))


Total computed volume is 19.999999999999993

Running example 2 on z=10


In [12]:
# Setting up points in 3D space
A = np.array([1.0,
              2.0,
              10.0])
B = np.array([3.0,
              5.0,
              10.0])
C = np.array([5.0,
              2.0,
              10.0])
# A = np.array([1.0,
#               2.0,
#               10.0])
# B = np.array([5.0,
#               1.0,
#               10.0])
# C = np.array([5.0,
#               5.0,
#               10.0])
# 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
a = NV[0]
b = NV[1]
c = NV[2]
xo = A[0]
yo = A[1]
zo = A[2]
print("The plane equation is {a}*(x-{xo})+{b}*(y-{yo})+{c}*(z-{zo})=0".format(a=a,b=b,c=c,xo=xo,yo=yo,zo=zo))
# With these values, you can fully define the plane equation
# computing plane equations
slopeAB, constantAB = computeLineEquation(A,B)
slopeBC, constantBC = computeLineEquation(B,C)
slopeAC, constantAC = computeLineEquation(A,C)
# Compute volumes
from scipy import integrate
volume1, error1 = integrate.dblquad(lambda y, x: ((-a*(x-xo)-b*(y-yo))/c)+zo, 
                                    A[0],
                                    B[0], 
                                    lambda x: slopeAC*x+constantAC,
                                    lambda x: slopeAB*x+constantAB)
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)
volumeTotal = volume1+volume2
print("Total computed volume is {v1}+{v2}={v}".format(v1=volume1,v2=volume2,v=volumeTotal))
# Volume Computed Correctly! :)
# Check if it works with line equations like x=3 -- parallel to y axis; Nope! Divide by 0 error. Analysing...
# -- when divide by zero occurs, corresponding volume is zero! :) find a way around this on code.


The plane equation is 0.0*(x-1.0)+0.0*(y-2.0)+-12.0*(z-10.0)=0
The line equation is y = 1.5x + 0.5
The line equation is y = -1.5x + 9.5
The line equation is y = 0.0x + 2.0
Total computed volume is 30.000000000000004+30.000000000000007=60.000000000000014

Running example with line parallel to y axis (slope=infinity)


In [46]:
# Setting up points in 3D space
A = np.array([1.0,
              3.0,
              10.0])
B = np.array([5.0,
              5.0,
              10.0])
C = np.array([11.0,
              3.0,
              10.0])
# 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
a = NV[0]
b = NV[1]
c = NV[2]
xo = A[0]
yo = A[1]
zo = A[2]
print("The plane equation is {a}*(x-{xo})+{b}*(y-{yo})+{c}*(z-{zo})=0".format(a=a,b=b,c=c,xo=xo,yo=yo,zo=zo))
# With these values, you can fully define the plane equation
# computing plane equations
slopeAB, constantAB = computeLineEquation(A,B)
slopeBC, constantBC = computeLineEquation(B,C)
slopeAC, constantAC = computeLineEquation(A,C)

# Compute volumes
from scipy import integrate
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
volumeTotal = abs(volume1)+abs(volume2)
print("Volume 1: {v}".format(v=volume1))
print("Volume 2: {v}".format(v=volume2))
print("Total computed volume is {v}".format(v=volumeTotal))
# Volume Computed Correctly! :)
# Check if it works with line equations like x=3 -- parallel to y axis; Nope! Divide by 0 error. Analysing...
# -- when divide by zero occurs, corresponding volume is zero! :) find a way around this on code.


The plane equation is 0.0*(x-1.0)+0.0*(y-3.0)+-20.0*(z-10.0)=0
The line equation is y = 0.5x + 2.5
The line equation is y = -0.3333333333333333x + 6.666666666666666
The line equation is y = 0.0x + 3.0
Volume 1: -40.0
Volume 2: 59.99999999999998
Total computed volume is 99.99999999999997

Testing volume calculation on known Pyramid


In [47]:
# Volume pyramid: 3x2x5/3 = 10

In [49]:
# Setting up points in 3D space
A = np.array([0.0,
              0.0,
              2.0])
B = np.array([0.0,
              3.0,
              2.0])
C = np.array([5.0,
              2.5,
              0.0])
# 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
a = NV[0]
b = NV[1]
c = NV[2]
xo = A[0]
yo = A[1]
zo = A[2]
print("The plane equation is {a}*(x-{xo})+{b}*(y-{yo})+{c}*(z-{zo})=0".format(a=a,b=b,c=c,xo=xo,yo=yo,zo=zo))
# With these values, you can fully define the plane equation
# computing plane equations
slopeAB, constantAB = computeLineEquation(A,B)
slopeBC, constantBC = computeLineEquation(B,C)
slopeAC, constantAC = computeLineEquation(A,C)

# Compute volumes
from scipy import integrate
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
volumeTotal = abs(volume1)+abs(volume2)
print("Volume 1: {v}".format(v=volume1))
print("Volume 2: {v}".format(v=volume2))
print("Total computed volume is {v}".format(v=volumeTotal))


The plane equation is -6.0*(x-0.0)+0.0*(y-0.0)+-15.0*(z-2.0)=0
Zero division. Line is parallel to y axis. x=0.0
The line equation is y = -0.1x + 3.0
The line equation is y = 0.5x + 0.0
Volume 1: 0.0
Volume 2: 10.000000000000002
Total computed volume is 10.000000000000002

In [50]:
# IT WORKED!

Computing volume of pyramid on top of cube


In [51]:
# Known volume: 25

In [52]:
# Setting up points in 3D space
A = np.array([0.0,
              3.0,
              4.0])
B = np.array([0.0,
              0.0,
              4.0])
C = np.array([5.0,
              2.5,
              2.0])
# 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
a = NV[0]
b = NV[1]
c = NV[2]
xo = A[0]
yo = A[1]
zo = A[2]
print("The plane equation is {a}*(x-{xo})+{b}*(y-{yo})+{c}*(z-{zo})=0".format(a=a,b=b,c=c,xo=xo,yo=yo,zo=zo))
# With these values, you can fully define the plane equation
# computing plane equations
slopeAB, constantAB = computeLineEquation(A,B)
slopeBC, constantBC = computeLineEquation(B,C)
slopeAC, constantAC = computeLineEquation(A,C)

# Compute volumes
from scipy import integrate
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
volumeTotal = abs(volume1)+abs(volume2)
print("Volume 1: {v}".format(v=volume1))
print("Volume 2: {v}".format(v=volume2))
print("Total computed volume is {v}".format(v=volumeTotal))


The plane equation is 6.0*(x-0.0)+0.0*(y-3.0)+15.0*(z-4.0)=0
Zero division. Line is parallel to y axis. x=0.0
The line equation is y = 0.5x + 0.0
The line equation is y = -0.1x + 3.0
Volume 1: 0.0
Volume 2: -25.0
Total computed volume is 25.0

So at the end of this chapter, we know how to calculate the volume of a solid object formed by a set of 3 points in space. Now it is a matter of ordering how this computations are to be made for a point cloud and sum all of them for total volume.