Use of hgtransform with axesm / mapping toolbox

I would like to plot a series of fields (these could be surfacem objects with zero-value ZData or pcolorm objects) in the same figure (could also be the same axes). The fields should be vertically offset from one another. The first figure created by the minimum working example below is what I'm interested in achieving. The code starting at '% with mapping toolbox' is just a starting point; I'm not sure why the single surface I'm plotting is not appearing.

[x,y,z] = peaks;    % data

% define a domain
domain.lon = [20 30];
domain.lat = [30 45];

% scale the data
lon = ((x + min(min(x))) ./ max(max(x))) .* ...
    (max(domain.lon) - min(domain.lon)) + min(domain.lon);
lat = ((y + min(min(y))) ./ max(max(y))) .* ...
    (max(domain.lat) - min(domain.lat)) + min(domain.lat);

% without mapping toolbox
figure;
hold on
view(3)
g1 = hgtransform('Matrix',makehgtform('translate',[0 0 -0.5]));
h1 = surface(lat, lon, zeros(size(z)), z, 'Parent', g1);
g2 = hgtransform('Matrix',makehgtform('translate',[0 0 0.5]));
h1 = surface(lat, lon, zeros(size(z)), z, 'Parent', g2);
zlim([-1 1])

% with mapping toolbox
figure;
view(3)
ax1 = axesm('mapprojection','miller', ...
    'frame','on');
setm(ax1, 'maplatlimit', domain.lat, ...
    'maplonlimit', domain.lon, ...
    'flatlimit', domain.lat, ...
    'flonlimit', domain.lon);
axis normal
tightmap
h1 = surfacem(lat, lon, zeros(size(z)), z);
zlim([-1 1])

Answer 1

[Prashant Arora]

Hi Chris,

That lat and lon variables in your code span from [0 30] and [0 20] respectively, as opposed to the actual limits set on the axes by domain.lat and domain.lon. The scaling equations are probably not correct. Also, in surfacem, the fourth argument should be "ZData" and third corresponds to "CData". You can verify the above statements by using the following code

figure;
view(3)
ax1 = axesm('mapprojection','miller', ...
    'frame','on');
setm(ax1, 'maplatlimit', [0 30], ...
    'maplonlimit', [0 20], ...
    'flatlimit', [0 30], ...
    'flonlimit', [0 20]);
axis normal
tightmap
h1 = surfacem(lat, lon,z,zeros(size(z)));
zlim([-1 1])

Given x and y are symmetric around origin, you can use the following equations to calculate lon and lat.

lat = y.*(max(domain.lat) - min(domain.lat))/(max(y(:)) - min(y(:))) + mean(domain.lat);
lon = x.*(max(domain.lon) - min(domain.lon))/(max(x(:)) - min(x(:))) + mean(domain.lon);