Finding 2D area defined by contour lines in Matlab

Question

I am having difficulty with calculating 2D area of contours produced from a Kernel Density Estimation (KDE) in Matlab. I have three variables:

X and Y = meshgrid which variable 'density' is computed over (256x256) density = density computed from the KDE (256x256)

I run the code

contour(X,Y,density,10)

This produces the plot that is attached. For each of the 10 contour levels I would like to calculate the area. I have done this in some other platforms such as R but am having trouble figuring out the correct method / syntax in Matlab.

C = contourc(density) 

I believe the above line would store all of the values of the contours allowing me to calculate the areas but I do not fully understand how these values are stored nor how to get them properly.

Answer 1

This little script will help you. Its general for contour. Probably working for contour3 and contourf as well, with adjustments of course.

[X,Y,Z] = peaks;   %example data
% specify certain levels
clevels = [1 2 3];
C = contour(X,Y,Z,clevels);
xdata = C(1,:);   %not really useful, in most cases delimters are not clear
ydata = C(2,:);   %therefore further steps to determine the actual curves:

%find curves
n(1) = 1;         %n: indices where the certain curves start
d(1) = ydata(1);  %d: distance to the next index
ii = 1;
while true

       n(ii+1) = n(ii)+d(ii)+1;    %calculate index of next startpoint

       if n(ii+1) > numel(xdata)   %breaking condition
           n(end) = [];            %delete breaking point
           break
       end

       d(ii+1) = ydata(n(ii+1));   %get next distance
       ii = ii+1; 
end

%which contourlevel to calculate?
value = 2;             %must be member of clevels
sel = find(ismember(xdata(n),value));
idx = n(sel);          %indices belonging to choice
L = ydata( n(sel) );   %length of curve array

% calculate area and plot all contours of the same level
for ii = 1:numel(idx)
    x{ii} = xdata(idx(ii)+1:idx(ii)+L(ii));
    y{ii} = ydata(idx(ii)+1:idx(ii)+L(ii));

    figure(ii)
    patch(x{ii},y{ii},'red');            %just for displaying purposes
    %partial areas of all contours of the same plot
    areas(ii) = polyarea(x{ii},y{ii});   

end

% calculate total area of all contours of same level
totalarea = sum(areas)

Example: peaks (by Matlab)

Level value=2 are the green contours, the first loop gets all contour lines and the second loop calculates the area of all green polygons. Finally sum it up.

---

If you want to get all total areas of all levels I'd rather write some little functions, than using another loop. You could also consider, to plot just the level you want for each calculation. This way the contourmatrix would be much easier and you could simplify the process. If you don't have multiple shapes, I'd just specify the level with a scalar and use contour to get C for only this level, delete the first value of xdata and ydata and directly calculate the area with polyarea