Approach for plotting with an indeterminate number of colormaps
2 次查看(过去 30 天)
显示 更早的评论
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
0 个评论
回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Convert Image Type 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!