Option 1: use map axes
If you can use map axes instead of geoaxes then you can use contourfm. This was the solution for a similar question by a user who wanted to combine a quiver plot with geoaxes.
Option 2: get contour coordinates and plot directly on geoaxes
This option added 26-May-2021.
Get the contour line coordinates using contourc which produces a matrix of contour line coordinates but the matrix arrangement is difficult to work with. I'm using my getContourLineCoordinates function from the file exchange to reorganize the matrix into a more useful table (there are several other similar functions on the file exchange). Now it's easy to plot the contour ridge lines directly to the geoaxes.
% Create geoaxes
figure()
ax = geoaxes();
geobasemap('satellite')
geolimits([20, 30], [-90, -80]);
% Produce contour matrix
z = peaks(15);
x = linspace(ax.LongitudeLimits(1),ax.LongitudeLimits(2), 15);
y = linspace(ax.LatitudeLimits(1), ax.LatitudeLimits(2), 15);
contdata = contourc(x,y,z);
cTbl = getContourLineCoordinates(contdata); % from the file exchange
% Plot contour lines
hold(ax,'on')
nContours = max(cTbl.Group);
colors = autumn(nContours);
for i = 1:nContours
gidx = cTbl.Group == i;
geoplot(ax, cTbl.Y(gidx), cTbl.X(gidx), ... % note: x & y switched
'LineWidth', 2, 'Color', colors(i,:))
end
% Add colorbar
colorbar(ax)
Option 3: axis overlay ❌
This is buggy and inefficient. Use option 1 or 2.
Alternatively you could create a second pair of transparent axes that hosts the contour plot directly on top of the geoaxes. An important step in this process is to link the axes positions and the axes limits between the two pairs of axes so that they move together when adding a colorbar, legend, or when panning, zooming, etc. This problem isn't too difficult to solve when the axes types are the same (example) but in your case, one pair of axes are geoaxes and the other will be regular axes and they have different properties (e.g. xlim/ylim vs LongitudeLimits/LatitudeLimits) and axis scales and this makes it difficult to use linkaxes & linkprop.
Here's a quick demo that overlays an invisible axes on top of the geoaxes and hosts the contour plot. It uses a listener to link the axis positions and another listener to update the LongitudeLimits/LatitudeLimits when the xlim/ylim changes.
If you're planning on panning/zooming, axis interaction is little laggy and this hasn't been stress testested too deeply but should get you started.
Problems
- As mentioned in the geolimits documentation, "the actual limits chosen by geolimits are greater in extent than the requested limits"
- The latitude scale is non-linear in geoaxes to account for earth curvature but the y-axis in the overlay is linear so that's a mess.
% Create geoaxes
figure()
ax1 = geoaxes();
geobasemap('satellite')
geolimits([20, 30], [-90, -80]);
% Create overlayed, regular axes with same position, xlim, and ylim.
% We'll turn off the axes at the end.
ax2 = axes('Units', ax1.Units, 'Position', ax1.Position, ...
'xlim', ax1.LongitudeLimits, 'ylim', ax1.LatitudeLimits, ...
'Color','none');
% Add listener that updates geoaxes position when regular axes position changes.
% ax2.UserData.linkprop = linkprop([ax1,ax2],{'Position','InnerPosition'}); % This didn't work as well as below
ax2.UserData.listener1 = addlistener(ax2, 'SizeChanged',@(~,~)set(ax1, 'Position', ax2.Position, 'InnerPosition', ax2.InnerPosition));
% Add listener that updates LongitudeLimits/LatitudeLimits when xlim/ylim changes
% ax2.UserData.listener = addlistener(ax2, {'XLim','YLim'},'PostSet',@(~,~)geolimits(ax1,ylim(ax2),xlim(ax2))); % fails when 'Restor View' is pressed in axis toolbar
ax2.UserData.listener2 = addlistener(ax2, 'MarkedClean',@(~,~)geolimits(ax1,ylim(ax2),xlim(ax2)));
% Plot contour
z = peaks(15);
[x,y] = meshgrid(linspace(ax2.XLim(1), ax2.XLim(2), 15), linspace(ax2.YLim(1), ax2.YLim(2), 15));
h = contour(ax2,x,y,z,'LineWidth', 2);
ax2.Colormap = autumn(255);
% Add colorbar
colorbar(ax2)
% After adding stuff to ax2, turn its visibility off
ax2.Visible = 'off';
