clear all
close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Jeffrey Massey
% 2/1/2015
% Edited by Chris Foster
% 11/1/2015
% University of Utah
%
% Purpose: index MODIS albedo to WRF lat/lon grid for later use
%
% download .hdf MODIS albedo product from here (adjust date):
% http://e4ftl01.cr.usgs.gov/MOTA/MCD43B3.005/2012.10.07/
%
% I opt to only consider the SW broadband albedo following:
% http://onlinelibrary.wiley.com/doi/10.1029/2008JD009969/full
% "In this study, the shortwave broadband albedo is used as it is the
% closest spectral match to the broadband pyranometers used by the field measurements."
%
% Unfortunately northern Utah lies between h09v05 and h09v04 so need to
% download both
%
% Therefore the script requires:
% 1. MCD43B3.AYYYYDDD.h09v04.005.YYYYXXXXXXXXX.hdf  
% 2. MCD43B3.AYYYYDDD.h09v05.005.YYYYXXXXXXXXX.hdf  
% 3. wrf met em file that albedo will be applied to
% 4. the MATLAB central script findIndexes.m available from here:
%      http://www.mathworks.com/matlabcentral/fileexchange/37033-find-pixel-indices-in-hdf-eos-files-based-on-latlon-coordinates/content/findIndexes.m
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%where modis files are located
path='/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/matlab/modis/';

%to create the gird (takes awhile...) (1) or just display the albedo (0).
%can only pick (0) if lat and lon has been created though
createlatlon=1;

%dates of modis files (directories in path)
events={'20110704'};

for event=1:length(events)
    
    modislist= dir([path char(events(event)) '/*.hdf']);
    
    for f=1:2
        if f==1
            filename=modislist(1,1).name;
        elseif f==2
            filename=modislist(2,1).name;
        end
        disp(filename)
        file=[path char(events(event)) '/' filename];
        
        sat=hdfsd('start',file,'read');
        for k=0:20
            ref1km =hdfsd('select',sat, k);
            [name,numdim,dimvector,type,numdescr] = hdfsd('getinfo', ref1km);
            %disp([num2str(k) '___' name])
        end
        
        albedo =hdfsd('select',sat, 19);
        [name,numdim,dimvector,type,numdescr] = hdfsd('getinfo', albedo);
        
        %%
        startvector = [0 0];
        endvector = dimvector;
        stride = [] ;
        
        albedovar = hdfsd('readdata', albedo, startvector, stride, endvector);
        
        albedovar=double(albedovar).*.001;
        albedovar(albedovar>1)=nan;
        %%
        
        % get lat/lon
        wrffile='/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/kingspeak/wrf-nlcd2011/WPS/met_em.d03.2011-01-01_06:00:00.nc';
        
        wrflat = double(ncread(wrffile,'XLAT_M'));
        wrflon = double(ncread(wrffile,'XLONG_M'));
        wrfalbedo = double(ncread(wrffile,'ALBEDO12M'));
        
        x_min = 150; 
        x_max = 250;
        y_min = 280;
        y_max = 410;
        
        wrflat = wrflat(x_min:x_max,y_min:y_max);
        wrflon = wrflon(x_min:x_max,y_min:y_max);
        wrfalbedo = wrfalbedo(x_min:x_max,y_min:y_max,:);
        
        wrfalbedo=squeeze(double(wrfalbedo(:,:,10)));
        y1=reshape(wrflat,numel(wrflat),1);
        x1=reshape(wrflon,numel(wrflon),1);
        if createlatlon==1
            LatLon=[y1 x1];
            fprintf('start')
            [X,Y,projType,mdist2point,xRes,yRes] = findIndexes(LatLon,file,1);
            fprintf('end')
            %%
            if f==1
                c=0;
                for i=1:length(X)
                    if mdist2point(i)>1000
                        albedoshow(i,1)=nan;
                        flag(i,1)=1;
                        continue
                    end
                    
                    try albedoshow(i,1)=albedovar(X(i),Y(i));
                        c=c+1;
                        flag(i,1)=0;
                        X_h09v04(c,1:2)=[X(i) y1(i)];
                        Y_h09v04(c,1:2)=[Y(i) x1(i)];
                    catch err
                        albedoshow(i,1)=nan;
                        flag(i,1)=1;
                    end
                end
            elseif f==2
                c=0;
                for i=1:length(X)
                    
                    if flag(i,1)==1
                        if mdist2point(i)>1000
                            albedoshow(i,1)=nan;
                            continue
                        end
                        try albedoshow(i,1)=albedovar(X(i),Y(i));
                            c=c+1;
                            X_h09v05(c,1:2)=[X(i) y1(i)];
                            Y_h09v05(c,1:2)=[Y(i) x1(i)];
                        catch err
                            albedoshow(i,1)=nan;
                        end
                    end
                end
            end
        elseif createlatlon==0
            
            load modis_mcd43B3_latlon_4dwxd03.mat
            latshow=reshape(y1,[size(wrflat,1),size(wrflat,2)]);
            lonshow=reshape(x1,[size(wrflat,1),size(wrflat,2)]);
            if f==1
                albedoshow=nan(size(latshow,1),size(latshow,2));
            end
            for i=1:size(latshow,1)
                for j=1:size(latshow,2)
                    if f==1
                        ind=find(X_h09v04(:,2)==latshow(i,j) & Y_h09v04(:,2)==lonshow(i,j));
                        if ~isempty(ind)
                            albedoshow(i,j)=albedovar(X_h09v04(ind,1),Y_h09v04(ind,1));
                        end
                    elseif f==2
                        ind=find(X_h09v05(:,2)==latshow(i,j) & Y_h09v05(:,2)==lonshow(i,j));
                        if ~isempty(ind)
                            albedoshow(i,j)=albedovar(X_h09v05(ind,1),Y_h09v05(ind,1));
                        end
                    end
                end
            end
        end
        
    end
    %%
    if createlatlon==1
        %save results for future scripts to make things run faster
        %save('modis_mcd43B3_latlon_4dwxd03.mat','X_h09v04','Y_h09v04','X_h09v05','Y_h09v05')
    end
    %%
    close all
    figure
    
    lat1=wrflat(1,:);
    lon1=wrflon(1,:);
    lat2=wrflat(:,end);
    lon2=wrflon(:,end);
    lat3=wrflat(end,:);
    lon3=wrflon(end,:);
    lat4=wrflat(:,1);
    lon4=wrflon(:,1);
    xmin=min(min(wrflon));
    xmax=max(max(wrflon));
    ymin=min(min(wrflat));
    ymax=max(max(wrflat));
    
    load conus

    hgt = double(ncread(wrffile,'HGT_M'));
    
    axesm('lambert','MapLatLimit',[ymin ymax],'MapLonLimit',[xmin xmax],'mapparallels',[30 60])
    gridm off; mlabel off; plabel off
    hold on
    tightmap
    if createlatlon==1
        geoshow(reshape(y1,[size(wrflat,1),size(wrflat,2)]),reshape(x1,[size(wrflat,1),size(wrflat,2)]),reshape(albedoshow,[size(wrflat,1),size(wrflat,2)]),'displaytype','texturemap')
        hold on
        contourm(wrflat,wrflon,hgt(x_min:x_max,y_min:y_max),10,'-w')
    elseif createlatlon==0
        geoshow(latshow,lonshow,albedoshow,'displaytype','texturemap')
    end
    
    caxis([0 0.5])
    colormap('jet')
    geoshow(statelat,statelon,'DisplayType','line','Color','k','linewidth',1);
    
    colorbar
    
    title([char(events(event)) ' MODIS Albedo'])
    
    imagedir=['~/public_html/modis/'];
    name=[imagedir char(events(event)) '.png'];
    saveas(gcf,name,'png')
    
    %display WRF
    figure
    axesm('lambert','MapLatLimit',[ymin ymax],'MapLonLimit',[xmin xmax],'mapparallels',[30 60])
    gridm off; mlabel off; plabel off
    hold on
    tightmap
    
    geoshow(reshape(y1,[size(wrflat,1),size(wrflat,2)]),reshape(x1,[size(wrflat,1),size(wrflat,2)]),wrfalbedo./100,'displaytype','texturemap')
    
    caxis([0 0.5])
    colormap('jet')
    
    geoshow(statelat,statelon,'DisplayType','line','Color','k','linewidth',1);
    %geoshow([S_military(167,1).Lat],[S_military(167,1).Lon],'color','k','LineWidth',1);
    
    colorbar
    
    title([char(events(event)) ' USGS'])
    
    imagedir=['~/public_html/modis/'];
    name=[imagedir 'WRF_BASE_' char(events(event)) '.png'];
    saveas(gcf,name,'png')
    
end



