Make a plot with gradient shaded confidence intervals

Dear all,

I would like to fill the red shaded area in the plot below using a gradient fill (ideally with a jet colormap), so that the area between the blue line and the upper red dashed line goes from red to blue and the area between the blue line and the lower red dashed line goes from red to blue.

x = Data(:,1);
y = Data(:,2);
y10 = Data(:,3);
y90 = Data(:,4);
figure
plot(x, y, 'b', 'LineWidth', 2);  
hold on;
plot(x, y10, '--r');  
plot(x, y90, '--r'); 
fill([x; flipud(x)], [y10; flipud(y90)], 'r',...
    'FaceAlpha', 0.1, 'EdgeColor', 'none');  
hold off;
legend('Median', '10th-90th Percentiles', 'Location', 'best');
grid on;

Just to give an idea, something like this (without color steps but with a gradient fill):

I have tried myself, but I can't really find a nice (and fast) solution.

I have attched the Data file.

Any help would be grately appreciated!

[Umar]

Hi Simone,

To achieve a gradient fill effect in the specified red shaded area of the plot, we can use a combination of surf and colormap functions in Matlab. Here's a step-by-step guide to implement this:

Create a Meshgrid: Generate a meshgrid that covers the area to be filled with a gradient. This meshgrid will define the coordinates for the gradient fill.

[X, Y] = meshgrid(x, linspace(min(y10), max(y90), 100));

Calculate Colors: Determine the colors for the gradient fill based on the jet colormap. You can adjust the colormap to suit your preferences.

colors = jet(100); % Use jet colormap with 100 colors

Plot the Gradient Fill: Use the surf function to plot the gradient fill on the meshgrid.

figure; surf(X, Y, zeros(size(X)), 'CData', colors, 'FaceColor', 'interp', 'EdgeColor', 'none'); view(2); % 2D view for a flat gradient fill

Adjust Plot Settings: Customize the plot appearance to match the original plot.

hold on; plot(x, y, 'b', 'LineWidth', 2); plot(x, y10, '--r'); plot(x, y90, '--r'); hold off; legend('Median', '10th-90th Percentiles', 'Location', 'best'); grid on;

By following these steps, you can create a gradient fill effect in the specified red shaded area of the plot, transitioning smoothly from red to blue based on the jet colormap. Feel free to adjust the colormap, shading, or other parameters to achieve the desired visual effect.

Answer 1

[DGM]

Look into how to configure patch() objects:

https://www.mathworks.com/help/matlab/ref/matlab.graphics.primitive.patch-properties.html

You could specify a specific color for each vertex with interpolation between:

% some fake data
load data.mat

% split the data array
x = Data(:,1);
y = Data(:,2);
y10 = Data(:,3);
y90 = Data(:,4);

% color endpoints
colorA = [1 0 0];
colorB = [0 0 1];
CData = [repmat(colorA,[numel(x),1]); repmat(colorB,[numel(x),1])];
alph = 0.5;

figure
hold on;
patch([x; flipud(x)], [y; flipud(y90)], permute(CData,[1 3 2]),...
    'FaceAlpha', alph, 'EdgeColor', 'none');  
patch([x; flipud(x)], [y; flipud(y10)], permute(CData,[1 3 2]),...
    'FaceAlpha', alph, 'EdgeColor', 'none');  
hp(1) = plot(x, y, 'b', 'LineWidth', 2);  
hp(2) = plot(x, y10, '--k', 'LineWidth', 2);  
plot(x, y90, '--k', 'LineWidth', 2); 
hold off;
legend(hp,'Median', '10th-90th Percentiles', 'Location', 'best');
grid on;

... or you could specify a linear value used as a relative position into the current colormap (with interpolation):

% some fake data
load data.mat

% split the data array
x = Data(:,1);
y = Data(:,2);
y10 = Data(:,3);
y90 = Data(:,4);

% color endpoints
CData = [zeros(numel(x),1); ones(numel(x),1)];
CT = jet(64);
alph = 0.5;

figure
hold on;
patch([x; flipud(x)], [y; flipud(y90)], CData,...
    'FaceAlpha', alph, 'EdgeColor', 'none');  
patch([x; flipud(x)], [y; flipud(y10)], CData,...
    'FaceAlpha', alph, 'EdgeColor', 'none');  
hp(1) = plot(x, y, 'b', 'LineWidth', 2);  
hp(2) = plot(x, y10, '--k', 'LineWidth', 2);  
plot(x, y90, '--k', 'LineWidth', 2); 
hold off;
legend(hp,'Median', '10th-90th Percentiles', 'Location', 'best');
grid on;

colormap(CT)
clim([0 1])

[Simone]

Hi @DGM this exactly what i was looking for, thanks a lot for it!

The only issue is if I try another dataset (find it attached), the output doesn't look nice and smooth as with the previous dataset. Do you know why?

Thanks a lot once again!

[DGM]

Oof. Maybe there's a way to wrangle simple interpolation to get that to work straight with patch(), but I'm not getting anywhere.

I'm just going to put my caveman hat on and do it with a raster image. I have no doubt that there are probably better ways to do this, and this example itself can probably be simplified.

% the example data
load Data2.mat
Data = Data2(1:100,:);

% split the data array
x = Data(:,1);
y = Data(:,2);
y10 = Data(:,3);
y90 = Data(:,4);

% color underlay parameters
CT = jet(64);
alph = 0.5;
reversegradient = false; 
sz = [1000 1000]; % size of the raster image

% the image space
xv = 1:sz(2);
yv = (1:sz(1)).';

% just rescale everything into image coordinates
xrange = imrange(x);
yrange = imrange([y10 y90]);
xim = rescale(x,xrange,[1 sz(2)]).';
y90im = rescale(y90,yrange,[1 sz(1)]); % lower
yim = rescale(y,yrange,[1 sz(1)]); % middle
y10im = rescale(y10,yrange,[1 sz(1)]); % upper

% interpolate onto the uniform image grid
y90im = interp1(xim,y90im,xv); % lower
yim = interp1(xim,yim,xv); % middle
y10im = interp1(xim,y10im,xv); % upper

% generate the underlay image
% this is just a piecewise line eqn with vector inputs
yv = repmat(yv,[1 sz(2)]);
mku = yv >= yim;
Z = ~mku.*(yv - y90im)./(yim - y90im) ...
    + mku.*(yv - y10im)./(yim - y10im);

% just plot the lines first so that xlim,ylim are defined
hold on;
hp(1) = plot(x, y, 'b', 'LineWidth', 2);  
hp(2) = plot(x, y10, '--k', 'LineWidth', 2);  
plot(x, y90, '--k', 'LineWidth', 2); 
legend(hp,'Median', '10th-90th Percentiles', 'Location', 'best');
xl = xlim; yl = ylim; % cache these

% show the image
% mask off all the external area and set opacity
hi = imagesc(xrange,yrange,Z);
hi.AlphaData = (Z >= 0)*alph;

% final plot setup
hold off;
grid on;
if reversegradient
    colormap(CT)
else
    colormap(flipud(CT))
end
clim([0 1])
xlim(xl); ylim(yl);

% purely out of spite for rescale()'s cumbersome syntax
function yy = rescale(xx,inrange,outrange)
    scalefactor = diff(outrange)/diff(inrange);
    yy = scalefactor*(double(xx) - inrange(1)) + outrange(1);
end

[Simone]

hi @DGM once again, thanks a lot for your time and for your help.

Your code is visually great, but unfortunatelly the colors get unconsistently distributed within the bounded region (see below).

As my code first computes the moving median of a given dataset (see code below), I need the colors to be representive of the that percentile. For instance, a given shade of green must be representative of the 70th percentile accross the whole figure. However, if I plot the actual 70th percentile, I can see that it touches upon different colors.

Can you think a way of obtaining the output I am looking for?

Thanks a lot!

Please find the original xy dataset attached, and see the code below:

%% My Code

clear
close all
clc

load('xy.mat')
x = xy(:,1);
y = xy(:,2);

windowSize = days(90);
prc_Min = 10;
prc_Max = 90;
prc_Stp = 0.1;
prc_Stps = prc_Min:prc_Stp:prc_Max;

Y_median = zeros(size(y));                                                 % Compute moving median and percentiles with variable time step
Y_10 = zeros(size(y));
Y_30 = zeros(size(y));
Y_70 = zeros(size(y));
Y_90 = zeros(size(y));

for i = 1:length(x)
    window_start = x(i);
    window_end = window_start + windowSize;
    
    idx = find(x >= window_start & x <= window_end);
    
    if ~isempty(idx)
        Y_median(i) = median(y(idx), 'omitnan');
        Y_10(i) = prctile(y(idx), prc_Min);
        Y_90(i) = prctile(y(idx), prc_Max);
        Y_70(i) = prctile(y(idx), 70);
        Y_30(i) = prctile(y(idx), 30);
    else
        Y_median(i) = NaN;
        Y_10(i) = NaN;
        Y_90(i) = NaN;
        Y_70(i) = NaN;
        Y_30(i) = NaN;
    end
end

figure;
plot(x, Y_median, 'b', 'LineWidth', 2);  
hold on;
plot(x, Y_10, '--r');  % 10th percentile
plot(x, Y_90, '--r');  % 90th percentile
fill([x; flipud(x)], [Y_10; flipud(Y_90)], 'r', 'FaceAlpha', 0.1, 'EdgeColor', 'none');  % Confidence interval fill
hold off;
legend('Median', '10th-90th Percentiles', 'Location', 'best');
grid on;



%% Your Code

keep x Y_median Y_10 Y_30 Y_70 Y_90 

% split the data array
y = Y_median;
y10 = Y_90; % Was set the other way around
y90 = Y_10; % Was set the other way around
y30 = Y_30;
y70 = Y_70;

% color underlay parameters
CT = jet(64);
alph = 0.5;
reversegradient = false; 
sz = [1000 1000]; % size of the raster image
% the image space
xv = 1:sz(2);
yv = (1:sz(1)).';
% just rescale everything into image coordinates
xrange = imrange(x);
yrange = imrange([y10 y90]);
xim = rescale(x,xrange,[1 sz(2)]).';
y90im = rescale(y90,yrange,[1 sz(1)]); % lower
yim = rescale(y,yrange,[1 sz(1)]); % middle
y10im = rescale(y10,yrange,[1 sz(1)]); % upper
% interpolate onto the uniform image grid
y90im = interp1(xim,y90im,xv); % lower
yim = interp1(xim,yim,xv); % middle
y10im = interp1(xim,y10im,xv); % upper
% generate the underlay image
% this is just a piecewise line eqn with vector inputs
yv = repmat(yv,[1 sz(2)]);
mku = yv >= yim;
Z = ~mku.*(yv - y90im)./(yim - y90im) ...
    + mku.*(yv - y10im)./(yim - y10im);
% just plot the lines first so that xlim,ylim are defined

figure
hold on;
hp(1) = plot(x, y, 'b', 'LineWidth', 2);  
hp(2) = plot(x, y10, '--k', 'LineWidth', 2);  
plot(x, y90, '--k', 'LineWidth', 2); 
plot(x, y90, '--k', 'LineWidth', 2); 
hp(3) = plot(x, y30, '-k', 'LineWidth', 2); 
plot(x, y70, '-k', 'LineWidth', 2); 
legend(hp,'Median', '10th-90th Percentiles','30th-70th Percentiles', 'Location', 'best');
xl = xlim; yl = ylim; % cache these
% show the image
% mask off all the external area and set opacity
hi = imagesc(xrange,yrange,Z);
hi.AlphaData = (Z >= 0)*alph;
% final plot setup
hold off;
grid on;
if reversegradient
    colormap(CT)
else
    colormap(flipud(CT))
end
clim([0 1])
xlim(xl); ylim(yl);

% purely out of spite for rescale()'s cumbersome syntax
function yy = rescale(xx,inrange,outrange)
scalefactor = diff(outrange)/diff(inrange);
yy = scalefactor*(double(xx) - inrange(1)) + outrange(1);
end
function [minval maxval] = imrange(inpict)
%   [min max] = IMRANGE(INPICT)
%       returns the global min and max pixel values in INPICT
%       
%   INPICT can be a vector or array of any dimension or class
x = inpict(:);
minval = double(min(x));
maxval = double(max(x));
if nargout < 2
    minval = cat(2,minval,maxval);
end
end

[DGM]

Ah. I'm just oversimplifying this whole thing, aren't I?

As you can see, I'm just mapping piecewise-linearly on y between the three curves. Obviously, that's not appropriate here. I didn't think about what the data meant. I was just slapping a fill in there.

I'm going to have to think about how to approach this. The solutions that come to mind only make the apparent inelegance of my prior example even worse. I don't know if I'll be able to finish this tonight. We'll see.

[DGM]

Sorry about not being at the computer. Try this. It's still a kludge on my part, but we can at least see if it's doing what you need. I'm just doing PWL interpolation for the gradient instead of trying to calculate the full distribution. Even with only 9 breakpoints, it's a challenge to notice any banding in the output, and the interpolation should be exact at the values in prc_Stps.

load('xy.mat')
x = xy(:,1);
y = xy(:,2);

windowSize = days(90);
prc_Min = 10; 
prc_Max = 90;

% idk if you were going to use this for something else
% but i'm going to use it
prc_Stp = 10;
prc_Stps = prc_Min:prc_Stp:prc_Max;


% Compute moving median and percentiles with variable time step
% don't really need Y_median so long as Y contains 50
npk = numel(prc_Stps);
Y_median = zeros(size(y)); 
Y = zeros(numel(y),npk);
for i = 1:length(x)
    window_start = x(i);
    window_end = window_start + windowSize;
    
    idx = x >= window_start & x <= window_end;
    
    if ~isempty(idx)
        Y_median(i) = median(y(idx), 'omitnan');
        for pk = 1:npk
            Y(i,pk) = prctile(y(idx), prc_Stps(pk));
        end
    else
        Y_median(i) = NaN;
        Y(i,pk) = NaN;
    end
end

figure
plot(x, Y_median, 'b', 'LineWidth', 2);  
hold on;
plot(x, Y(:,1), '--r');  % 10th percentile
plot(x, Y(:,end), '--r');  % 90th percentile
fill([x; flipud(x)], [Y(:,1); flipud(Y(:,end))], 'r', 'FaceAlpha', 0.1, 'EdgeColor', 'none');  % Confidence interval fill
hold off;
legend('Median', '10th-90th Percentiles', 'Location', 'best');
grid on;
xlim(imrange(x))

%% Your Code
% idk what keep() is, but i don't have it
%keep x Y_median Y_10 Y_30 Y_70 Y_90 

% split the data array
y = Y_median;

% color underlay parameters
CT = jet(64);
alph = 0.5;
reversegradient = false; 
sz = [1000 1000]; % size of the raster image

% the image space
xv = 1:sz(2);
yv = (1:sz(1)).';

% just rescale everything into image coordinates
xrange = imrange(x);
yrange = imrange(Y);
xim = rescale(x,xrange,[1 sz(2)+1]).';
Yim = rescale(Y,yrange,[1 sz(1)]); % lower
% interpolate onto the uniform image grid
Yim = interp1(xim,Yim,xv); % lower

% generate the underlay image
% this could probably also be done with interp2()
pctrange = [prc_Min prc_Max];
pctswing = max(abs(50 - pctrange));
z = 1 - abs(prc_Stps - (pctswing + prc_Min))/pctswing;
Z = zeros(sz);
for k = 1:sz(2)
    Z(:,k) = interp1(Yim(k,:),z,yv.','linear','extrap');
end

% just plot the lines first so that xlim,ylim are defined
figure
hold on;
hp(1) = plot(x, y, 'b', 'LineWidth', 2);  
hp(2) = plot(x, Y(:,1), '--k', 'LineWidth', 2);  
plot(x, Y(:,end), '--k', 'LineWidth', 2); 

hp(3) = plot(x, Y(:,3), '-k', 'LineWidth', 2); 
plot(x, Y(:,7), '-k', 'LineWidth', 2); 
legend(hp,'Median', '10th-90th Percentiles','30th-70th Percentiles', 'Location', 'best');
yl = ylim; % cache these

% show the image
% mask off all the external area and set opacity
hi = imagesc(xrange,yrange,Z);
hi.AlphaData = (Z >= 0)*alph;

% final plot setup
hold off;
grid on;
if reversegradient
    colormap(CT)
else
    colormap(flipud(CT))
end
clim([0 1])
xlim(imrange(x)); ylim(yl);

% purely out of spite for rescale()'s cumbersome syntax
function yy = rescale(xx,inrange,outrange)
    scalefactor = diff(outrange)/diff(inrange);
    yy = scalefactor*(double(xx) - inrange(1)) + outrange(1);
end

function [minval maxval] = imrange(inpict)
    %   [min max] = IMRANGE(INPICT)
    %       returns the global min and max pixel values in INPICT
    %       
    %   INPICT can be a vector or array of any dimension or class
    x = inpict(:);
    minval = double(min(x));
    maxval = double(max(x));
    if nargout < 2
        minval = cat(2,minval,maxval);
    end
end

[Simone]

Hi @DGM please don't be sorry. Your help has been essential (and so much appreciated).

Indeed, your code works great and does exactly what I need.

However, the only problem with imagesc is that if I want to plot it with a log scale, it won't work.

To overcome this issue, I came up with another piece of code that can deal with log scales (see code below for future reference).

That said, without your help it would have not been possible to solve my problem.

Thank you very much indeed

clear 
close all
clc

load('xy.mat')

x = xy(:,1);
y = xy(:,2);

windowSize = days(90);
prc_Min = 10;
prc_Max = 90;
prc_Stp = 1;

prc_Stps = prc_Min:prc_Stp:prc_Max;

npk = numel(prc_Stps);
Y_Median = zeros(size(y));
YY = zeros(numel(y),npk);

for i = 1:length(x)
    window_start = x(i);
    window_end = window_start + windowSize;
    idx = x >= window_start & x <= window_end;
    if ~isempty(idx)
        Y_Median(i) = median(y(idx), 'omitnan');
        for pk = 1:npk
            YY(i,pk) = prctile(y(idx), prc_Stps(pk));
        end
    else
        Y_Median(i) = NaN;
        YY(i,pk) = NaN;
    end
end

figure
plot(x, Y_Median, 'b', 'LineWidth', 2);
hold on;
plot(x, YY(:,1), '--r');    % 10th percentile
plot(x, YY(:,end), '--r');  % 90th percentile
fill([x; flipud(x)], [YY(:,1); flipud(YY(:,end))], 'r', 'FaceAlpha', 0.1, 'EdgeColor', 'none');  % Confidence interval fill
hold off;
legend('Median', '10th-90th Percentiles', 'Location', 'best');
grid on;

j_Up = find(YY(1,:)>Y_Median(1));
j_Down = find(YY(1,:)<Y_Median(1));
sz = (numel(YY(1,:))-1)/2;
CT_Down = jet(sz);
CT_Up = flipud(jet(sz));
FaceAlpha_Transp = 1;

figure
hold on

Patch_X = [x; flipud(x)];

k = 1;
for i_Up = 1 : numel(j_Up)
    Patch_Y_Up = [YY(:, j_Up(i_Up)); flipud(YY(:, j_Up(i_Up)-1))];
    fill(Patch_X, Patch_Y_Up, CT_Up(k,:), 'FaceAlpha',...
        FaceAlpha_Transp,'EdgeColor','none');
    k = k + 1;
end

k = 1;
for i_Down = 1 : numel(j_Down)
    Patch_Y_Down = [YY(:, j_Down(i_Down)); flipud(YY(:, j_Down(i_Down)+1))];
    fill(Patch_X, Patch_Y_Down, CT_Down(k,:), 'FaceAlpha',...
        FaceAlpha_Transp,'EdgeColor','none');
    k = k + 1;
end

hp(1) = plot(x, Y_Median, '-k', 'LineWidth', 2); % Median
hp(2) = plot(x, YY(:,21), '-k', 'LineWidth', 1); % 30th Percentile
hp(3) = plot(x, YY(:,61), '-k', 'LineWidth', 1); % 70th Percentile
legend(hp,'Median', '30th-70th Percentiles', 'Location', 'best');
grid on