%addpath('/data/srb1/wchen/tools/m_fx/','-end')
addpath('/data/srb1/wchen/tools/m_fx//','-end')
%comparisons between MODIS5km and SURFRAD ground-truth for instantaneous
%cases by wchen 10.13.2020
pth       = './';
StOrg     = 'st_8sts';
nStd      = 1.5;
modSat    = 'terra';%'aqua';%
yyyy      = 2015;
mm        = 7;
%gridR     = 5;
FigName1   = ['ssda_modis_',modSat,'_C07_vs_st_',...
				num2str(yyyy),num2str(mm,'%02d')];         

if nStd == 0
    optOut = 'no';
else
    optOut = 'yes';   
end
xlim	 = [0 1500];ylim     = xlim;
nxbin	 = 100;nybin    = nxbin;

StLon       = [-116.018,-88.37,-105.10,-89.97,-105.237,-77.93,-96.62,-97.49];%
StLat       = [  36.626, 40.06,  48.31,34.25,   40.125, 40.72, 43.73,36.61];%
nSt         = length(StLat);
stid  		= {'dra', 'bon','fpk', 'gwn','tbl','psu','sxf', 'SGPC1'};%
subti  		= {'(a)', '(b)','(c)', '(d)','(e)','(f)','(g)', '(h)'};%

nmin        = 1440;
nday        = 365;

timeRs=[0,5,10,15,30];
t0=datenum(2015,1,1);tn=datenum(2015,12,31);t=t0:tn;dates=datevec(t);
it0=find(dates(:,1)==2015&dates(:,2)==7&dates(:,3)==1);
it1=find(dates(:,1)==2015&dates(:,2)==7&dates(:,3)==31);
count=0;
%++++++++++++++
%+Ground-truth+
%++++++++++++++
 load('ssda_st_8sts.mat','st_flx');st_flx= squeeze(st_flx(:,:,3,:));%[1440 365 3 8]
 for gridR = [5 10 15 25  50]
    load(['modis.',modSat,'.4.st.8sts.R',num2str(gridR,'%02d'),'.201507.mat'],'stFlx') ;mod_est = reshape(stFlx,[288 365 8]);
    XX=mod_est;XX(:)=nan;
    for itimeR=1:5
        timeR=timeRs(itimeR);
         for ist=1:nSt
            for idy=it0:it1
                tmp=mod_est(:,idy,ist);
                indt=find(~isnan(tmp));
                if ~isempty(indt)               
                    imin=(indt-1)*5+1;
                    %disp(imin)
                    if timeR~=0
                            XX(indt,idy,ist)=nanmean(st_flx(imin-timeR:imin+timeR,idy,ist));
                    else
                            XX(indt,idy,ist)=st_flx(imin,idy,ist);
                    end
                        clear indt
                    end %isempty
                    clear tmp
            end %idy
        end %ist
 
%%
        for ist=1:9
            %disp(upper(stid{ist}))
            %(1440/5,nday,nSt)    
            if ist<9
                XXx=XX(:,:,ist);
                YY = mod_est(:,:,ist);
                %figcap = [subti{ist},' SURFRAD/',upper(stid{ist}),'.vs.','MODIS\_EST'];
                figcap = [subti{ist},' ',upper(stid{ist})];
            else
                XXx=XX(:,:,:);
                YY = mod_est(:,:,:);
                figcap = 'On-site.vs.MODIS\_EST';
            end
x=XXx(:);
y=YY(:);
%++++++++++++++++++
 x(isnan(y))=NaN;
%+++++++++++
% Statistics
%+++++++++++
% disp({['MODIS ',modSat],...
%     ['ST = ',stid{ist}],...
%     ['TimeR = ', num2str(timeR, '%02d')],...
%     ['girdR = ', num2str(gridR, '%02d')]})
 [Bias,avg,RMSE,Std,stdbfr,Corr,npts,N,outline]=CalculateStatistics_v5(x,y,1,nStd,optOut);  
 statistics =  {sprintf('Corr: %.2f', Corr),...
                [sprintf('Bias: %.2f', Bias), ' (',num2str(abs(Bias/avg)*100,'%.2f'),'%)'],...
                [sprintf('Std: %.2f', Std),' (',num2str(abs(Std/avg)*100,'%.2f'),'%)'],...
                 [sprintf('RMSE: %.2f', RMSE),' (',num2str(abs(RMSE/avg)*100,'%.2f'),'%)'],...
                 sprintf('N: %d', N)};
 if count==0
    sitename= ist;
    tR=timeR;gR=gridR;
    R= Corr;Mean_diff=Bias;Stds=Std;Stdp=abs(Bias/avg)*100;RMSEs=RMSE;RMSEp=abs(RMSE/avg)*100;Ncase=N;
 else
    sitename=[sitename;ist];
    tR=[tR;timeR];gR=[gR;gridR];     
    R=[R; Corr];Mean_diff=[Mean_diff;Bias];Stds=[Stds;Std];Stdp=[Stdp;abs(Bias/avg)*100];
    RMSEs=[RMSEs;RMSE];RMSEp=[RMSEp;abs(RMSE/avg)*100];Ncase=[Ncase;N]; 
 end
 inds = find(~isnan(x));
 xx   = x(inds);
 yy   = y(inds);
 %disp([xx yy (yy-xx)])
 if timeR==30 && gridR==5 
 [i,j]=ind2sub([1440/5 nday],inds);
 %disp([i j])
 hr=floor(i*5/60);
 min=mod(i*5,60);
 dates=datevec(datenum(2015,1,0)+j);
 dates(:,4)=hr;
 dates(:,5)=min;
 dates=datevec(datenum(dates)-datenum(0,0,0,0,5,0));
 timstr = datestr(dates,'mm/dd/yyyy HH:MM')
 tbl1=table(timstr);
 writetable(tbl1,[modSat,'_datetime_',num2str(ist),'.txt'])
 %disp(tbl1)
% disp(dates)
%  keyboard 

%  disp('-----min/max x/y')
%  disp([min(xx),max(xx);min(yy),max(yy)])
 %disp([length(xx),length(yy)])

%++++++++++
% Plotting
%++++++++++ 
tistrs  = {'Observed SW{\downarrow} (W m^{-2})','Estimated SW{\downarrow} (W m^{-2})',''};

xlim	 = [0 1200];ylim     = xlim;
nxbin	 = 300;nybin    = nxbin;
if ist<9
figure('Units','inches','Position',[0.08 3 18.5 10])
pltdata=nan(length(xx),2);
pltdata(:,1)=xx;
pltdata(:,2)=yy;
bar(pltdata,'linewidth',2.5)
box off
set(gca,'linewidth',2.125,'FontSize', 22);
if ist==1 
legend(tistrs,'Orientation','horizontal','FontSize',18);
end
cap1=annotation('textbox', [0.15 0.80 0.1 0.1],...
 'String',figcap,...
  'EdgeColor', 'none', 'FontSize', 18);
% keyboard
print(gcf,'-djpeg','-r300',['Bar_',FigName1,'_gridR',num2str(gridR),'_timeR',num2str(timeR),'_st',num2str(ist)])


else
 
heatscatters(xx,yy,'100',tistrs,statistics,'100')
 

line(xlim,ylim,'LineStyle','--','color', 'k','linewidth',2.5)

set(gca,'xlim',xlim,'xtick',xlim(1):nxbin:xlim(2),'ylim',ylim,'ytick',ylim(1):nybin:ylim(2),...
'linewidth',2.125,'FontSize', 16);
  
cap1=annotation('textbox', [0.55 0.16 0.1 0.1],...
 'String',figcap,...
  'EdgeColor', 'none', 'FontSize', 16);
text(50,1000,statistics,'FontSize', 16);
% cap2=annotation('textbox', [0.2 0.77 0.1 0.1],...
%  'String','UMD .vs. SURFRAD',...
%   'EdgeColor', 'none', 'FontSize', 16);
ax1=gca;
print(gcf,'-djpeg','-r300',['scatter_',FigName1,'_timeR',num2str(timeR),'_st',num2str(ist)])
end
clear i j hr min dates
end
count=count+1;
clear x y XXx YY inds 
 end %ist
end % timeR
clear XX  mod_est
 end %gridR
tbl2=table(sitename,tR,gR,R,Mean_diff,Stds,Stdp,RMSEs,RMSEp,Ncase)
writetable(tbl2,[modSat,'_stats_std',num2str(nStd),'.txt'])
 %{
%  disp('-----min/max x/y')
%  disp([min(xx),max(xx);min(yy),max(yy)])
 %disp([length(xx),length(yy)])

%++++++++++
% Plotting
%++++++++++ 

tistrs  = {'Observed SW{\downarrow} (W/m^{2})','Estimated SW{\downarrow} (W/m^{2})',''};

 
heatscatters(xx,yy,'100',tistrs,statistics,'100')
 

line(xlim,ylim,'LineStyle','--','color', 'k','linewidth',2.5)

set(gca,'xlim',xlim,'xtick',xlim(1):nxbin:xlim(2),'ylim',ylim,'ytick',ylim(1):nybin:ylim(2),...
'linewidth',2.125,'FontSize', 18,'fontweight','bold','FontName', 'Times');
  
cap1=annotation('textbox', [0.5 0.16 0.1 0.1],...
 'String','UMD .vs. SURFRAD',...
  'EdgeColor', 'none', 'FontSize', 18,'fontweight','bold','FontName', 'Times');
cap2=annotation('textbox', [0.175 0.77 0.1 0.1],...
 'String','(a)',...
  'EdgeColor', 'none', 'FontSize', 18,'fontweight','bold','FontName', 'Times');
ax1=gca;
set(gcf,'PaperUnits', 'inches');
set(gcf,'PaperSize', [2 1]);
% print(gcf,'-djpeg','-r300',[FigName,'_1'])

%----- compare with LIS3 ----------------------
ytmp = estdata(2,:,:);
 y    = ytmp(:);
 disp([min(y),max(y)])
 %++++++++++++++++++
 x(isnan(y))=NaN;
 y(isnan(x))=NaN;
%+++++++++++
% Statistics
%+++++++++++
 [Bias,avg,RMSE,Std,stdbfr,Corr,npts,N,outline]=CalculateStatistics_v5(x,y,1,nStd,optOut);  
%  statistics =  {sprintf('Corr: %.2f', Corr),sprintf('Bias: %.2f', Bias), ...
%  sprintf('Std: %.2f', Std),sprintf('RMSE: %.2f', RMSE),...
%  sprintf('N: %d', N)};  
 statistics =  {sprintf('Corr: %.2f', Corr),...
                [sprintf('Bias: %.2f', Bias), '(',num2str(abs(Bias/avg)*100,'%3.1f'),'%)'],...
                [sprintf('Std: %.2f', Std),'(',num2str(abs(Std/avg)*100,'%3.1f'),'%)'],...
                 [sprintf('RMSE: %.2f', RMSE),'(',num2str(abs(RMSE/avg)*100,'%3.1f'),'%)'],...
                 sprintf('N: %d', N)}
 if strcmpi(optOut,'yes') 
                x(outline)=NaN;
                y(outline)=NaN;
                clear outline
 end 
 inds = find(~isnan(x));
 xx   = x(inds);
 yy   = y(inds);
 disp('-----min/max x/y')
 disp([min(xx),max(xx);min(yy),max(yy)])
heatscatters(xx,yy,'100',tistrs,statistics,'100')
 

line(xlim,ylim,'LineStyle','--','color', 'k','linewidth',2.5)

set(gca,'xlim',xlim,'xtick',xlim(1):nxbin:xlim(2),'ylim',ylim,'ytick',ylim(1):nybin:ylim(2),...
'linewidth',2.125,'FontSize', 18,'fontweight','bold','FontName', 'Times');
  
cap1=annotation('textbox', [0.5 0.16 0.1 0.1],...
 'String','LIS3 .vs. SURFRAD',...
  'EdgeColor', 'none', 'FontSize', 18,'fontweight','bold','FontName', 'Times');
cap2=annotation('textbox', [0.175 0.77 0.1 0.1],...
 'String','(b)',...
  'EdgeColor', 'none', 'FontSize', 18,'fontweight','bold','FontName', 'Times');
%{
ax2=gca;

fnew = figure;
ax1_copy = copyobj(ax1,fnew);
subplot(1,2,1,ax1_copy)

ax2_copy =  copyobj(ax2,fnew);
subplot(1,2,2,ax2_copy)
set(gcf,'PaperUnits', 'inches');
set(gcf,'PaperSize', [4 2]);
%}

print(gcf,'-djpeg','-r300',[FigName,'_2'])
% print(fig,'-depsc2','-r300',['./figs/heatscatters_',FigName])
% print(fig,'-dpng','-r300',['heatscatter_',FigName])
                

%}
function outfile = heatscatters(X, Y, numbins, tistrs,stats,markersize, marker, plot_colorbar, plot_lsf )

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% heatscatter(X, Y, outpath, outname, numbins, markersize, marker, plot_colorbar, plot_lsf, xlab, ylab, title)
% mandatory:
%            X                  [x,1] array containing variable X
%            Y                  [y,1] array containing variable Y
%            outpath            path where the output-file should be saved.
%                                leave blank for current working directory
%            outname            name of the output-file. if outname contains
%                                filetype (e.g. png), this type will be used.
%                                Otherwise, a pdf-file will be generated
% optional:
%            numbins            [double], default 50
%                                number if bins used for the
%                                heat3-calculation, thus the coloring
%            markersize         [double], default 10
%                                size of the marker used in the scatter-plot
%            marker             [char], default 'o'
%                                type of the marker used in the scatter-plot
%            plot_colorbar      [double], boolean 0/1, default 1
%                                set whether the colorbar should be plotted
%                                or not
%            plot_lsf           [double], boolean 0/1, default 1
%                                set whether the least-square-fit line
%                                should be plotted or not (together with
%                                the correlation/p-value of the data
%            xlab               [char], default ''
%                                lable for the x-axis
%            ylab               [char], default ''
%                                lable for the y-axis
%            title              [char], default ''
%                                title of the figure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


    %%%% mandatory
    if ~exist('X','var') || isempty(X)
        error('Param X is mandatory! --> EXIT!');
    end
    if ~exist('Y','var') || isempty(Y)
        error('Param Y is mandatory! --> EXIT!');
    end
%     if ~exist('outpath','var')
%         error('Param outpath is mandatory! --> EXIT!');
%     end
%     if ~exist('outname','var') || isempty(outname)
%         error('Param outname is mandatory! --> EXIT!');
%     end
    
    %%%% optional
    if ~exist('numbins','var') || isempty(numbins)
        numbins = 20;
    else
        % force number, not char input
        numbins = str2double(numbins);
    end
    if ~exist('markersize','var') || isempty(markersize)
        markersize = 50;
    else
        % force number, not char input
        markersize = str2double(markersize);
    end
    if ~exist('marker','var') || isempty(marker)
        marker = '.';
    end
    if ~exist('plot_colorbar','var') || isempty(plot_colorbar)
        plot_colorbar = 1;
    end
    if ~exist('plot_lsf','var') || isempty(plot_lsf)
        plot_lsf = 1;
    end
    if ~exist('tistrs','var') || isempty(tistrs)
        tistrs ={'','',''} ;
    end
    if ~exist('stats','var') || isempty(stats)      
		[r,p] = corr(X, Y);
		str = {sprintf('corr: %.3f', r), sprintf('pval: %d', p)};
		else
		str=stats;
	 end
    [values, centers] = hist3([X Y], [numbins numbins]);
%     disp([min(values(:)),max(values(:))])

    centers_X = centers{1,1};
    centers_Y = centers{1,2};
%     disp([min(centers_X),max(centers_X)])
%     disp([min(centers_Y),max(centers_Y)])

    
    binsize_X = abs(centers_X(2) - centers_X(1)) / 2 ;
    binsize_Y = abs(centers_Y(2) - centers_Y(1)) / 2 ;
    bins_X = zeros(numbins, 2);
    bins_Y = zeros(numbins, 2);

    for i = 1:numbins
        bins_X(i, 1) = centers_X(i) - binsize_X;
        bins_X(i, 2) = centers_X(i) + binsize_X;
        bins_Y(i, 1) = centers_Y(i) - binsize_Y;
        bins_Y(i, 2) = centers_Y(i) + binsize_Y;
    end
%     disp([bins_X;bins_Y])
    scatter_COL = zeros(length(X), 1);

    onepercent = round(length(X) / 100.);
    
    fprintf('Generating colormap...\n');
    
    for i = 1:length(X)

        if (mod(i,onepercent) == 0)
            fprintf('.');
        end            

        last_lower_X = NaN;
        last_higher_X = NaN;
        id_X = NaN;

        c_X = X(i);
        last_lower_X = find(c_X >= bins_X(:,1));
        if (~isempty(last_lower_X))
            last_lower_X = last_lower_X(end);
        else
            last_higher_X = find(c_X <= bins_X(:,2));
            if (~isempty(last_higher_X))
                last_higher_X = last_higher_X(1);
            end
        end
        if (~isnan(last_lower_X))
            id_X = last_lower_X;
        else
            if (~isnan(last_higher_X))
                id_X = last_higher_X;
            end
        end

        last_lower_Y = NaN;
        last_higher_Y = NaN;
        id_Y = NaN;

        c_Y = Y(i);
        last_lower_Y = find(c_Y >= bins_Y(:,1));
        if (~isempty(last_lower_Y))
            last_lower_Y = last_lower_Y(end);
        else
            last_higher_Y = find(c_Y <= bins_Y(:,2));
            if (~isempty(last_higher_Y))
                last_higher_Y = last_higher_Y(1);
            end
        end
        if (~isnan(last_lower_Y))
            id_Y = last_lower_Y;
        else
            if (~isnan(last_higher_Y))
                id_Y = last_higher_Y;
            end
        end
%         disp([num2str(id_X),' ',num2str(id_Y),' ',num2str(values(id_X, id_Y))])
        scatter_COL(i) = values(id_X, id_Y);
        
    end
    %disp([min(scatter_COL),max(scatter_COL)])
    %fprintf(' Done!\n');
    
    fprintf('Plotting...');
    
    f = figure();
    %colormap(jet)
     load BlAqGrYeOrReVi200.rgb %BlGrYeOrReVi200.rgb
     colormap(BlAqGrYeOrReVi200/255.) %(BlGrYeOrReVi200/255.)
    
    %scatter(X, Y, markersize, scatter_COL, marker);
    scatter(X, Y, markersize, 'b', marker);
    plot_colorbar=0;
	ax = gca;
	box(ax,'on')
    
    if (plot_colorbar)
        colorbar;      
    end
    
    if (plot_lsf)
        %[r,p] = corr(X, Y);
        %str = {sprintf('corr: %.3f', r), sprintf('pval: %d', p)};
        l = lsline;
        set(l, 'Color', 'r','LineStyle','-','linewidth',2.5);
        %txt=annotation('textbox', [0.14 0.80 0.1 0.1], 'String', str, 'EdgeColor', 'none');
        %set(txt, 'FontSize', 12.5,'fontweight','bold')
    end
         
        xlabel(tistrs(1));        
        ylabel(tistrs(2));     
        title(tistrs(3));
  
end
%}
 %statistics =  cellstr(num2str([Bias;avg;RMSE;Std;Corr;N]))
 %axis([xlim,ylim],'linewidth',2.5,'FontSize', 16,'fontweight','bold')
 %{
x=random('Normal',282,14,5000*4,1);
y=random('Normal',278,12,5000*4,1);
heatscatters(x,y)
%}
