Compute that portion of a surface above a threshold
6 views (last 30 days)
Hi. I have a variable Z define over a grid of points X,Y. I would like to compute the area of my domain where Z>threshold.
Right now I tried to compute it as
A = trapz(X_vector, trapz(Y_vector, Z> threshold));
This gives acceptable results, however this approach set the border between the "in" and the "out" areas always in the middle of the last point outside and the first point inside.
To better explain the problem, consider the 1-dimensional case where
x = 1:10;
y = ones(1,10);
y(5) = 2;
threshold = 1.99;
plot(xlim,[1 1]*threshold ,'--k')
Analytically, the portion of the curve above the threshold is equal to 0.2 x-units. However if we plot
and we compute
L = trapz(x,y>threshold)
This is due to the fact that this approach ignores that y(5) is way closer to the threshold than y(4) and y(6).
On the 2D case You can visualize the "correct" solution using
The problem is that if the contour line is open the computation of the area becomes tricky.
Ahmet Cecen on 19 Dec 2016
The problem here is that you appear to be committing a fundamental fallacy in how discrete data work. In your example case, you are plotting this:
and claiming the "real" line segment under the threshold is 0.2 units long. However what you are actually inputting as data looks more like this:
in which case you can clearly see why MATLAB thinks the line segment is 1 units long. You are treating the graphical consequence of MATLAB's plotting by connecting the dots with lines as if it is the actual nature of the data.
I am assuming you are here because finding the analytical definition and integrating for your actual problem is a headache or not possible. I would suggest one of the following:
1) If you have a rough idea about the underlying nature of the surface, you can try to fit an analytical form to it, then find the area above the threshold via integration. For this you should be able to get away with `cftool` function.
2) If your data is just plain ugly, you can use an interpolant to estimate the area. Assuming your data is gridded like your example, use `griddedinterpolant` (`scatteredinterpolant` for ungridded) with an interpolation scheme appropriate for your case. You have treated the above example as if the data points are connected by straight lines, so you can choose `linear` (default) for that. Now you sample your data much denser, for the above example -> `Denser = F(4:0.01:6)`. F is the interpolant, and the step 0.01 you can change based on your accuracy needs and computation power. Now you can threshold `Denser` and use trapz, simple pixel count, a triangulation method or something like marching squares/cubes to get the area.