在“平面"下的体积为0.由数据点定义-python [英] Volume under "plane" defined by data points - python

问题描述

我有一个很大的网格网格，这些网格是我从模拟中生成的，并且与xy平面中的每个点相关联的是z值(模拟的结果).

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

我将x，y，z值转储到纯文本文件中，我想做的是测量xy平面(即z = 0)和由定义的平面"之间的体积数据点.尽管应该在模拟运行完成后立即使用数据点，但这些数据点目前尚未均匀分布.

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.

我一直在浏览scipy文档，不确定scipy.integrate是否提供我需要的功能-似乎只能在2d中执行此操作，而不是我需要的3d.

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.

感谢您的帮助.

谢谢

下面介绍的两种解决方案都可以正常工作.就我而言，事实证明，使用样条曲线可能会导致平面上的尖锐最大值周围出现涟漪"，因此Delaunay方法效果更好，但我建议人们同时检查一下.

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.

推荐答案

如果您要严格遵守梯形规则，则可以执行以下操作:

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.

这篇关于在“平面"下的体积为0.由数据点定义-python的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！