Volume under "plane" defined by data points - python

Question

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 the x,y,z values dumped into a plain text file, and what I'd like to do is measure the volume between the xy plane (ie. z=0) and the "plane" defined by the data points. The data points are not currently uniformly spaced, although they SHOULD be once the simulations have finished running.

I've been looking through the scipy documentation, and I'm uncertain whether scipy.integrate provides the functionality I need - it seems that there is only the ability to do this in 2d, not 3d as I need.

To begin with, unless necessary, I can do without interpolation, integration based purely on the "trapezium rule" or similar approximation is a good basis to start from.

Any help is appreciated.

thanks

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

Answer 1

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

Answer 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.

Answer 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()