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

Question

我有一个大型的数据网格，由模拟产生的。在xy平面上的每个点都与一个z值相关联（模拟结果）。我将x、y、z值转储到纯文本文件中，我想做的是测量xy平面（即z=0）和由数据点定义的“平面”之间的体积。目前数据点不是均匀分布的，尽管一旦模拟运行完成，它们应该均匀分布。我查看了scipy文档，并不确定它是否提供我所需的功能——看来只能在2d中进行操作，而我需要的是3d。除非必要，否则首先可以不使用插值，基于“梯形法则”或类似近似的积分是一个很好的起点。感谢任何帮助。编辑：下面描述的两种方法都有效。在我的情况下，使用样条会导致平面上的尖峰周围出现波动，因此Delaunay方法更好，但我建议人们都试一下。

- samb8s

Answer 1

你可以尝试使用 integral() 方法，该方法属于 scipy.interpolate.RectBivariateSpline()。

- jfs

---

我尝试了一下，它得到的答案和我非常粗略的近似值差不多，所以我很满意。谢谢。 - samb8s

Answer 2

如果你想严格遵循梯形法则，可以尝试类似于以下方法：

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

无论是样条插值还是线性（梯形）插值，哪种更适合取决于您的问题。

- Joe Kington

Answer 3

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

- simonh_2