I have 2D data sampled at irregular points. I want to estimate the integral of this data over a finite x-y region. Is there a standard way to do this in python/scipy, or should I roll my own integration? The algorithm I imagine is to use Delaunay triangulation to construct a surface f(x, y), and then sum the volume under this approximate surface. The form my data currently takes: x = [1.1, 1.3, 1.31, 1.33, 2, 2.05] #x-coordinates at which data is known y = [0.15, 0.7, 0.01, 0.01, 0.9, 0] #y-coordinates at which data is known f = [1.1, 1.3, 1.15, 1.2, 1.18, 1.3] # Data at these x-y coordinates (This isn't real data, just meant to concretely illustrate the type of data I have) I have to process roughly a million of these datasets to produce a single image, so a vectorized solution would be nice. A version of my question is described here: http://math.stackexchange.com/questions/187730/efficiently-estimate-a-2d-int...