% Phoenix ASTER Script

% This script will load an ASTER scene of the Phoenix area using remote 
% sensing data to formulate four figures:  CIR, SWIR, TIR and NDVI.

% Created by Christy Baker on:  10/24/11
% Modified on:  10/25/11
% ----------------------------------------------------------------------

%--------------------------VNIR or CIR----------------------------------
% Pull bands 3N, 2, and 1 from the HDF file into their own matrices.  These
% bands are in the VNIR (aka CIR)range and will highlight actively
% photosynthesizing vegetation (shown as red), undisturbed bedrock and 
% soils (browns, greens and grays), and the built environment (blues,
% greens, reddish-purples and whites)
% ----------------------------------------------------------------------
Band_3 = hdfread('PhoenixASTERData.hdf','VNIR_Swath','Fields','ImageData3N');
Band_2 = hdfread('PhoenixASTERData.hdf','VNIR_Swath','Fields','ImageData2');
Band_1 = hdfread('PhoenixASTERData.hdf','VNIR_Swath','Fields','ImageData1');

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

% Perform histogram equalization stretch on the image
for i=1:3
    Phoenix_321_stretch(:,:,i)=histeq(Phoenix_321(:,:,i));
end
% This imagery will be used to heighten images seen in the visible light
% spectrum.

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

% ---------------------------- SWIR ---------------------------------------
% Pull bands 8, 6, and 4 from the HDF file into their own matrices. These
% bands are in the SWIR range and will highlight spectral features
% diagnostic for iron oxides, illite, kaolinite, and carbonates.
Band_8 = hdfread('PhoenixASTERData.hdf','SWIR_Swath','Fields','ImageData8');
Band_6 = hdfread('PhoenixASTERData.hdf','SWIR_Swath','Fields','ImageData6');
Band_4 = hdfread('PhoenixASTERData.hdf','SWIR_Swath','Fields','ImageData4');

% Construct 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

% This type of remote sensing is best illustrated by use of determining
% between sedimentary rock units.

% Display the image
figure(2)
imshow(Phoenix_864_stretch)
xlabel('Longitude (pixels)')
ylabel('Latitude (pixels)')
title('ASTER Bands 8, 6, 4')

% -------------------------- TIR ------------------------------------------
% Pull bands 13, 12, and 10 from the HDF file into their own matrices.
% These bands are in the TIR range and will highlight spectral features
% diagnostic for silicates, iron- and magnesium-bearing minerals and rocks,
% and carbonates. Quartzites are bright red, basaltic rocks are blue,
% granitoids are purple-violet, and carbonates are greenish yellow.
Band_13 = hdfread('PhoenixASTERData.hdf','TIR_Swath','Fields','ImageData13');
Band_12 = hdfread('PhoenixASTERData.hdf','TIR_Swath','Fields','ImageData12');
Band_10 = hdfread('PhoenixASTERData.hdf','TIR_Swath','Fields','ImageData10');

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

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

% This type of imagery is best utilized in determining difference in
% composition of varying bedrock units.

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

% ---------------------------- NDVI ---------------------------------------
% Compute the Normalized Difference Vegetation Index (NDVI) using the following
% formula: [(Band_3 - Band_2)./(Band_3 + Band_2)]. The NDVI shows the relative
% abundance of actively photosynthesizing vegetation. In other words, it
% assesses how healthy the vegetation is!

% Pull the bands.
Band_3 = hdfread('PhoenixASTERData.hdf','VNIR_Swath','Fields','ImageData3N');
Band_2 = hdfread('PhoenixASTERData.hdf','VNIR_Swath','Fields','ImageData2');
Band_1 = hdfread('PhoenixASTERData.hdf','VNIR_Swath','Fields','ImageData1');

% Calculate NDVI.
% Before we can calculate this, let's look at the format (aka "class" of
% Bands 1, 2, and 3 in the MATLAB Workspace. You will notice that they are
% "uint8", which means they are 8-bit integers! This isn't useful to us,
% especially if we want to do division because division of integers can
% only result in whole numbers. So, we first need to convert these bands to
% a class known as "single" so that we can do some band ratioing!
Band_2 = im2single(Band_2);
Band_3 = im2single(Band_3);
Phoenix_NDVI = (Band_3 - Band_2) ./ (Band_3 + Band_2);

% Display the image.
figure(4)
% Note that because the range of numbers in our NDVI ratio is between
% 0 and 1, we need to display the image's colors using an appropriate scale
% [0 1]
imshow(Phoenix_NDVI,'DisplayRange',[0 1])
xlabel('Longitude (pixels)')
ylabel('Latitude (pixels)')
title('NDVI computed from ASTER Image')