Approach for plotting with an indeterminate number of colormaps

The following code allows multiple landsat tiles to be plotted using the m_map package. The variable IMAGERY_PATH is a cell which is intended to hold an indeterminate number of paths in which the landsat band .tif files reside. Presently, only the first (or last?) tile is rendered with a sensible color. It would seem that a new colormap should be used for each tile.

Any comments on suggested approach are greatly appreciated. I'm happy to provide a link to a set of relatively small data files.

% plot_multiple_landsat_tiles.m

addpath 'C:\matlab\L8read\'
addpath 'C:\matlab\m_map1_4\m_map'

% constants
proj = 'UTM';
datum = 'wsg84';

%%{
LON = [-81.5, -74.5];    
LAT = [-13.5, -8.5];
%}

%{
LON = [-80, -77];    
LAT = [-13.5, -12.5];
%}


% landsat data

% IMAGERY_PATH = {'...\LC08_L1TP_007069_20170408_20170414_01_T1\'; ...
%     '...\LC08_L1TP_009067_20170406_20170414_01_T1\'};

% IMAGERY_PATH = {'...\LC08_L1TP_007069_20170408_20170414_01_T1\'; ...
%     '...\LC08_L1TP_006071_20170401_20170414_01_T1\'};

IMAGERY_PATH = {'...\LC08_L1TP_007069_20170408_20170414_01_T1\'};

num_source_images = length(IMAGERY_PATH);

% instantiate figure before loop

ax1 = axes;                 % instantiate axes
hold(ax1,'on')              % hold on

set(gcf,'Color','w')
set(ax1,'Color','none')
set(ax1,'Box','on')         % box on

m_proj(proj, ...
    'long', LON, ...
    'lat', LAT, ...
    'ellipse', datum);              % map projection

m_grid;                             % show grid
set(ax1, 'Color', 'none')           % set axis background color to none
m_gshhs_f('patch', [0.8 0.8 0.8]);  % coastline patch


for i = 1: num_source_images
    
    imagery_path = IMAGERY_PATH{i};
    
    [I_B2, ~, ~, info] = L8read(imagery_path, LAT, LON, 'band', 2, 'stretch');
    [I_B3, ~, ~, ~] = L8read(imagery_path, LAT, LON, 'band', 3, 'stretch');
    [I_B4, ~, ~, ~] = L8read(imagery_path, LAT, LON, 'band', 4, 'stretch');
    
    % composite
    I = cat(3, I_B4, I_B3, I_B2);
    
    % convert the image to an indexed image
    n = 65536;
    [I_indexed, color_map] = rgb2ind(I, n);
    
    % for support of nan's, convert to type double
    I_double = im2double(I_indexed, 'indexed');
    
    % where the image pixel is 1, let it be nan
    I_double(I_double == 1) = NaN;
    
    % get the corner longitude, lattitude coordinates
    image_corner_lon = [info.CornerCoords.Lon(1,1) info.CornerCoords.Lon(1,2)];
    image_corner_lat = [info.CornerCoords.Lat(1,1) info.CornerCoords.Lat(1,3)];
    
    % image corners
    [image_corner_mmapx, image_corner_mmapy] = m_ll2xy(image_corner_lon, image_corner_lat);
    
    % define colormap
    colormap(ax1, color_map);
    
    % display the composite landsat imagery
    im = image(image_corner_mmapx, image_corner_mmapy, I_double);
    
    % not sure what this does, if anything
    im.AlphaData = double(~isnan(I_double));
    
    clear I_B4 I_B3 I_B2 x y info I I_double
    
end