% Nathan Hanchey
% 10/28/15
% ASTER image processing Phoenix assignment

% This script contains an uploaded ASTER file of the greater Phoenix area
% and four figures (CIR, SWIR, TIR, and NDVI)

% HDF GUI fire up
hdftool('PhoenixASTERData (1).hdf');

% ------------------ VNIR or CIR ----------------------
Band_1 = hdfread('PhoenixASTERData (1).hdf', '/VNIR/VNIR_Swath/Data Fields/ImageData1', 'Index', {[1  1],[1  1],[4200  4980]});
Band_2 = hdfread('PhoenixASTERData (1).hdf', '/VNIR/VNIR_Swath/Data Fields/ImageData2', 'Index', {[1  1],[1  1],[4200  4980]});
Band_3N = hdfread('PhoenixASTERData (1).hdf', '/VNIR/VNIR_Swath/Data Fields/ImageData3N', 'Index', {[1  1],[1  1],[4200  4980]});

% Constructed an RGB matrix from the three bands into one 3D matrix.
Phoenix_321 = cat(3,Band_3N,Band_2,Band_1);

% image
figure(1)
imshow(Phoenix_321)
xlabel('Longitude (pixels)')
ylabel('Latitude (pixels)')
title('ASTER Bands 3, 2, 1')

% --------------- SWIR -----------------
Band_4 = hdfread('PhoenixASTERData (1).hdf', '/SWIR/SWIR_Swath/Data Fields/ImageData4', 'Index', {[1  1],[1  1],[2100  2490]});
Band_6 = hdfread('PhoenixASTERData (1).hdf', '/SWIR/SWIR_Swath/Data Fields/ImageData6', 'Index', {[1  1],[1  1],[2100  2490]});
Band_8 = hdfread('PhoenixASTERData (1).hdf', '/SWIR/SWIR_Swath/Data Fields/ImageData8', 'Index', {[1  1],[1  1],[2100  2490]});

% Constructed an RGB matrix from the three bands into one 3D matrix.
Phoenix_864 = cat(3 , Band_8 , Band_6 , Band_4);

% Perform histogram equalization stretch on the image
for i=1:3
    Phoenix_864_stretch(:,:,i)=histeq(Phoenix_864(:,:,i));
end
% image
figure(2)
imshow(Phoenix_864_stretch)
xlabel('Longitude (pixels)')
ylabel('Latitude (pixels)')
title('ASTER Bands 8, 6, 4')

% -------------- TIR ------------------
Band_10 = hdfread('PhoenixASTERData (1).hdf', '/TIR/TIR_Swath/Data Fields/ImageData10', 'Index', {[1  1],[1  1],[700  830]});
Band_12 = hdfread('PhoenixASTERData (1).hdf', '/TIR/TIR_Swath/Data Fields/ImageData12', 'Index', {[1  1],[1  1],[700  830]});
Band_13 = hdfread('PhoenixASTERData (1).hdf', '/TIR/TIR_Swath/Data Fields/ImageData13', 'Index', {[1  1],[1  1],[700  830]});

% Constructed an RGB matrix from the three bands into one 3D matrix.
Phoenix_131210 = cat(3 , Band_13 , Band_12 , Band_10);

% Decorrelation stretch to enhance differences in colors.
Phoenix_131210_decorr = decorrstretch(Phoenix_131210,'Tol',0.01);

% image
figure(3)
imshow(Phoenix_131210_decorr)
xlabel('Longitude (pixels)')
ylabel('Latitude (pixels)')
title('ASTER Bands 13, 12, 10')

% -------------- NDVI ------------------
Band_1 = hdfread('PhoenixASTERData (1).hdf', '/VNIR/VNIR_Swath/Data Fields/ImageData1', 'Index', {[1  1],[1  1],[4200  4980]});
Band_2 = hdfread('PhoenixASTERData (1).hdf', '/VNIR/VNIR_Swath/Data Fields/ImageData2', 'Index', {[1  1],[1  1],[4200  4980]});
Band_3N = hdfread('PhoenixASTERData (1).hdf', '/VNIR/VNIR_Swath/Data Fields/ImageData3N', 'Index', {[1  1],[1  1],[4200  4980]});

% Calculating NDVI.
Band_2 = im2single(Band_2);
Band_3N = im2single(Band_3N);
Phoenix_NDVI = (Band_3N - Band_2) ./ (Band_3N + Band_2);

% image.
figure(4)
imshow(Phoenix_NDVI,'DisplayRange',[0 1])
xlabel('Longitude (pixels)')
ylabel('Latitude (pixels)')
title('NDVI computed from ASTER Image')