Объем под «плоскостью» определяется точками данных - питон - PullRequest
5 голосов
/ 09 января 2012

У меня есть большая сетка с точками данных, которые я сгенерировал из симуляций, и с каждой точкой в ​​плоскости xy связано значение z (результат симуляции).

У меня есть значения x, y, z, выгруженные в простой текстовый файл, и я хотел бы измерить громкость между плоскостью xy (т. Е. Z = 0) и плоскостью, определяемой Точки данных. Точки данных в настоящее время не распределены равномерно, хотя они ДОЛЖНЫ быть после завершения моделирования.

Я просматривал документацию по scipy, и я не уверен, предоставляет ли scipy.integrate необходимую мне функциональность - кажется, что есть возможность делать это только в 2D, а не в 3D, как мне нужно.

Начнем с того, что без необходимости я могу обойтись без интерполяции. Хорошей основой для начала является интеграция, основанная исключительно на «правиле трапеции» или подобном приближении.

Любая помощь приветствуется.

спасибо

РЕДАКТИРОВАТЬ: оба решения, описанные ниже, работают хорошо. В моем случае выясняется, что использование сплайна может вызвать «рябь» вокруг резких максимумов в плоскости, поэтому метод Делоне работает лучше, но я бы посоветовал людям проверить оба.

Ответы [ 2 ]

4 голосов
/ 09 января 2012

Вы можете попробовать integral() метод scipy.interpolate.RectBivariateSpline().

3 голосов
/ 10 января 2012

Если вы хотите строго придерживаться правила трапеции, вы можете сделать что-то похожее на это:

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

Подойдет ли сплайн или линейная (трапециевидная) интерполяция, будет сильно зависеть от вашей проблемы.

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...