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

小编典典

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

python

我有一个很大的网格数据网格，是通过模拟生成的，并且与xy平面中的每个点相关联的是z值（模拟的结果）。

我将x，y，z值转储到纯文本文件中，我想做的是测量xy平面（即z =
0）和数据点定义的“平面”之间的体积。数据点目前还没有均匀地间隔开，尽管它们应该在模拟运行完成后才出现。

我一直在浏览scipy文档，不确定scipy.integrate是否提供我需要的功能-似乎只有2d的能力，而我没有3d的能力。

首先，除非有必要，否则我可以不用插值，而纯粹基于“梯形规则”或类似近似的积分是一个很好的起点。

任何帮助表示赞赏。

谢谢

编辑：下面介绍的两种解决方案都可以正常工作。以我为例，事实证明使用样条曲线可能会在平面上的尖锐最大值附近造成“波纹”，因此Delaunay方法效果更好，但我建议人们同时检查一下。

阅读 108

2021-01-20

共1个答案

小编典典

如果要严格遵守梯形法则，可以执行以下操作：

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

样条曲线或线性（梯形）插值是否更合适，在很大程度上取决于您的问题。

2021-01-20