Your data points lie on a regular grid, right? Than you can apply the trapezoidal rule directly to the data set:
First, I have to generate a fictive data set: pts
are the coordinates in the plane and z
denotes the elevation of the surface (yes, I presume that your surface is actually a graph).
xmin = 0.; xmax = 1.; xn = 1000;
ymin = 0.; ymax = 1.; yn = 1000;
pts = Tuples[{Subdivide[xmin, xmax, xn], Subdivide[ymin, ymax, yn]}];
f = {x, y} \[Function] x^2 x Sin[5 x + 3 y] + 1/2 Sin[7 x + 13 Pi y];
z = f @@ Transpose[pts];
This computes the integral with the two-dimensional trapezoidal rule (utilizing the fact that the data lies on a tensor product grid). If your surface is of class $C^2$, then the error should be proportional to the square of the diagional of the grid quadrilaterals.
xω = ConstantArray[1., xn + 1]; xω[[1]] = 0.5;
xω[[-1]] = 0.5;
yω = ConstantArray[1., yn + 1]; yω[[1]] = 0.5;
yω[[-1]] = 0.5;
int1 = (xmax - xmin)/xn (ymax - ymin)/yn (yω.Partition[z, xn + 1].xω);
-0.0742535
For checking the accuracy, here is the same integral computed with NIntegrate
:
int2 = NIntegrate[ f[x, y], {x, xmin, xmax}, {y, xmin, xmax}, AccuracyGoal -> 10, Method -> "GaussKronrodRule" ]
-0.0742538
There is definitely some error, but for inaccurate input data, this is really negligble:
Abs[int1 - int2]/Abs[int2]
3.96486*10^-6
General hint
Make sure that your data z = data[[All,3]]
is a packed array of machine precision numbers, for example with z = Deverloper`ToPackedArray[data[[All,3]]]
.