How to combine geoaxes and contourf plot?

117 views (last 30 days)
I would like to combine geoaxes with satellite view and a contourf plot. For example:
geoaxes(); geolimits([20, 30], [-90, -80]); geobasemap('satellite')
And now, given longitude, latitude and sea surface temperature, plot the sea surface temperature as a contour plot. For instance:
contourf(lon, lat, SST)
But this throws this error:
Error using newplot (line 81)
Adding Cartesian plot to geoaxes is not supported.
Error in contourf (line 75)
cax = newplot(cax);
Trying the corresponding function from the Mapping Toolbox:
contourfm(lon, lat, SST)
Error using hggroup
Group cannot be a child of GeographicAxes.
Error in internal.mapgraph.HGGroupAdapter (line 62)
g = hggroup('Parent',ax);
Error in internal.mapgraph.ContourGroup (line 282)
h = h@internal.mapgraph.HGGroupAdapter(args{:});
Error in internal.mapgraph.GeographicContourGroup (line 55)
h = h@internal.mapgraph.ContourGroup(varargin{:});
Error in contourm (line 111)
h = internal.mapgraph.GeographicContourGroup(ax, Z, R, levelList);
Error in contourfm (line 39)
contourm(varargin{:},'Fill','on','DefaultLineColor','black');
I am wondering if it is possible to combine the satellite land view offered by geoaxes together with a contour plot of a 2D geographical array. Or perhaps is it possible to get a satellite land view with the Mapping Toolbox and then use contourfm?
Thanks in advance

Accepted Answer

Adam Danz
Adam Danz on 26 Mar 2021
Edited: Adam Danz on 27 May 2021
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';
  3 Comments
Justino Martinez
Justino Martinez on 27 May 2021
Thanks @Adam Danz, this was my first option but I haven't be capable to use "patchesm" or "patchm" function for geoaxes... Then, I use "patch" in a regular axes to plot line that changes its color depending on a third parameter
:(

Sign in to comment.

More Answers (0)

Products


Release

R2020b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!