Volume sotto "piano" definito da punti di dati - python

Question

Ho una griglia di punti di dati di grandi dimensioni che ho prodotto da simulazioni e associato a ciascun punto nel piano xy è un valore z (il risultato della simulazione).

Ho i valori x, y, z riversati in un file di testo normale e quello che mi piacerebbe fare è misurare il volume tra il piano xy (cioè z = 0) e il "piano" definito da i punti dati. I punti dati non sono attualmente distribuiti uniformemente, sebbene DOVREBBE essere una volta terminate le simulazioni.

Ho esaminato la documentazione scipy e sono incerto se scipy.integrate fornisca la funzionalità di cui ho bisogno: sembra che ci sia solo la possibilità di farlo in 2d, non in 3d come mi serve.

Per cominciare, se non necessario, posso fare a meno dell'interpolazione, l'integrazione basata esclusivamente sulla "regola del trapezio" o un'approssimazione simile è una buona base da cui partire.

Qualsiasi aiuto è apprezzato.

grazie

EDIT: entrambe le soluzioni descritte di seguito funzionano bene. Nel mio caso, risulta che l'uso di una spline può causare "increspature" intorno ai massimi affilati nel piano, quindi il metodo Delaunay funziona meglio, ma consiglierei alle persone di controllarli entrambi.

Answer 1

Se si vuole attaccare rigorosamente la regola trapezoidale che si può fare qualcosa di simile a questo:

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

Sia spline o lineare (guida trapezoidale) interpolazione è una misura migliore dipenderà in larga misura il problema.

Answer 2

Si potrebbe provare il metodo integral() di scipy.interpolate.RectBivariateSpline().