How to extract areas of contourf-plots below a given level

Hello,

i try to extract the data of the areas of a contourf-Plot below a certain level.

My approach:

1) get contourlines of the given level

2) merge Datas of contourlines and Z-Data

3) get binary Matrix of Z-Data

4) Find all areas with bwboundaries

L =

I have the Data of the red, green and blue crosses, but i have no idea how to fill the areas surrounded by the green crosses.

5) save these areas for later operations (e.g. as polygons)

Furthermore i would like to save these areas. Can someone help me?

Maybe theres something wrong with my approach?

Thanks a lot!

My Code:

% Data-Example
X = [8, 12, 16, 20, 24, 28, 32];
Y = [8, 20, 32, 44, 56, 68, 80];
Z1 = [0.2795, 0.2799, 0.2814, 0.2817, 0.2800, 0.2809, 0.2808;
    0.2786, 0.3040, 0.3093, 0.3113, 0.3130, 0.3131, 0.3118;
    0.2824, 0.2908, 0.2975, 0.2997, 0.2992, 0.2922, 0.2930;
    0.2864, 0.2847, 0.2871, 0.2798, 0.2753, 0.2714, 0.2686;
    0.2850, 0.2776, 0.2802, 0.2699, 0.2684, 0.2680, 0.2681;
    0.2851, 0.2830, 0.2732, 0.2684, 0.2678, 0.2684, 0.2688;
    0.2857, 0.2824, 0.2747, 0.2688, 0.2668, 0.2690, 0.2687];

% Given contourlevel
contourLevel = 0.28;
% Get contourlines
contourData = contourc(X, Y, Z1, [contourLevel contourLevel]);

% Get X-Y-Data of contourlines
k = 1;
allContourPoints = [];
while k < size(contourData, 2)
    level = contourData(1, k);
    numPoints = contourData(2, k);
    points = contourData(:, k+1:k+numPoints);

    if level == contourLevel
        allContourPoints = [allContourPoints; points(1, :)', points(2, :)'];
    end

    k = k + numPoints + 1;
end

% Remove Duplicates
allContourPoints = unique(allContourPoints, 'rows'); 

% Extend X and Y (in preparation to merge contourline-Data and Z1-Data)
extendedX = sort(unique([X(:); allContourPoints(:, 1)]));
extendedY = sort(unique([Y(:); allContourPoints(:, 2)]));
[gridX, gridY] = meshgrid(extendedX,extendedY);
% Extend Z1 by interpolation
extendedZ1 = griddata(X, Y, Z1, gridX, gridY, 'linear');
% Make sure that Z-Data of contourlines are = contourlevel
for i = 1:size(allContourPoints, 1)
    [~, closestXIdx] = min(abs(gridX(1, :) - allContourPoints(i, 1)));
    [~, closestYIdx] = min(abs(gridY(:, 1) - allContourPoints(i, 2)));
    disp(['Original:',num2str(extendedZ1(closestYIdx, closestXIdx))])
    extendedZ1(closestYIdx, closestXIdx) = contourLevel;
end

% Creat Binary MAtrix of extendedZ1
extendedBinaryZ1 = extendedZ1 <= contourLevel;

% Find connected areas in extendedBinaryZ1
[B, L] = bwboundaries(extendedBinaryZ1, 'noholes');

% Show contourf and Boundary-Data
figure
% plot contourf at contourlevel
contourf(X,Y,Z1,[contourLevel,contourLevel])
hold on
% plot Boundary-Data
for k = 1:length(B)
    b = B{k};

    x_values = extendedX(b(:, 2));
    y_values = extendedY(b(:, 1));

    switch k
        case 1
            x_color = 'rX';
        case 2
            x_color = 'gX';
        case 3
            x_color = 'bX';
        otherwise
            x_color = 'mX';
    end
    plot(x_values, y_values, x_color, 'MarkerSize', 10, 'LineWidth', 2);
end

Answer 1

[Star Strider]

It is a bit difficult to follow your code.

In my additions, these lines:

idx = find(contourData(1,:) == contourLevel);

clen = contourData(2,idx);

find and return the indices of the starting colums of the contours, and the lengths of each contour (in index units), corresponding to the contour matrix conventions. The code then uses contourf too fill the middle patch, and a loop:

axis([min(x1) max(x2) min(y1) max(y2)])

cm = 'rb';

for k = 1:numel(idx)

idxrng = (1:clen(k)-1)+idx(k);

yval = min(y1)*(k==1) + max(y2)*(k==2);

patch([contourData(1,idxrng), flip(contourData(1,idxrng))], [contourData(2,idxrng), ones(size(contourData(2,idxrng)))*yval],cm(k))

end

to fill the others.

Try this —

% Data-Example
X = [8, 12, 16, 20, 24, 28, 32];
Y = [8, 20, 32, 44, 56, 68, 80];
Z1 = [0.2795, 0.2799, 0.2814, 0.2817, 0.2800, 0.2809, 0.2808;
    0.2786, 0.3040, 0.3093, 0.3113, 0.3130, 0.3131, 0.3118;
    0.2824, 0.2908, 0.2975, 0.2997, 0.2992, 0.2922, 0.2930;
    0.2864, 0.2847, 0.2871, 0.2798, 0.2753, 0.2714, 0.2686;
    0.2850, 0.2776, 0.2802, 0.2699, 0.2684, 0.2680, 0.2681;
    0.2851, 0.2830, 0.2732, 0.2684, 0.2678, 0.2684, 0.2688;
    0.2857, 0.2824, 0.2747, 0.2688, 0.2668, 0.2690, 0.2687];



figure
% Given contourlevel
contourLevel = 0.28;
% Get contourlines
contourData = contourf(X, Y, Z1, [1 1]*contourLevel);
colormap(eye(3));

idx = find(contourData(1,:) == contourLevel);
clen = contourData(2,idx);

% figure
hold on
for k = 1:numel(idx)
    idxrng = (1:clen(k)-1)+idx(k);
    % plot(contourData(1,idxrng), contourData(2,idxrng), '.-')
    [x1(k),x2(k)] = bounds(contourData(1,idxrng));
    [y1(k),y2(k)] = bounds(contourData(2,idxrng));
end

axis([min(x1) max(x2) min(y1) max(y2)])
cm = 'rb';
for k = 1:numel(idx)
    idxrng = (1:clen(k)-1)+idx(k);
    yval = min(y1)*(k==1) + max(y2)*(k==2);
    patch([contourData(1,idxrng), flip(contourData(1,idxrng))], [contourData(2,idxrng), ones(size(contourData(2,idxrng)))*yval],cm(k))
end
hold off

return

% % Get X-Y-Data of contourlines
% k = 1;
% allContourPoints = [];
% while k < size(contourData, 2)
%     level = contourData(1, k);
%     numPoints = contourData(2, k);
%     points = contourData(:, k+1:k+numPoints);
% 
%     if level == contourLevel
%         allContourPoints = [allContourPoints; points(1, :)', points(2, :)'];
%     end
% 
%     k = k + numPoints + 1;
% end
% 
% % Remove Duplicates
% allContourPoints = unique(allContourPoints, 'rows'); 
% 
% % Extend X and Y (in preparation to merge contourline-Data and Z1-Data)
% extendedX = sort(unique([X(:); allContourPoints(:, 1)]));
% extendedY = sort(unique([Y(:); allContourPoints(:, 2)]));
% [gridX, gridY] = meshgrid(extendedX,extendedY);
% % Extend Z1 by interpolation
% extendedZ1 = griddata(X, Y, Z1, gridX, gridY, 'linear');
% % Make sure that Z-Data of contourlines are = contourlevel
% for i = 1:size(allContourPoints, 1)
%     [~, closestXIdx] = min(abs(gridX(1, :) - allContourPoints(i, 1)));
%     [~, closestYIdx] = min(abs(gridY(:, 1) - allContourPoints(i, 2)));
%     disp(['Original:',num2str(extendedZ1(closestYIdx, closestXIdx))])
%     extendedZ1(closestYIdx, closestXIdx) = contourLevel;
% end
% 
% % Creat Binary MAtrix of extendedZ1
% extendedBinaryZ1 = extendedZ1 <= contourLevel;
% 
% % Find connected areas in extendedBinaryZ1
% [B, L] = bwboundaries(extendedBinaryZ1, 'noholes');
% 
% % Show contourf and Boundary-Data
% figure
% % plot contourf at contourlevel
% contourf(X,Y,Z1,[contourLevel,contourLevel])
% hold on
% % plot Boundary-Data
% for k = 1:length(B)
%     b = B{k};
% 
%     x_values = extendedX(b(:, 2));
%     y_values = extendedY(b(:, 1));
% 
%     switch k
%         case 1
%             x_color = 'rX';
%         case 2
%             x_color = 'gX';
%         case 3
%             x_color = 'bX';
%         otherwise
%             x_color = 'mX';
%     end
%     plot(x_values, y_values, x_color, 'MarkerSize', 10, 'LineWidth', 2);
% end

.