Get quiver plot arrows to connect up contour lines

I have the following code here that plots a contour plot for a given function and also plots a color-indicated directional field of the gradient on top of it. The arrows are all the same length but the difference of color indicates the magnitude. I am trying to do two things here:

1. Get the arrows to start at a contour and end at the next contour. So there should be no arrows "crossing" a contour line
2. Get the colormap of the contours to be different than the colormap of the directional field

I think in order to do (1) I need to get the coordinates of the contours themselves so to provide the starting points of the arrows, but I don't know how to do that. As for (2), I better there is clever way of doing this, but my only thought is to create two entirely different plots and draw them on top of each other which seems too brute force. Any help is appreciated, code is below:

clf

f = @(x,y) 0.5*sin(pi/2*(2*x-y)) + 2*y + 1

f = function_handle with value:
    @(x,y)0.5*sin(pi/2*(2*x-y))+2*y+1

[X,Y] = meshgrid(-1:0.05:1,-1:0.05:1);

Z = f(X*cos(pi/6)-Y*sin(pi/6),Y*cos(pi/6)+X*sin(pi/6));

% Create axes
ax = axes('Color', [1 1 1]); 
hold(ax,'on')

contour(ax, X,Y,Z,10)

% Demo data
scale_size = 0.25;

[X,Y] = meshgrid(-1:0.2:1,-1:0.2:1);

Z = f(X*cos(pi/6)-Y*sin(pi/6),Y*cos(pi/6)+X*sin(pi/6));
[U,V] = gradient(Z,0.2,0.2);

% Define the colormap for the quiver arrows.
% cmap can have any number of rows.
cmap = spring(255); 
ax.Colormap = cmap; 
% Assign colors based on magnitude of vectors
vectorMagnitude = hypot(U(:),V(:)); 
% Scale magnitudes to rows of colormap
vecMagNorm = (vectorMagnitude-min(vectorMagnitude))./range(vectorMagnitude);
vecColorIdx = round(vecMagNorm * (size(cmap,1)-1)) + 1; 
% Plot the quiver data
dU = U./sqrt(U.^2 + V.^2);
dV = V./sqrt(U.^2 + V.^2);
for i = 1:numel(Z)
    quiver(ax, X(i),Y(i),dU(i),dV(i), 0.1, ...
        'Color', cmap(vecColorIdx(i),:), 'LineWidth',0.2,'MaxHeadSize',2)
end
% Set properties for the main axes
axis equal
xlim(ax, [-1 1])
ylim(ax, [-1 1])
% Add colorbar
cb = colorbar(ax); 
% Set colorbar range
caxis(ax, [floor(vectorMagnitude(1)), ceil(vectorMagnitude(2))])
% Label the colorbars
ylabel(cb,'Vector magnitude')
hold off

EDIT:

I was able to produce something for the vectors being spaced exactly between contours pertaining to (1) which uses the File Exchange function intersections found here: Fast and Robust Curve Intersections. The code and plot are below

clf

f = @(x,y) 0.5*sin(pi/2*(2*x-y)) + 2*y + 1;
[X,Y] = meshgrid(-1:0.05:1,-1:0.05:1);

Z = f(X*cos(pi/6)-Y*sin(pi/6),Y*cos(pi/6)+X*sin(pi/6));

syms x y
dz = gradient(f(x*cos(pi/6)-y*sin(pi/6),y*cos(pi/6)+x*sin(pi/6)),[x,y]);
dz = matlabFunction(dz);

% Create axes
ax = axes('Color', [1 1 1]); 
hold(ax,'on')
NumOfContours = 10;
[M,~] = contour(ax, X,Y,Z,NumOfContours);

i = 1;
Ns = zeros(1,NumOfContours);
counter = 0;
while i<=length(M)
    counter = counter + 1;
    Ns(counter) = M(2,i);
    i = i + M(2,i)+1;
end

X = nan(1,1000); % Pre-allocation
Y = X;

Ms = zeros(1,NumOfContours-1);

startIdx = 2;
overallCounter = 1;
for j=1:NumOfContours-1
    subX = M(1,startIdx:(startIdx + Ns(j) - 1));
    subY = M(2,startIdx:(startIdx + Ns(j) - 1));
    startIdx = startIdx + Ns(j) + 1;

    dX = gradient(subX);
    dY = gradient(subY);

    if j==1
        PrimaryLength = sum(hypot(dX,dY));
        NumVecPerUnit = 4;
        dl = PrimaryLength/(NumVecPerUnit+1);
    end

    Lens = cumsum(hypot(dX,dY));

    counter = 1;
    gatePass = 1;
    while gatePass      
        [~,idx] = min(abs(dl*counter - Lens));
        if idx+1 > length(Lens)
            gatePass = 0;
            Ms(j) = counter-1;
        else
            counter = counter + 1;
            X(overallCounter) = subX(idx);
            Y(overallCounter) = subY(idx);
            overallCounter = overallCounter + 1;
        end
    end
end

X(isnan(X)) = [];
Y(isnan(Y)) = [];

U = zeros(1,length(X));
V = U;

for i=1:length(X)
    d = dz(X(i),Y(i));
    U(i) = d(1);
    V(i) = d(2);
end

% Define the colormap for the quiver arrows.
% cmap can have any number of rows.
cmap = parula(255);
ax.Colormap = cmap; 

scale_size = zeros(1,length(X)); % Pre-allocate scaling array

dU = U./sqrt(U.^2 + V.^2);
dV = V./sqrt(U.^2 + V.^2);

startIdx = Ns(1)+3;
j = 0;
for i=1:length(X)    
    j_gate = i > cumsum([0,Ms]);
    j_switch = find(j_gate == max(j_gate),1,"last");
    if j_switch > j
        j = j_switch;
        x1 = M(1,startIdx:(startIdx + Ns(j+1) - 1));
        y1 = M(2,startIdx:(startIdx + Ns(j+1) - 1));
        startIdx = startIdx + Ns(j+1) + 1;
    end

    x2 = [X(i),X(i)+sqrt(2)*dU(i)];
    y2 = [Y(i),Y(i)+sqrt(2)*dV(i)];

    [x0,y0] = intersections(x1,y1,x2,y2,1);

     false_switch = (isempty(x0) || isempty(y0));
     x0(false_switch) = inf;
     y0(false_switch) = inf;

    scale_size(i) = hypot(X(i)-x0(1),Y(i)-y0(1));
end
dU = scale_size.*dU;
dV = scale_size.*dV;

remove_id = false(1,length(X));

for i=1:length(X)
    if norm([X(i)+dU(i),Y(i)+dV(i)],"inf") > 1
        remove_id(i) = true;
    end
end

X(remove_id) = [];
Y(remove_id) = [];
dU(remove_id) = [];
dV(remove_id) = [];

quiver(ax, X+0.05*dU,Y+0.05*dV,0.925*dU,0.925*dV, ...
    'Color', [0.1 0.1 0.1], 'LineWidth',0.8,'MaxHeadSize',2,'AutoScale','off')

% scatter(X,Y)

% Set properties for the main axes
axis equal
xlim(ax, [-1 1])
ylim(ax, [-1 1])

hold off

Answer 1

[Cris LaPierre]

Getting the coordinates of your contour lines is fairly straightforward, as it's an output of the contour function. There's a little coding work needed to turn that output into something you can use, but doable.

Getting the quiver arrows to appear as you describe is not giong to be an easy task. Doable, but not pretty, as that is not what a quiver plot does. I by no means attempted to arriave at a final solution, but here is a working example (simplified)

f = @(x,y) 0.5*sin(pi/2*(2*x-y)) + 2*y + 1;
[X,Y] = meshgrid(-1:0.05:1,-1:0.05:1);

Z = f(X*cos(pi/6)-Y*sin(pi/6),Y*cos(pi/6)+X*sin(pi/6));

% Create axes
ax = axes('Color', [1 1 1]); 
hold(ax,'on')

%% Extract contour line coordinates
m = contour(ax, X,Y,Z,10);

i = 1;
lvl = [];
x = [];
y = [];
while i < length(m)
    clvl = m(1,i);
    num = m(2,i);
    
    lvl = [lvl;clvl];

    % extracted (x,y) coordinates
    tmpx = m(1,i+1:i+num);
    tmpy = m(2,i+1:i+num);

    % convert contour coordinates to 10 equally spaced coordinates
    newx = linspace(tmpx(1),tmpx(end),10);
    newy = interp1(tmpx,tmpy,newx);
    x = [x,newx'];
    y = [y,newy'];

    i = i+num+1;
end

% Added this to visualize the points
scatter(x,y,10,"magenta","filled")

%% Calculate U as difference of x coordinates and V as differeing of y coordinates
U = diff(x,[],2);
V = diff(y,[],2);

% Create quiver plot
quiver(ax,x(:,1:end-1),y(:,1:end-1),U,V,"off")

% Set properties for the main axes
axis equal
xlim(ax, [-1 1])
ylim(ax, [-1 1])
% Add colorbar
cb = colorbar(ax); 
hold off

As for colormap, that is an axiss property. The solution is to either manually set the color of one of your plots, or overlay 2 axes, and specify the colormap of each. Place your plots on the axes whose colormap you want it to use.

Here are some posts about this

• https://www.mathworks.com/matlabcentral/answers/194554-how-can-i-use-and-display-two-different-colormaps-on-the-same-figure
• https://www.mathworks.com/matlabcentral/answers/101346-how-do-i-use-multiple-colormaps-in-a-single-figure-in-r2014a-and-earlier

[Cris LaPierre]

Here's a different example that shows how you could use 2 different colormaps in a single figure. This is meant to address question 2, and does not make an attempt at incorporating a solution for your first question.

f = @(x,y) 0.5*sin(pi/2*(2*x-y)) + 2*y + 1;
[X,Y] = meshgrid(-1:0.05:1,-1:0.05:1);

Z = f(X*cos(pi/6)-Y*sin(pi/6),Y*cos(pi/6)+X*sin(pi/6));

% Create axes
t = tiledlayout(1,1);
ax = axes(t,'Color', [1 1 1]); 
contour(ax,X,Y,Z,10)

% Set properties for the main axes
axis equal
xlim(ax, [-1 1])
ylim(ax, [-1 1])

% Demo data
scale_size = 0.25;

[X,Y] = meshgrid(-1:0.2:1,-1:0.2:1);

Z = f(X*cos(pi/6)-Y*sin(pi/6),Y*cos(pi/6)+X*sin(pi/6));
[U,V] = gradient(Z,0.2,0.2);

% Define the colormap for the quiver arrows.
% cmap can have any number of rows.
ax2 = axes(t);
colormap(ax2,'spring'); 
cmap = colormap(ax2);
hold on

% Assign colors based on magnitude of vectors
vectorMagnitude = hypot(U(:),V(:)); 
% Scale magnitudes to rows of colormap
vecMagNorm = (vectorMagnitude-min(vectorMagnitude))./range(vectorMagnitude);
vecColorIdx = round(vecMagNorm * (size(cmap,1)-1)) + 1; 
% Plot the quiver data
dU = U./sqrt(U.^2 + V.^2);
dV = V./sqrt(U.^2 + V.^2);
for i = 1:numel(Z)
    quiver(ax2, X(i),Y(i),dU(i),dV(i), 0.1, ...
        'Color', cmap(vecColorIdx(i),:), 'LineWidth',0.2,'MaxHeadSize',2)
end


% Add colorbar
cb = colorbar(ax2); 
% Set colorbar range
caxis(ax2, [floor(vectorMagnitude(1)), ceil(vectorMagnitude(2))])
% Label the colorbars
ylabel(cb,'Vector magnitude')



% Set properties for the main axes
axis equal
xlim(ax2, [-1 1])
ylim(ax2, [-1 1])
ax2.Color = 'none';
hold off

Answer 2

[Adam Danz]

Let's break this down into steps. This solution doesn't address all of the requirements in the question.

Get the arrows to start at a contour

Step 1 is getting the contour line coordinates. This isn't easy so there are several functions on the file exchange that solve this problem. I'll use getContourLineCoordinates which you'll need to download and add to the MATLAB path.

The contour lines coordinates are not sampled at a uniform interval but you probably don't want an arrow at each coordinate anyway. You could resample each curve to get coordinates along the curves at a uniform interval but I used an easier approach. The variable qint=3 means use every 3 points along the curve.

To get the horizontal and vertical gradient components at each of those coordinates, use interp2 to interpolate the gradient matrices. Then use the result to plot the quiver arrows.

I also recommending turning auto scaling off (5th input in quiver, below).

f = @(x,y) 0.5*sin(pi/2*(2*x-y)) + 2*y + 1;
[X,Y] = meshgrid(-1:0.05:1,-1:0.05:1);

Z = f(X*cos(pi/6)-Y*sin(pi/6),Y*cos(pi/6)+X*sin(pi/6));
[U,V] = gradient(Z);

ax = axes(); 
hold(ax,'on')

[~, h] = contour(ax, X,Y,Z,10);

% https://www.mathworks.com/matlabcentral/fileexchange/74010-getcontourlinecoordinates
CT = getContourLineCoordinates(h); 

qint = 3;  % quiver interval along contour lines (non-zero integer)
for i = unique(CT.Group)'
    cx = CT.X(CT.Group == i); 
    cy = CT.Y(CT.Group == i); 
    idx = mod(1:numel(cx),qint)==0; 
    uInterp = interp2(X,Y,U,cx(idx),cy(idx));
    vInterp = interp2(X,Y,V,cx(idx),cy(idx));
    quiver(cx(idx),cy(idx),uInterp,vInterp,'off','k')
end

cmap = spring(255); 
ax.Colormap = cmap; 

Get the arrows to start at a contour and end at the next contour

🤔 this is a head scratcher. Closer spacing of contours indicate that elevation is changing faster over a smaller distance. In other words, there's a steeper gradient. A steeper gradient would be depicted as longer quiver arrows. So, there's an inverse releationship between the distance between contours and the length of quiver arrows. I'll skip this step.

So there should be no arrows "crossing" a contour line

If this is what you're concerned with, you could try scaling the quiver magnitudes using the 5th input argument that I currently have set to 'off', although it may be difficult to algorithmically choose a scale factor.

Color the quiver arrows based on magnitude

This we can do. But it will require plotting each arrow individually. In this case, it's critical to turn off auto-scaling so we can be sure maintain the relative differences between each arrow. Plotting single quiver arrows results in different arrow head sizes so this solution also increases the maxHeadSize to compensate for smaller arrow heads.

f = @(x,y) 0.5*sin(pi/2*(2*x-y)) + 2*y + 1
[X,Y] = meshgrid(-1:0.05:1,-1:0.05:1);

Z = f(X*cos(pi/6)-Y*sin(pi/6),Y*cos(pi/6)+X*sin(pi/6));
[U,V] = gradient(Z);

ax = axes(); 
hold(ax,'on')

[~, h] = contour(ax, X,Y,Z,10);

% https://www.mathworks.com/matlabcentral/fileexchange/74010-getcontourlinecoordinates
CT = getContourLineCoordinates(h); 

UVMagnitudes = hypot(U,V); 
minMag = min(UVMagnitudes,[],'all'); 
maxMag = max(UVMagnitudes,[],'all');
nColors = 256;  
quivColors = cool(nColors);  
magnitudes = linspace(minMag, maxMag, nColors)'; 

qint = 3;  % quiver interval along contour lines (non-zero integer)
for i = unique(CT.Group)'
    cxg = CT.X(CT.Group == i); 
    cyg = CT.Y(CT.Group == i); 
    idx = mod(1:numel(cxg),qint)==0; 
    cx = cxg(idx);
    cy = cyg(idx); 
    uInterp = interp2(X,Y,U,cx,cy);
    vInterp = interp2(X,Y,V,cx,cy);
    for j = 1:sum(idx)
        colr = interp1(magnitudes, quivColors, hypot(uInterp(j),vInterp(j)));
        quiver(cx(j),cy(j),uInterp(j),vInterp(j),'off','k','color',colr,'MaxHeadSize',.8)
    end
end

cmap = spring(255); 
ax.Colormap = cmap; 

You can even add two colorbars by following this doc example.