How to create 3D surface with different layers that each have their own colormap in matlab?

Question

How can I create a 3D map of the glacier (bedrock + glacier surface + debris surface) with each a different colormap in Matlab? I tried several things but it does not work. Only 1 colormap is used and not a separate one for each distinct layer. How can I solve this?

close all;
clear;
clc;

% Load data
bedrock = double(geotiffread('bedrock_oct2023.tif'));
icethick = ncread(file_name, 'icethick');
h_debris = ncread(file_name, 'h_debris');

% Define grid for plotting
[X, Y] = meshgrid(1:size(bedrock, 2), 1:size(bedrock, 1));

% Compute initial elevations
elevation_glacier = bedrock + icethick(:, :, 1);
elevation_debris = bedrock + icethick(:, :, 1) + h_debris(:, :, 1);

% Create figure
figure('Position', [100 100 1200 800]);
ax = axes;
hold(ax, 'on');
% Define colormaps
bedrock_cmap = copper(256);
glacier_cmap = winter(256);
debris_cmap = parula(256);

% Plot initial surfaces with different colormaps manually set using CData
bedrock_alpha = double(icethick(:, :, 1) == 0); % Transparent where glacier exists

% Plot bedrock
bedrock_surf = surf(ax, X, Y, bedrock, 'FaceColor', 'interp', 'FaceAlpha', 1, 'AlphaData', bedrock_alpha);
set(bedrock_surf, 'CData', bedrock, 'EdgeColor', 'none'); % Assign bedrock elevation as CData for coloring

% Plot glacier
glacier_surf = surf(ax, X, Y, elevation_glacier, 'FaceColor', 'interp', 'FaceAlpha', 1);
set(glacier_surf, 'CData', elevation_glacier, 'EdgeColor', 'none'); % Assign glacier elevation as CData for coloring

% Plot debris
debris_surf = surf(ax, X, Y, elevation_debris, 'FaceColor', 'interp', 'FaceAlpha', 1);
set(debris_surf, 'CData', elevation_debris, 'EdgeColor', 'none'); % Assign debris elevation as CData for coloring

% Set different color limits and colormaps
% Bedrock surface color limits
caxis([2000 4500]);
colormap(ax, bedrock_cmap);
% Freeze the colormap for bedrock so that changes to colormap do not affect it
freezeColors(bedrock_surf);

% Glacier surface color limits
caxis([2000 4500]);
colormap(ax, glacier_cmap);
freezeColors(glacier_surf);

% Debris surface color limits
caxis([2000 4500]);
colormap(ax, debris_cmap);
freezeColors(debris_surf);

% Add glacier mask
glacier_mask = double(icethick(:, :, 1) > 0); % Convert mask to double for coloring
glacier_mask_elev = bedrock + icethick(:, :, 1);
mask_surf = surf(ax, X, Y, glacier_mask_elev, 'FaceAlpha', 1, 'EdgeColor', 'w');
set(mask_surf, 'FaceColor', 'white'); % Set mask color to white directly
% Set transparency where elevation equals bedrock (no ice)
alpha_data = (glacier_mask_elev ~= bedrock);
set(mask_surf, 'AlphaData', alpha_data, 'FaceAlpha', 'interp');

% Add colorbars
cb1 = colorbar('Position', [0.92 0.1 0.02 0.2]);
colormap(cb1, bedrock_cmap);
title(cb1, 'Bedrock');
cb2 = colorbar('Position', [0.92 0.4 0.02 0.2]);
colormap(cb2, glacier_cmap);
title(cb2, 'Glacier');
cb3 = colorbar('Position', [0.92 0.7 0.02 0.2]);
colormap(cb3, debris_cmap);
title(cb3, 'Debris');

% Set axis properties
xlabel('X (m)');
ylabel('Y (m)');
zlabel('Elevation (m)');
set(gca, 'XDir', 'reverse', 'FontSize', 14);
title('3D View of Glacier, Debris, Bedrock, and Glacier Mask Surface');
set(gcf, 'color', 'white');

% Set view and lighting
view(45, 35);
shading interp;
camlight('left');
camlight('right');
material dull;  % Adjust material properties
lighting phong;

I have tried the code above but the plot generates only 1 colormap for all 3 layers (bedrock + glacier + debris). I want the plot to have a distinct colormap for each of the layers.