Volumen en "plano" definido por puntos de datos - python

Question

Tengo una malla grande de puntos de datos que he producido a partir de simulaciones, y asociada a cada punto en el plano xy es un valor z (el resultado de la simulación).

Tengo los valores x, y, z vertidos en un archivo de texto plano, y lo que me gustaría hacer es medir el volumen entre el plano xy (es decir, z = 0) y el "plano" definido por los puntos de datos. Los puntos de datos no están actualmente espaciados uniformemente, aunque DEBERÍAN ser una vez que las simulaciones hayan terminado de ejecutarse.

He estado buscando a través de la documentación scipy, y no estoy seguro de si scipy.integrate proporciona la funcionalidad que necesito, parece que solo hay la posibilidad de hacer esto en 2d, no en 3D, como lo necesito.

Para empezar, a menos que sea necesario, puedo hacerlo sin interpolación, la integración basada únicamente en la "regla del trapecio" o una aproximación similar es una buena base para comenzar.

Cualquier ayuda es apreciada.

gracias

EDIT: tanto las soluciones descritas a continuación funcionan bien. En mi caso, resulta que usar una spline puede causar "ondulaciones" alrededor de máximos agudos en el avión, por lo que el método de Delaunay funciona mejor, pero aconsejo a las personas que comprueben ambas.

Answer 1

Si desea atenerse estrictamente a la regla trapezoidal se puede hacer algo similar a esto:

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

Si spline o lineal (dirección trapezoidal) La interpolación es un mejor ajuste dependerá en gran medida de su problema.

Answer 2

Puede probar el método integral() de scipy.interpolate.RectBivariateSpline().