According to my understanding, you want to calculate the surface area of a certain portion of your graph. I have made an assumption that you will be identifying this portion by limiting the z - values of the plot as you have stated taking Z < 0.01 as an example. In that case, we can divide the plot into triangles and then sum the areas of the smaller triangles to find the area of the plot.
To do this, I would suggest plotting your data in the form of a 3 - D triangular mesh and create a triangulation object to plot. This is because we can leverage some properties of the triangulation object to calculate the areas of the triangles that form the plot. I am sharing a simple example to demonstrate my point.
[X,Y] = meshgrid(-8:.05:8);
R = sqrt(X.^2 + Y.^2) + eps;
Z = sin(R)./R;
T = delaunay(X,Y);
TO = triangulation(T,X(:),Y(:),Z(:));
The triangulation object TO has two properties TO.Points and TO.ConnectivityList.
Triangulation connectivity list, specified as an m-by-3 matrix, where m is the number of triangles in the plot. Each row of T contains the vertex IDs that define a triangle.
Points, specified as a matrix whose columns are the x-coordinates, y-coordinates, and z-coordinates of the triangulation points. The row numbers of P are the vertex IDs in the triangulation.
The code below calculates the surface area of the graph of the function defined above but ignores those where the functional value is greater than 0.4 (you can replace this with your threshold value and function).
area = 0;
for i = 1:size(TO.ConnectivityList,1)
a = TO.Points(TO.ConnectivityList(i,:),:);
if any(a(:,3) > 0.4)
p1 = a(1,:);
p2 = a(2,:);
p3 = a(3,:);
area = area + 0.5 * norm(cross(p2-p1,p3-p1));
You may use these ideas to calculate the area of the graph that satisfies certain conditons on the function range.
See , and  to learn more about the concepts used above.