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.

score 5 · Answer 1 · answered Jan 09 '12 at 18:03

5

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

answered Jan 09 '12 at 18:03

jfs

399,953
195
994
1,670

tried it out, and it gets ~ the same answer as my very rough approximation, so I'm happy. Thanks – samb8s Jan 09 '12 at 19:01

score 3 · Accepted Answer · answered Jan 09 '12 at 23:10

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.

score 2 · Answer 3 · answered Apr 06 '21 at 13:01

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

Volume under "plane" defined by data points - python

3 Answers3

Linked