Extract the information of inside and outside of a contourf

Hello everyone,

I have a single line of code, quite simple, thanks to which I get (visually) exactly what I need:

h=figure; [pippo1, pippo2]=contourf(AZ,EL,FF, [0, 0]);

where the input arguments have size 301x301. The plot I get is like this:

As you can see, the outside of the contour is not only the big cyan area, but also some small areas in the white (not filled) polygon. I need the coordinates of the outside polygon (the whole cyan region) but I couldn't extract them from either "h" or "pippo2". Please note, I DON'T NEED the coordinates of the single contours, because if I extract XDATA and YDATA they lose the information about "the inside" and the "outside". I need to extract, in some way, a single contour that represents the whole cyan area, and another one that represents the complementary (white) one.

Thanks a lot!

Answer 1

[Garmit Pant]

Hi Roberto

To extract the regions based on the output of the “contourf” function, you can use the contour matrix output of the “contour” function.

“contour” outputs a contour matrix that defines the contour lines. It is of the following form:

You can extract the line data from the matrix and then use “inpolygon” function to create binary masks to extract the regions inside and out the contour lines. The following code snippet demonstrates how to extract regions inside the contour lines:

% Generate sample data
[X, Y, Z] = peaks;

% Generate contour data
M = contourf(X, Y, Z, [0, 0]);

% Extract regions from contour data
r1 = M(:, 2:M(2, 1) + 1);
r2 = M(:, M(2, 1) + 1 + 2:M(2, 1) + 1 + 2 + M(2, M(2, 1) + 1 + 1) - 1);
r3 = M(:, M(2, 1) + 1 + 2 + M(2, M(2, 1) + 1 + 1) + 1:end);

% Plot the regions
figure;
plot(r1(1, :), r1(2, :), 'o');
hold on;
plot(r2(1, :), r2(2, :), 'o');
plot(r3(1, :), r3(2, :), 'o');
xlim([-3 3]);
ylim([-3 3]);
hold off;

% Create masks for each region
in_mask1 = inpolygon(X, Y, r1(1, :), r1(2, :));
in_mask2 = inpolygon(X, Y, r2(1, :), r2(2, :));
in_mask3 = inpolygon(X, Y, r3(1, :), r3(2, :));

% Combine the masks
in_maskmain = in_mask1 | in_mask2 | in_mask3;

% Visualize the combined mask
figure;
surf(X, Y, int8(in_maskmain));
title('Combined Mask');
xlabel('X');
ylabel('Y');
zlabel('Mask');

% Initialize the masked Z matrix with zeros (or another value like NaN)
Z_masked = zeros(size(Z)); % or Z_masked = NaN(size(Z));

% Assign the heights from Z within the region to Z_masked
Z_masked(in_maskmain) = Z(in_maskmain);

% Visualize the masked Z data
figure;
M2 = contourf(X, Y, Z_masked);
title('Masked Z Data');
xlabel('X');
ylabel('Y');

For further understanding, kindly refer to the following MathWorks Documentation:

• Refer to the ‘Output Arguments’ section to understand the outputs generated: https://www.mathworks.com/help/releases/R2023b/matlab/ref/contourf.html
• Refer to the ‘Examples’ section: https://in.mathworks.com/help/releases/R2023b/matlab/ref/inpolygon.html

I hope you find the above explanation and suggestions useful!

[Steve]

Hi Garmit,

thanks for this useful solution! I have a pretty similar problem to solve.

Unfortunately, the inpolygon-function seems to "forget" about the edges ofthe plot. So the bottom right area isn't correctly masked. Like polyshape it seems to connect the first and the last point of the contourline.

Do you know a nice solution to fix that?

Thanks a lot!

[DGM]

If the goal is to create a mask, why are we even using contour()?

% Generate sample data
[X, Y, Z] = peaks;

% there's the mask.
mask = Z > 0;

There's no need for a second mask. Just use ~mask.

[Steve]

Hello,

thank you for your reply.

I have the goal to show the contourf-plot and afterwards i want to intersect it with other contourf-plots. (the part below a certain level)

Using masks is not giving me the same results, because i dont know exactly how contourf is interpolating and smoothing Data.

One example - Code:

% Give Data
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
level = 0.28;
[C,h] = contourf(X,Y,Z1,[level level]);
hold on
% Interpolation
[Xq, Yq] = meshgrid(8:0.05:32, 8:0.05:80); % Feinere Auflösung
% Zq = interp2(X, Y, Z1, Xq, Yq, 'linear'); 
% Zq = interp2(X, Y, Z1, Xq, Yq, 'cubic'); 
% Zq = interp2(X, Y, Z1, Xq, Yq, 'makima'); 
% Zq = interp2(X, Y, Z1, Xq, Yq, 'spline'); 
Zq = griddata(X, Y, Z1, Xq, Yq, 'linear'); 
% Zq = griddata(X, Y, Z1, Xq, Yq, 'cubic'); 
% Zq = griddata(X, Y, Z1, Xq, Yq, 'v4'); 
% Zq = griddata(X, Y, Z1, Xq, Yq, 'natural'); 
mask_Zq = Zq <= level;
% Indizes der 1en finden
[row, col] = find(mask_Zq == 1);
% Plot mask_Zq
figure;
hold on;
scatter(Xq(mask_Zq), Yq(mask_Zq), 10, 'k', 'filled');
contour(X,Y,Z1,[level level])
hold on
% Adjust Axes
xlim([min(min(X)),max(max(X))])
ylim([min(min(Y)),max(max(Y))])
% Find separate Areas
[B, L] = bwboundaries(mask_Zq, 'noholes');
% Plot separate Areas as Polygons
figure;
hold on;
for k = 1:length(B)
    boundary = B{k};
    fill(boundary(:,2), boundary(:,1), rand(1,3), 'FaceAlpha', 0.5);
end
hold off;