I am trying to find the volume of a region (defined by X,Y and Z coordinates) enclosed below a (Z='constant') plane. The data has peaks (positive Z) and a valley (negative Z), and the mean surface is assigned z=0. I have tried with the following piece of code, but I doubt it gives me the total volume bound by the surface against z=0, including the peak volumes as well.
[X,Y,Z] = xyzread("data.xyz");
[fitobject, gof, output] = fit([X,Y],Z, 'biharmonicinterp');
a = min(X);
b = max(X);
c = min(Y);
d = max(Y);
volume_under_fit = quad2d(fitobject,a,b,c,d)