Volume Under "plane" Defined By Data Points - Python

I have a large mesh grid of data points which I have produced from simulations, and associated with each point in the xy plane is a z value (the outcome of the simulation). I have

Solution 1:

You could try integral() method of scipy.interpolate.RectBivariateSpline().

Solution 2:

If you want to strictly stick to the trapezoidal 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()

Whether spline or linear (trapezodial) interpolation is a better fit will depend heavily on your problem.

Solution 3:

Joe Kington's answer is almost good (and highly performant) but not quite correct. Here is the correct code, using the @ operator to keep the operations at the correct level with the full numpy performance.

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]
    vol = np.cross(a, b) @ tri[:,:,2]
    return vol.sum() / 6.0

main()