In MATLAB, how to shade a small non overlapping area when two plots and a threshold line involved, as shown in figure attached

Question

I need to shade the area in Matlab plots, above the threshold, the small area as highlighted, under the blue curve (DC) but outside red one (PV), above the threshold; th=120

load('LP_DCHEMS_PV.mat','DC'); % DCCHEMS loaded
load('FifteenPercentHotDay.mat','PV') % Hot Day Case I, 15%
th=120;

plot(DC);
hold on
plot(PV);
yline(th);

When I use the following commands, it shades the area as shown in fig 2 attached. something has to be corrected in this line : h2=area([y1(:) , (PV(:)-DC(:)).* (DC(:)>th )]); If anyone can please guide,

y1=zeros(1,144);
y1(1,:)=th;
x=1:144;

    h2=area([y1(:) ,  (PV(:)-DC(:)).* (DC(:)>th )]);

    h2(1).FaceColor=[1 1 1];  % white
    h2(2).FaceColor=[0 0 1 ]; % blue
    hold on;
    plot(x,DC,'b',x,PV,'r');
    title('overlapped area ');
    
    xlim([1 144]);

Answer 1

Consider the script below:

clear; close all; clc;

% Generate dummy data
x  = 1:150;
y1 = exp(-(x(:)-35).^2/750)+1.4*exp(-(x(:)-100).^2/1000)+0.1*rand(150,1);
y2 = -1e-3*(x(:)-75).^2 + 1.85 + 0.1*rand(150,1); y2(y2<0) = 0;

% Set some threshold
threshold = 0.85;

% Create the plot
fig = figure(1); hold on; box on;

plot(x, y1, 'LineWidth', 1.5);
plot(x, y2, 'LineWidth', 1.5);
plot([x(1) x(end)], threshold*ones(1,2), 'k:', 'LineWidth', 1.5);

% Mark the areas and customise them (by changing FaceAlpha, as an example)
areas = MarkAreasInPlot(x, y1, y2, threshold, 'g');
for i = 1:length(areas)
    areas{i}.FaceAlpha = 0.5;
end

legend('y_1', 'y_2', 'threshold', 'Location', 'NorthWest');

At the top, I generate some dummy data. The generation of the y1 values is loosely based on what bla did in this answer, while y2 is simply a parabolic function (with some distortion) where all values below zero are filtered out. The next thing we do is define some threshold, which indicates the lowest possible boundary. After this, we simply plot the data and the threshold.

Then, I call the MarkAreasInPlot function:

function areas = MarkAreasInPlot(x, y1, y2, threshold, color)
% This function fills the areas for which the following two statements
% are true:
%   (1) y1 > y2
%   (2) y1 > threshold

    % Reshape arrays (to column vectors) if necessary
    x  = reshape(x, 1, []);
    y1 = reshape(y1, 1, []);
    y2 = reshape(y2, 1, []);
    
    % Create filter and find indices
    filter = (y1 > y2) & (y1 > threshold);
    idx = find(filter==1);
    
    % If nothing was found, return empty array
    if isempty(idx)
        areas = {};
        return;
    end

    % Determine the number of areas by looping over the indices (idx) and
    % checking if the next index follows up the previous one. If this is
    % not the case, it means that we are dealing with separate areas
    num_areas = 0;
    idx_areas = idx(1);
    
    for i = 1:length(idx)-1
        if(idx(i+1) ~= idx(i)+1)
            num_areas = num_areas + 1;
            idx_areas = [idx_areas, idx(i), idx(i+1)];
        end
    end
    idx_areas = [idx_areas, idx(end)];
    
    % Draw all areas
    areas = cell(num_areas+1, 1);
    
    for i = 1:2:2*(num_areas+1)
        % Determine span of indices for area
        idx_span = idx_areas(i):idx_areas(i+1);
        
        % Determine x- and y-coordinates
        x_area = x(idx_span);
        x_area = [x_area, flip(x_area)];
        
        y_top = y1(idx_span);
        y_bot = y2(idx_span);
        y_bot(y_bot < threshold) = threshold;
        
        y_area = [y_top, flip(y_bot)];
        
        % Fill the area
        areas{i} = fill(x_area, y_area, color, 'EdgeColor', 'None');
    end
end

This function works as follows. First, we make sure that all data arrays (x, y1 and y2) are column vectors (this is to prevent possible issues in case they are row vectors). Subsequently, we determine at which indices of the data arrays y1 > y2 and y1 > threshold. Note that filter is a logical array where an element is 1 only when the aforementioned two statements are true. Then, we determine the indices of the elements in filter were we have a value of 1. At these indices of the data arrays, we know that the two statements are true. Of course, we check for the case that the statements are false at all positions, in which case we don't draw any areas and return an empty cell array. After finding the indices, we loop over them and check for a case were the element at the next index is not equal to the value at the current index plus 1. If we find such a case, we now that we are dealing with two separate areas. Any time this happens, we store the two values at these two indices in idx_areas and increase the counter for the number of areas (num_areas) by 1. After doing this, idx_areas contains an even amount of elements. We know that we need to fill the area between y1 and max(y2, threshold) at the indices given by the range idx_areas(n):idx_areas(n+1) for n = 1, 3, 5, 7, ... (the values for n of course depending on the value of num_areas). This is exactly what we do at the end of the function, using the fill function. Note that we use a for loop over the number of areas we should fill and store the handle for each fill call in the areas cell array, which is returned.

If I run the script at the top, I get the following output: As you can see, it properly detects and fills the areas. However, due to the discrete nature of the data, it can happen that there is a small part at the side of an area that ends up not being completely filled: This is not an error! Instead, this is due to the fact that the next data point for y1 is located below the threshold, meaning that the second requirement for the MarkAreasInPlot function does not hold! You could try to 'solve' this using an interpolation procedure (to find at what x-value the y1-value is exactly equal to the threshold), but I think it is not worth the effort (as it is only a very small piece that is not filled).

To show you some more results, I now set threshold to 1.6 such that no area will be detected:

And as a final example, I swap y1 and y2 in the call to MarkAreasInPlot:

areas = MarkAreasInPlot(x, y2, y1, threshold, 'g');

As expected, it now fills the area where y2 > y1 and y2 > threshold: