Geek Answers Handbook

Volume under a “plane” defined by data points - python

I have a large grid of data point grid that I created from simulations, and associated with each point in the xy plane is the value z (simulation result).

I have x, y, z values ​​dumped into a text file, and what I would like to do is measure the volume between the xy plane (ie z = 0) and the “plane” defined by the data. Data points are currently not evenly spaced, although they MUST be after the simulations have ended.

I look through the scant documentation and I'm not sure that scipy.integrate provides the functions I need - it seems that there is only the ability to do this in 2d, and not in 3d, as I need.

To begin with, if necessary, I can do without interpolation, integration based solely on the "trapezoid rule" or a similar approximation is a good basis for a start.

Any help is appreciated.

thank

EDIT: Both of the solutions described below work well. In my case, it turns out that the spline can cause ripples around sharp maxima in the plane, so the Delaunay method works better, but I would advise people to check both.

+5
source share
2 answers

If you want to strictly adhere to the trapezoid rule, you can do something similar to this:

import numpy as np
import scipy.spatial

def main():
    xyz = np.random.random((100, 3))
    area_underneath = trapezoidal_area(xyz)
    print area_underneath

def trapezoidal_area(xyz):
    """Calculate volume under a surface defined by irregularly spaced points
    using delaunay triangulation. "x,y,z" is a <numpoints x 3> shaped ndarray."""
    d = scipy.spatial.Delaunay(xyz[:,:2])
    tri = xyz[d.vertices]

    a = tri[:,0,:2] - tri[:,1,:2]
    b = tri[:,0,:2] - tri[:,2,:2]
    proj_area = np.cross(a, b).sum(axis=-1)
    zavg = tri[:,:,2].sum(axis=1)
    vol = zavg * np.abs(proj_area) / 6.0
    return vol.sum()

main()

If spline or linear (trapezoidal) interpolation is better suited, it will greatly depend on your problem.

+3
source

More articles:


All Articles