由数据点定义的“平面”下的体积 - python

Question

我有一个很大的网格数据，都是通过模拟生成的。每个在xy平面上的点都有一个对应的z值（就是模拟的结果）。

我把这些x、y、z值都放在了一个普通的文本文件里。现在我想要测量xy平面（也就是z=0）和这些数据点所定义的“平面”之间的体积。现在这些数据点的分布不均匀，不过等模拟完成后，它们应该会均匀分布。

我查了一下scipy的文档，但不太确定scipy.integrate是否能满足我的需求——看起来它只能处理二维的，而我需要的是三维的。

一开始，如果不必要的话，我可以不使用插值，单纯用“梯形法则”或类似的近似方法来进行积分，这样也是一个不错的起点。

任何帮助都非常感谢。

谢谢

补充：下面提到的两种解决方案都很好。在我的情况下，使用样条插值会在平面上尖锐的最大值附近产生“波纹”，所以Delaunay方法效果更好，但我建议大家都试试这两种方法。

Answer 1

Joe Kington的回答几乎是正确的（而且性能很高），但还是有点问题。这里是正确的代码，使用了@运算符，这样可以保持操作在正确的层级，同时还能充分发挥numpy的性能。

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

Answer 2

你可以试试 integral() 这个方法，它是属于 scipy.interpolate.RectBivariateSpline() 的。

Answer 3

如果你想严格按照梯形法则来做，可以尝试类似下面的做法：

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

使用样条插值还是线性插值（梯形法）哪个更合适，主要还是要看你的具体问题。