How can I generate this plot?

Hi,

I am trying to generate similar looking plots to the one attached. With this, do you know what kind of plot this is? And how to get the black bars showing where the orifice opening is?

[Adam Danz]

The heat map can be created with imagesc().

It looks like there may be small arrows along a grid. That can be done with quiver().

The black lines are easy to plot with the plot() or line() functions.

The colorbar can be added with colorbar().

Check out the documentation for each of those functions for examples to get you started.

[Monique Embury]

Thank you

[Adam Danz]

Sure thing. If you get stuck, show us the code you've got and explain where you're stuck. I'd be happy to help as needed.

[Monique Embury]

Hi Adam,

Going back to the code I shared with you on another post, I am trying to plot the heat map with the quiver function overlaid. I am a bit confused on how to set up imagesc() after reading through the different tutorials. Could you help?

function [x,y,u_avg,v_avg,CHC_tot] = piv_averages(prefix,suffix,Nstart,Nfinish,interp)

% This program reads in a series of instantaneous PIV vector fields from
% Insight and averages them.  The user has the option of excluding
% interpolated vectors, which have CHC > 1.  (interp = 0 means do not 
% interpolate, while interp = 1 means interpolate).

% Create file name for each image
c_exit=2.776; %speed of sound
x0=1040; %origin of the jet in pixels  
y0=53.8; %origin of the jetin pixels
D=923.71; %diameter of jet exit in pixels 
v_shift=0;
for i = Nstart:Nfinish
    Nstring = int2str(i);       % Convert iteration number to a character string
    if i < 10
        filename_inst = strcat(prefix,'0000',Nstring,suffix);
    elseif i < 100
        filename_inst = strcat(prefix,'000',Nstring,suffix);
    elseif i < 1000
        filename_inst = strcat(prefix,'00',Nstring,suffix);
    elseif i < 10000
        filename_inst = strcat(prefix,'0',Nstring,suffix);
    else
        filename_inst = strcat(prefix,Nstring,suffix);
    end

    % Read file name
    A_inst = csvread(filename_inst,1,0);

    x = A_inst(:,1)-32;             % x-position (mm)
    y = A_inst(:,2)-198.2;             % y-position (mm)
    u = A_inst(:,3);             % x-velocity (m/s)
    v = A_inst(:,4);             % y-velocity (m/s)
    chc = A_inst(:,5);           % number of good vectors at this location

    N = size(x,1);               % Length of entire vector array
    
    % Initialize output variables if this is the first file
    if i == Nstart
        u_tot = zeros(N,1);
        v_tot = zeros(N,1);
        CHC_tot = zeros(N,1);
    end
     
    for j = 1:N
        if interp == 0
            if chc(j,1) == 1
                u_tot(j,1) = u_tot(j,1) + u(j,1);
                v_tot(j,1) = v_tot(j,1) + v(j,1);
                CHC_tot(j,1) = CHC_tot(j,1) + 1;
            end
        elseif interp == 1
            if chc(j,1) > 0
                u_tot(j,1) = u_tot(j,1) + u(j,1);
                v_tot(j,1) = v_tot(j,1) + v(j,1);
                CHC_tot(j,1) = CHC_tot(j,1) + 1;
            end
        end
    end
end

for j = 1:N
    u_avg(j,1) = u_tot(j,1)/CHC_tot(j,1);
    v_avg(j,1) = v_tot(j,1)/CHC_tot(j,1);
end

% define orifice opening in pixels 
orificepl=linspace(0, x0-(D/2));
orificepr=linspace(x0+(D/2), 2016);
yorifice=2016*ones(size(orificepl));

% define orifice opening in mm
orificeml=linspace(0, (x0-(D/2))/48.11);
orificemr=linspace((x0+(D/2))/48.11, 2016/48.11);

% define orifice opening non-dimensional
orificendl=(orificepl-x0)/D;
orificendr=(orificepr-x0)/D;

% Set origin to jet exit centerline
x_c = x - (x0);
y_c = y - y0;

% Shift by convective velocity
v = v - v_shift;

% Nondimensionalize variables
x_non = x_c/D;          % Nondimensionalize using jet diameter
y_non = y_c/D;          % Nondimensionalize using jet diameter
u_non = u_avg/c_exit;       % Nondimensionalize using sonic speed
v_non = v_avg/c_exit;       % Nondimensionalize using sonic speed


% Plot nondimensional vector field
figure(2)
quiver(x_non,y_non,u_non,v_non,4)
line(orificendl,yorifice/1008, 'linewidth', 4, 'color', 'k')
line(orificendr,yorifice/1008, 'linewidth', 4, 'color', 'k')
axis([-1 1 0 2])
xticks([-1 0 1])
yticks([0 1 2])
colormap(jet)
set(gca,'YtickLabel',2:-1:0)
set(gca,'XAxisLocation','top','YAxisLocation','left');
xlabel('z/D')
ylabel('x/D')
title('Nondimensional velocity field')
camroll(90)

[Adam Danz]

Hi Monique,

How can I help with the imagesc implementation? It's not clear where you're stuck.

Perhaps running your code could help me out and speed up the process but first I need the inputs and necessary files. Could you save the inputs and files in zip file and attach it?

Alternatively, perhaps the problem you're having can be reduced to a simpler example.

[Monique Embury]

I am stuck on how to take my (x,y,u,v) and make the correct matrix for the imagesc command.

Perhaps we can use this:

x=[10 9 8 7 6]

y=[5 4 3 2 1]

u=[0 0.5 1 1.5 2]

v=[0 0.4 0.8 1.2 1.6]

How would I make the imagesc() command work here with a quiver overlay?

[Adam Danz]

Is that what your input variables look like? Vectors (not matricies)? The image you shared contains quiver arrows along a grid and usually gridded data are matrices. Sometimes they can be vectors, too, but it's important to describe your variable correctly so I want to make sure your data are indeed vectors. If they are vectors, do they have repeated values (again, as a grid would)?

[Monique Embury]

My data is in a matrix. An avergae output matrix is as listed: https://drive.google.com/open?id=1Lk1zYrgMKFVlWNVuELOkrU85PJ5HN7xF

[Adam Danz]

I'm guessing the 1st and 2nd cols are the x and y vals and the 3rd and 4th cols are the u and v value and I am to ignore the 5th col. Correct? Or is the 5th column the heatmap values?

output_avg(1:10,:)
ans =
          232         2016          NaN          NaN            0
          248         2016      0.20987      0.24137           21
          264         2016      0.31538      0.47129           25
          280         2016      0.95483      0.66786           18
          296         2016      0.56917       1.5054           20
          312         2016      0.54912      0.77364           18
          328         2016       0.5497       1.1534           17
          344         2016      0.59935       1.1871           34
          360         2016      0.90434       1.3109           46
          376         2016      0.58487      0.78972           53

[Monique Embury]

5th col. is number of good vectors and can be ignored.

[Adam Danz]

What defines the color mapping? Another way of asking that, what defines the colobar scale?

Is it the magnitude of the quiver vectors? What defines positive and negatve values? If pos & neg values are defined by the direciton of the quiver vectors, which directions define pos and neg?

[Monique Embury]

I would say that it is defined by the magnitude. Where positive is defined from moving left to right and negitive by right to left.

Without applying the camera roll as I did in the figure, the positive would be from top to bottom and negative would be bottom to top.

Answer 1

[Adam Danz]

The quiver arrows all have fairly small magnitudes which makes it difficult to see the arrows but you can reproduce the figure and zoom in inspect parts of the image.

See inline comments for details.

% Load the data
data = load('average_vels.mat');

% Extract quiver variable for readability
x = data.output_avg(:,1); 
y = data.output_avg(:,2); 
u = data.output_avg(:,3); 
v = data.output_avg(:,4); 

% Compute signed magnitude of the quiver vectors where 
% positive magnitudes point upward and neg mags point downward
z = hypot(u,v) .* sign(v);

% Reshape the z vector into a matrix where 
% zMat(i,j) is the value at x=xVec(i) and y=yVec(j).
[xVec,~, zRow] = unique(x); 
[yVec,~, zCol] = unique(y);
zMat = nan([numel(xVec),numel(yVec)]); 
zIdx = sub2ind(size(zMat),zCol, zRow); 
zMat(zIdx) = z;

% Plot the "heatmap" first
clf()
imagesc(xVec, yVec, zMat)
set(gca, 'YDir', 'Normal') 
cbh = colorbar();
ylabel(cbh, 'Vector magnitude (positive is upward')
axis equal

% add the quiver data
hold on
quiver(x,y,u,v)

[Monique Embury]

Thank you very much!

[Adam Danz]

Glad I could help.

PS - don't use camroll(). Intead, just switch around the x&y, u&v values as needed if you want those values associated with a different axis.

[Monique Embury]

Thank you for the tip!

[Adam Danz]

Monique, one last thought. If you want white to represent 0 velocity and two other colors (red & blue in this example) to represent upward and downward quiver vectors, you can apply this colormap at the end of the code.

cmap = [linspace(0,1,128)'*[1,1], ones(128,1)]; % The blue range
cmap = [cmap; rot90(cmap,2)];                   % add the red range
set(gca, 'Colormap', cmap)                      % set the colormap
caxis(max(abs(cbh.Limits)) * [-1,1]);  %center the color range so white is in the middle