clear;close all;clc;

re = 6.37122e6;% Earth radius

addpath('/Users/xudo627/donghui/mylib/m/');
addpath('/Users/xudo627/donghui/CODE/dengwirda-inpoly-355c57c/');
addpath('/Users/xudo627/donghui/CODE/Setup-E3SM-Mac/matlab-scripts-for-mosart');
addpath('/Users/xudo627/developments/getPanoply_cMap/');
addpath('functions/');

days_of_month = [31;28;31;30;31;30;31;31;30;31;30;31];

plot_dem  = 1;
plot_glwd = 0;
plot_soil = 1;
% Outlet;
[lon_gsim,lat_gsim,S,yr_gsim,mo_gsim,da,mu_gsim,sd,cv] = get_GSIM_discharge('BR_0002835',1);
aream2 = areaint(S.Y,S.X)*4*pi.*(re^2); 

glwd_labels = {'Lake','Reservoir','River','Freshwater Marsh, Floodplain', ...
               'Swamp Forest, Flooded Forest','Coastal Wetland','Pan, Brackish/Saline Wetland',...
               'Bog, Fen, Mire','Intermittent Wetland/Lake','50-100% Wetland',...
               '25-50% Wetland','Wetland Complex (0-25% Wetland)'};
cols = [38	34	142;
        199	41	28 ;
        170	43	246;
        55	129	246;
        183	165	53 ;
        242	175	61 ;
        255	254	84 ;
        141	75	76 ;
        183	246	252;
        70	102	27 ;
        151	172	69 ;
        191	221	144]./255;

if exist('data/glwd_SA.mat','file')
    load('data/glwd_SA.mat');
else
    dx = 360/43200;
    dy = 360/43200;
    data = double(imread('~/projects/land-river-two-way-coupling/Writing/submit/major_revision/glwd/glwd3.tif'));
    data(data == 255) = NaN;
    lon = -180+dx/2 : dx : 180-dx/2;
    lat = 90-dy/2 : -dy : -90+dy/2;
    i1 = min(find(lon >= -82)); i2 = max(find(lon <= -34));
    j1 = min(find(lat <= 10));  j2 = max(find(lat >= -56));

    [x,y] = meshgrid(lon,lat);
    data = data(j1:j2,i1:i2);
    lon  = x(j1:j2,i1:i2);
    lat  = y(j1:j2,i1:i2);
    clear x y;
    save('data/glwd_SA.mat','lon','lat','data');
end

figure; set(gcf,'Position',[10 10 1400 800]);
ax1 = axes('Position',[0.05 0.05 0.425 0.9]);
plot_globalspatial(lon,lat,data);  hold on; colormap(cols); axis equal;
set(colorbar,'visible','off');
set(gca,'xtick',[],'ytick',[]);
colorbar('westoutside','XTickLabel',glwd_labels, 'XTick', [1 : 12], ...
         'FontSize',14,'FontWeight','bold');
%xlim([-82 -34]); ylim([-56 10]);
xlim([-82 -32]); ylim([-56 10]);
caxis([0.5 12.5]);
pos = get(gca,'Position');
plot(S.X,S.Y,'k-','LineWidth',2.5);
set(gca,'Position',[pos(1)+0.2 pos(2) pos(3) pos(4)]);
ax1.Position(1) = 0.2;

in = inpoly2([lon(:) lat(:)],[S.X' S.Y']);
data(~in) = NaN;

ax2 = axes('Position',[0.25 0.04 0.3 0.34]);
plot_globalspatial(lon,lat,data); hold on; colormap(cols); axis equal;
caxis([0.5 12.5]);set(colorbar,'visible','off');
plot(S.X,S.Y,'k-','LineWidth',2.5);
xlim([-65 -53]); ylim([-22 -14]);
s= scatter(lon_gsim,lat_gsim,150,'s','MarkerFaceColor','g','MarkerEdgeColor','k');
set(s,'LineWidth',2);
% ax2.Position(1) = x0;
% ax2.Position(2) = y0;
ax2.XAxis.LineWidth = 2;ax2.XAxis.Color = 'r'; ax2.XAxis.FontSize = 13;ax2.XAxis.FontWeight = 'bold';
ax2.YAxis.LineWidth = 2;ax2.YAxis.Color = 'r'; ax2.YAxis.FontSize = 13;ax2.YAxis.FontWeight = 'bold';

x3  = ax1.Position(1) + ax1.Position(3) + 0.025;
y3  = ax1.Position(2)+ax1.Position(4) - 0.34;
ax3 = axes('Position',[x3 y3 0.3 0.34]);

if plot_dem
    cmap = getPyPlot_cMap('terrain');
    dem = double(imread('/Users/xudo627/projects/wetland_inundation/HYDROSHED_CON_3s/CON_3S_Pantanal.tif'));
    dx = 15/18000;
    dy = 15/18000;
    lon = -65 + dx/2 : dx : -50 -dx/2;
    lat = -10 - dy/2 : -dy : -25 + dy/2;
    [lon,lat] = meshgrid(lon,lat);
    in = inpoly2([lon(:) lat(:)],[S.X' S.Y']);
    dem(~in) = NaN;

    cb3 = plot_globalspatial(lon,lat,dem,1,1); hold on; colormap(ax3,cmap); axis equal;
    plot(S.X,S.Y,'k-','LineWidth',2.5); caxis([0 800]);
    xlim([-65 -53]); ylim([-22 -14]);
    set(gca,'FontSize',13);
    cb3.Position = [x3 + 0.3 + 0.005 y3 0.01 0.34];
end
x4  = ax1.Position(1) + ax1.Position(3) + 0.025;
y4  = ax1.Position(2);
ax4 = axes('Position',[x4 y4 0.3 0.48]);

if plot_soil
end

[tmon_sim,yr_sim,mos_sim,da_sim] = get_time_number(1979,2010);

yr1 = 1998;
yr2 = 2010;

flds = {'trmm','gpcc','mswep','udel','merra','precl','gswp','crun'};
prec = struct([]);
yr_prec = struct([]);
mo_prec = struct([]);
for i = 1 : 6
    disp(['Read ' upper(flds{i})]);
    if exist(['data/' upper(flds{i}) 'in.mat'],'file')
        load(['data/' upper(flds{i}) 'in.mat']);
    else
        [precipin, yr_prec, mo_prec, xc_prec, yc_prec, xv_prec, yv_prec] = get_precip(S,upper(flds{i}));
        save(['data/' upper(flds{i}) 'in.mat'],'precipin', 'yr_prec', 'mo_prec', 'xc_prec', 'yc_prec', 'xv_prec', 'yv_prec');
    end
    ind  = find(yr_prec  >= yr1 & yr_prec <= yr2);
    prec(1).(flds{i}) = nanmean(precipin(:,ind),1).*24;
    yr(1).(flds{i})   = yr_prec;
    mo(1).(flds{i})   = mo_prec;
end
rain1 = load('data/outputs_sub_CLMCRUNCEPv7.mat','rain');
rain2 = load('data/outputs_sub_CLMGSWP3v1.mat','rain');
ind_sim   = find(yr_sim   >= yr1 & yr_sim   <= yr2);
ind_sim2  = find(yr_sim   >= 1979 & yr_sim  <= 2010);
prec(1).crun = nanmean(rain1.rain(:,ind_sim),1).*86400;
prec(1).gswp = nanmean(rain2.rain(:,ind_sim),1).*86400;
pr_crun = nanmean(rain1.rain(:,ind_sim2),1).*86400;

[GLEAMin,GLEAMlon,GLEAMlat,tmon_in2] = get_GLEAM_ET(S,yr1,yr2); % [mm/s]

ind_gsim  = find(yr_gsim  >= yr1 & yr_gsim  <= yr2);
ind_gsim2 = find(yr_gsim  >= 1979 & yr_gsim  <= 2010);

yyaxis left;
bar(nanmean(reshape(pr_crun,[12,32]),2),'FaceAlpha',0.5);
set(gca,'FontSize',13);
ylabel('CRUNCEPv7 Precipitation [mm/day]','FontSize',15,'FontWeight','bold'); hold on;
yyaxis right;
plot(nanmean(reshape(mu_gsim(ind_gsim2),[12,32]),2),'-','LineWidth',2); grid on;
ax4.XTick = 1 : 12;
ax4.XTickLabel = {'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'};
set(gca,'FontSize',14);
ylabel('Discharge [m^{3}/s]','FontSize',15,'FontWeight','bold');
xlim([0.5 12.5]);


add_title(ax1,'(a). GLWD');
add_title(ax3,'(b). DEM [m]');
add_title(ax4,'(c). Precipitation and streamflow seasonality');

cols = distinguishable_colors(10);
figure; set(gcf,'Position',[10 10 1200 1000]);
ax(1) = subplot_tight(2,1,1,[0.06 0.08]);
for i = 1 : length(flds)
    if i == 3
        bar(nanmean(reshape(prec.(flds{i}),[12,yr2-yr1+1]),2),'FaceColor',cols(i,:),'FaceAlpha',0.5);hold on; grid on;
    elseif i == 4
        bar(nanmean(reshape(prec.(flds{i}),[12,yr2-yr1+1]),2),'FaceColor',cols(i,:),'FaceAlpha',0.5);hold on; grid on;
    elseif i == 5
        scatter(1:12,nanmean(reshape(prec.(flds{i}),[12,yr2-yr1+1]),2),200,'d','MarkerFaceColor',cols(i,:),'MarkerEdgeColor','none');
    elseif i == 6
        scatter([1:12]-0.1,nanmean(reshape(prec.(flds{i}),[12,yr2-yr1+1]),2),100,'o','MarkerFaceColor',cols(i,:),'MarkerEdgeColor','none');
    elseif i == 7
        scatter([1:12]+0.1,nanmean(reshape(prec.(flds{i}),[12,yr2-yr1+1]),2),200,'s','MarkerFaceColor',cols(i,:),'MarkerEdgeColor','none');
    elseif i == 8
        plot(nanmean(reshape(prec.(flds{i}),[12,yr2-yr1+1]),2),'--','Color','k','LineWidth',3);hold on; grid on;
    else
        plot(nanmean(reshape(prec.(flds{i}),[12,yr2-yr1+1]),2),'-','Color',cols(i,:),'LineWidth',3);hold on; grid on;
    end
    
    mu_yr(i) = nanmean(prec.(flds{i})).*365.25;
end
leg = legend('TRMM','GPCC','MSWEP','UDEL','MERRA','PRECL','GSWP','CRUN');
leg.FontSize = 16;leg.FontWeight = 'bold';
xlim([0.5 12.5]);
set(gca,'FontSize',13);
ylabel('Precipitation [mm/day]','FontSize',15,'FontWeight','bold');

ax(2) = subplot_tight(2,1,2,[0.06 0.08]);
X = categorical({'TRMM','GPCC','MSWEP','UDEL','MERRA','PRECL','GSWP','CRUN'});
X = reordercats(X,{'TRMM','GPCC','MSWEP','UDEL','MERRA','PRECL','GSWP','CRUN'});
b = bar(X,mu_yr-nanmean(mu_yr)); grid on;
set(gca,'FontSize',13);
ylabel('Discrepency [mm/yrear]','FontSize',15,'FontWeight','bold');

ax(2).XAxis.FontSize = 15;
ax(2).XAxis.FontWeight = 'bold';
discharge = (mu_yr-nanmean(mu_yr))./1000.*aream2./365.25./86400;

ax(1).XTick = 1 : 12;
ax(1).XTickLabel = {'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'};

add_title(ax(1),'(a) Precipitation seasonality');
add_title(ax(2),'(b) Annual precipitation');
% figure;
% for i = 1 : length(flds)
%     tmp = reshape(prec.(flds{i}),[12,yr2-yr1+1]);
%     [~,imax] = max(tmp,[],1);
%     imax(imax == 12) = 0;
%     tmp = NaN(4,1);
%     for j = 0 : 3
%         tmp(j+1) = length(find(imax == j));
%     end
%     disp(nanmean(imax));
%     plot(tmp,'d-','Color',cols(i,:),'LineWidth',2); hold on; grid on;
%     %histogram(imax,'FaceColor',cols(i,:),'FaceAlpha',0.5); hold on; grid on;
% end

figure; set(gcf,'Position',[10 10 1200 1000]);
ax(1) = subplot_tight(2,1,1,[0.06 0.08]);
X = categorical({'TRMM','GPCC','MSWEP','UDEL','MERRA','PRECL','GSWP','CRUN'});
X = reordercats(X,{'TRMM','GPCC','MSWEP','UDEL','MERRA','PRECL','GSWP','CRUN'});
b = bar(X,mu_yr-nanmean(mu_yr)); grid on;
set(gca,'FontSize',13);
ylabel('Deviation [mm/yrear]','FontSize',15,'FontWeight','bold');

ax(1).XAxis.FontSize = 15;
ax(1).XAxis.FontWeight = 'bold';
discharge = (mu_yr-nanmean(mu_yr))./1000.*aream2./365.25./86400;

ax(2) = subplot_tight(2,1,2,[0.06 0.08]);
Pr = NaN(12,8);
for i = 1 : length(flds)
    Pr(:,i) = nanmean(reshape(prec.(flds{i}),[12,yr2-yr1+1]),2);
    %plot(Pr(:,i),'-','Color',[0.5 0.5 0.5],'LineWidth',1);hold on; grid on;
    mu_yr(i) = nanmean(prec.(flds{i})).*365.25;
end
x = 1 : 12;
y1 = min(Pr,[],2)';
y2 = max(Pr,[],2)';
p(1) = patch([x fliplr(x)], [y1 fliplr(y2)], 'g','EdgeColor','none','FaceColor',[0.5 0.5 0.5],'FaceAlpha',0.5); hold on;
for i = 1 : length(flds)
    plot(Pr(:,i),':','Color','k','LineWidth',1);hold on; grid on;
end
p(2) = plot(nanmean(Pr(:,1:7),2),'-','Color','k','LineWidth',3);
p(3) = plot(Pr(:,8),'--','Color','r','LineWidth',3);hold on; grid on;

xlim([1 12]);
ax(2).XTick = 1 : 12;
ax(2).XTickLabel = {'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'};

set(gca,'FontSize',13);
ylabel('Precipitation [mm/day]','FontSize',15,'FontWeight','bold');
leg = legend(p,{'Multi-model ensemble','Multi-model mean','CRUNCEPv7'});
leg.FontSize= 14;
leg.FontWeight = 'bold';

add_title(ax(1),'(a) Annual precipitation',18);
add_title(ax(2),'(b) Precipitation seasonality',18);