How to combine geoaxes and contourf plot?

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 = [e-mail](args{:});

Error in internal.mapgraph.GeographicContourGroup (line 55)
            h = [e-mail](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

Answer 1

[Adam Danz]

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';

[Justino Martinez]

Beware! I think that the adjust between both axes is not complete. The limits are well adjusted as well as the longitude grid, but bot the latitude grid... At least in my case... Look this example. I have put as "visible" both axis. The correspondence in longitude and for upper and lower latitudes is correct but not for central latitudes. I think that the differences come from the spherical Earth... The distance in latitude is about 333.585 km considering an spherical earth and 3 degrees in latitude (35-38) but considering a flat earth the distance is 333.547 km The difference is small and I don't know if is the cause of the diffference....

[Adam Danz]

@Justino Martinez, Thanks for pointing that out. I just updated the answer and added option 2. It fixes the problem and is much cleaner and more efficient than option 3 which is the one you pointed out.

@Diego Pereiro Rodríguez, sorry for not thinking of this solution earlier. Consider updating your code.

Note that option 2 requires a function from the file exchange.

[Justino Martinez]

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

:(

[bayan]

I want to fit contour and surface on geomap in matlab. Could I fit I want to fit contour and surface on geomap in matlab?

[John Kalogiros]

For Option 3 to work correctly (at least for small Earth curvature effects, i.e. small regions), just add a listener for ax1 after the listeners for ax2:

ax1.UserData.listener=addlistener(ax1,'MarkedClean',@(~,~)set(ax2,'XLim',ax1.LongitudeLimits,'YLim',ax1.LatitudeLimits));

Also, hide the toolbar of ax2: axtoolbar(ax2,'Visible','off');

Answer 2

[Shai Katz]

Hi @Adam Danz,

I looked at your suggestion and I have a similar qestion with the error:

Adding Cartesian plot to geoaxes is not supported.

I would like to plot wind vectors on geobasemap:

S= wind_speed;
D = 270 - wind_direction;
u = S.*cosd(D);
v = S.*sind(D);
figure(1)
quiver(lat, lon, u, v);
axis equal
grid on 
hold on;
gx = geoaxes;
geobasemap(gx,'topographic');

I want in this way because geobasemape has maps with better visualtion of the surfac such as topography, soil, etc.

Do you have any idea how to make it?

Thanks,

Shai

[Adam Danz]

I've come across this problem in the past and just used map axes. That's probably not the answer you're looking for 🙁 .

https://www.mathworks.com/matlabcentral/answers/628408-topography-in-a-pcolor-graph#answer_526288

Also see quiverm.