
function Temp_Albedo_Processor
%%  This function takes time series data of calculated albedo from MODIS
%  TERRA data product MCD43B and computes a graph of albedo of each pixel
% over a period of seven years and a seperate graph of the averages. This
% is for custom (long lat coordinates) sites from ORNAL DAAC.
%
%
%% -------------------- RAW DATA COLLECTION --------------------------- %%
%************************************************************************%
%
%   So first go to
%%   http://daac.ornl.gov/cgi-bin/MODIS/GLBVIZ/GENERIC.email.pass.time.4a.1.pl
%
%%  Input the long lat coordinates for the site you want.  On the next
% page, choose [MCD43B} Terra and Aqua BRDF and Calculated Albedo.  For the
% kilometer dimensions of the plot, leave it at the default 3, 3 (making it
% a 7x7 km plot).  Want bigger dimensions?
%    Can't help you with that one.  For the m-file to process the data, the
% files from this website have to be imported into the workspace.  MATLAB
% for some reason cannot import that data with dimensions other than
% default (although it can do a 1x1 pixel plot also). I still haven't found
% the reason behind it, the problem is with the MATLAB importation tool and
% not the m-file.
%
%%    For the next page, leave everything on defaults (unless you want
% different optical depth and solar zenith angle values) enter your email
% address, click 'continue' and then 'create subest' on the next page and
% wait for the data to be processed.
%
%%  Once the email is recieved click on the link to enter the processed
% data page.  Now you have to download a bunch of stuff to your computer.
%** 1. Make a folder with this convention for the name:
%    [site #] [Desciptor] [lat] [long]
%  For example:
%    1 Hoffman Forest 1 34.861647 -77.395536
%    -- if you are using a non default optical depth and solar zenith angle
%  include those also
%
%** 2. Scroll down to below the graphs where you will find the land cover
% images.  Download the local and regional landcover images and the legend
% into the fold you just created.  Then go down to the bottom where you
% have a series of download buttons.  Download "MODIS Land Product ASCII
% Data" (named MCD43B.asc) and "Land Cover in ASCII Grid Format" (named 
% "GRID.MOD12Q1.STD.asc").  Both of these files need to be imported into
% MATLAB, but the land cover in asci grid format need to be modified
% because there is a whole lot of metadata that will confuse MATLAB.  All 
% MATLAB needs from that file is the landcover matrix.
%  So open GRID.MOD12Q1.STD.asc with Wordpad- around 15 lines down it list
%  the dimensions and coordiates for the corner of the matrix, and below it
%  is the matrix.  Copy just the matrix below it, should be 7x7, not that
%  bigger regional land cover matrix below it.  Create a new Worpad file,
%  paste the matrix, and save it as "landcover.txt".
%
%  Congrats!  You have just gathered all of the raw data needed for that
% site.
%
%
%% --------------------  DATA PROCESSING   ------------------------------%%
%************************************************************************%
%
%   Now open matlab, set the workspace to where this m-file is (Note: there
% must be a folder named "Data" in the same directory as this m-file,
% that is where all of the processed data files are exported to).
%   Now go to File>Import data, and navigate to where the raw data folder is
% and import MCD43B.asc and landcover.txt, import using the default
% options (this should create the workspace variables 'data', landcover',
% and 'textdata').
%   Now go open up the image file showing the landcover that you downloaded
% from the website (named zoom.12q1.png).  Open it up and choose your
% landcover type you want to process.
%  Run this m-file, a dialog box should appear requesting meta data info.
% Enter in all of the metadata information (this is used for the names of 
% the output data so you know where it comes from).
%   Press ok and another box comes up asking you to input a number for the
% land cover type, input the number that correspond to the landcover type
% (this is used for the pixel filterer) and the number that corresponds 
% to the spectrum type.  Press ok and you should be done.
%   This m-file produces a graph with the average albedo and
% the annual range of albedo shown with error bars.  A .fig file of the
% graph is exported, as well as a .jpeg, and a .mat file with all of the
% statistics for the data.  They are all exported to the /Data directory
% with a unique name for each site and spectrum.
% 
%**************************************************************************
%
%     NOTE: There must be a folder titled 'Data' in the same directory as
%this m-file.  The data will be exported there.  If you want the data
%exported to the same directory as the workspace delete './Data/' on line
%113.
%
%**************************************************************************
%
% Author: Bruce Corliss
% Dr. Rob Jackson
% Duke University
% Version: .8 (beta)
% 7/26/2007



%% --------------------       USER INPUT PROMPT      ------------------- %%
%************************************************************************%

prompt = {'Enter data variable', 'Enter metadata variable', ...
    'Enter Location Name', 'Enter Optical Depth',...
    'Enter solar Zenith Angle', 'Enter Longitude', 'Enter Latitude', ...
    'Enter Projection (Resolution)', 'Enter Land Cover Variable'};
title = 'Specify data table and metadata';
lines = 1;
def = {'data', 'textdata', '', '0.2',  'local', '', '', 'SIN (1km)', 'landcover'};
answer = inputdlg(prompt,title,lines,def);
data= evalin ('base', answer{1});
textdata= evalin ('base', answer{2});
Location= answer{3};
%type= answer{4}; %%%%%
od= answer{4};
szn = answer{5};
long= answer{6};
lat= answer{7};
proj= answer{8};
landcover= evalin ('base', answer{9});

pstring= ['What type of forest is being sampled? \n\n' ...
    '(0) Water \n' '(1) Evergreen Needleleaf Forest \n'...
    '(2) Evergreen Broadleaf Forest \n' '(3) Deciduous Needleleaf Forest \n'...
    '(4) Deciduous Broadleaf Forest \n' '(5) Mixed Forest \n'...
    '(6) Closed Shrublands \n' '(7) Open Shrublands \n'...
    '(8) Woody Savannas \n' '(9) Savannas\n' '(10) Grasslands\n' '(11) Permenent Wetlands \n'...
    '(12) Croplands \n' '(13) Urban \n' '(14) Cropland/natural Mosiac \n'];
p2string= ['Select Spectrum Type: \n 1. Visible\n 2. Shortwave\n 3. NIR'];
prompt = {sprintf(pstring), sprintf(p2string)};
title= 'Forest Extraction Type';
line=1;
def= {'2', '3'};
answer= inputdlg(prompt,title,lines,def);
cover_type= str2num( answer{1});
spectrum= str2num( answer{2});

% Assign the proper values to spectrum so it extracts the right data
switch spectrum
    case 1
        spectrum = 8;
        spec_type= 'VIS';
    case 2
        spectrum = 5;
        spec_type= 'SW';
    case 3
        spectrum = 2;
        spec_type= 'NIR';
end

% the variable type is used for labeling purposes (included in graph title
% and name of exported files).  describes landcover type.
switch cover_type
    case 0
        type= 'Water';
    case 1
        type= 'Evergreen Needleleaf';
    case 2
        type= 'Evergreen Broadleaf';
    case 3
        type= 'Deciduous Needleleaf';
    case 4
        type= 'Deciduous Broadleaf';
    case 5
        type= 'Mixed Forest';
    case 6
        type= 'Closed Shrublands';
    case 7
        type= 'Open Shrublands';
    case 8
        type= 'Woody Savannas';
    case 9
        type= 'Savannas';
    case 10
        type= 'Grasslands';
    case 11
        type= 'Permanant Wetlands';
    case 12
        type= 'Croplands';
    case 13
        type= 'Urban';
    case 14
        type= 'Cropland/Natural Vegetation Mosaic';
end
       
        

%% ---------------------     DATA Processing    ------------------------ %%
%*************************************************************************%

% Only pixels of a certain homogenous forest type are filtered
% This makes an index of all pixels matching the landcover chosen
pixels=find(landcover'==cover_type)';

k=0;
for i=pixels
    %Extracts every 10th value from table to get it all actual_albedo values
    %for ONE square of the grid
    k=k+1;
    test(:,k)=data(spectrum:10:end,i); %replace 'spectrum' with the albedo variable that you want
end                                %if manual selection is needed.  see 'textdata' variable


%Changes all null data values to a flag that matlab recognizes
%Gets array of index coordinates
for i= 1:k
    aa=find(test(:,i)==32767);
    %Matlab's null data value
    test(aa,i)=NaN;
end

%Extract the time data (constant dt) 
timetest=textdata(2:10:end,3);

%Change code to yearday
for i=1:length(timetest)
    yearday(i)=str2num(timetest{i}(2:8));
    day(i)=str2num(timetest{i}(6:8));
    year(i)=str2num(timetest{i}(5));
end

test(test > 0.20)= NaN;  %turn values that are outliers into NAN

% Temporal stats (not to be confused with yearly stats: albedo_mean, albedo_median)
mean_albedo= nanmean(test');   %DIFFERENT FROM albedo_mean
median_albedo= nanmedian(test');


% Annual data
k=1;
for i=unique(year)
   aa=find(year==i);
   
   % Max/Min/range 95 and 5 percentile
   a_max95(k,:)= prctile(test(aa,:),95);
   albedo_max95(k)= nanmax(a_max95(k,:));
   
   a_min5(k,:)= prctile(test(aa,:), 5);
   albedo_min5(k)= nanmin(a_min5(k,:));
   
   albedo_range95(k)= albedo_max95(k)- albedo_min5(k);
   
   % Max/Min/range
   a_max(k,:)= nanmax(test(aa,:));
   albedo_max(k)= nanmax(a_max(k,:));
   
   a_min(k,:)= nanmin(test(aa, :));
   albedo_min(k)= nanmin(a_min(k,:));
   
   albedo_range(k)= albedo_max(k)-albedo_min(k);
   
   % Mean/Median
   a_mean(k,:)= nanmean(test(aa,:));
   albedo_mean(k)= nanmean(a_mean(k,:));
   
   a_median(k,:)= nanmedian(test(aa,:));
   albedo_median(k)= nanmedian(a_median(k,:));
   
   mean_range(k) = range(mean_albedo(aa));
   
   median_range(k) =range(median_albedo(aa));

   k=k+1;
end


%% -----------------------     DATA PLOTS    -------------------------- %%
%*************************************************************************%

% First setup some stats
year= 2000+year;
day= day/365;
yearday= year+day;
stat_time= (unique(year)+.5);

albedo_range = albedo_max-albedo_min;


plot(yearday,mean_albedo,'.k')
hold on
%plot(yearday, test, 'x')  % Plot command for all data points
plot(yearday,mean_albedo,'r')
plot(yearday, median_albedo)
%plot(stat_time, albedo_max)
%plot(stat_time, albedo_min)
h = errorbar(stat_time, albedo_median, (albedo_median-albedo_min5),...
    (albedo_max95-albedo_median));
set(h, 'linestyle', 'none', 'marker', 's', 'color', 'b')
hold off

xlabel('Time')
ylabel('Albedo')
Title(['Blue-Sky Albedo (' spec_type ') ' 'of ' Location ' (' type ') ' ' '...
    num2str(length(pixels)) ' Pixels'])
axis([2003 2007 0 .3])
grid on



%% -----------------------     DATA EXPORTATION    --------------------- %%
%*************************************************************************%

filename= ['./Data/' Location ' ' type ' ' spec_type ' '  od ' ' szn ' ' long...
    ' ' lat ' ' proj];

%Reassign values to columns for export
stat_export(:,1)= stat_time(2:5)';
stat_export(:,2)= albedo_mean(2:5)';
stat_export(:,3)= albedo_median(2:5)';
stat_export(:,4)= albedo_range(2:5)';  %not made yet
stat_export(:,5)= albedo_max95(2:5)';
stat_export(:,6)= albedo_min5(2:5)';
stat_export(:,7)= albedo_range95(2:5); %not made yet
stat_export(:,8)= mean_range(2:5)';
stat_export(:,9)= median_range(2:5)';

stat_hdr= ['Time Mean Median Range Max95 Min5 Range95 RangeOfMean RangeOfMedian'];
pixel_tot= length(pixels);



save([filename '.mat'], 'stat_hdr', 'stat_export', 'pixel_tot', '-mat')
%save([filename 'export' '.mat'], 'final' , '-mat')
saveas(gcf, [filename '.fig'], 'fig')
saveas(gcf, [filename '.jpg'], 'jpg')
display('Data Processed')



%% Extraneous INFO _______________________________________________________
% Error bar extraneous
% %subplot(2,2,2)
% plot(yearday, test, 'x')
% hold on

% hold off
% grid on

% xlabel('Time')
% ylabel('Albedo')
% Title(['Blue-Sky Albedo of ' Location ' (' type ') '])