clear;close all;clc;

% Test the trend of GRUN runoff in Susquehanna river basin.

re = 6.37122e6;% Earth radius

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

yr1 = 1980;
yr2 = 2012;

%[lon_gsim,lat_gsim,Sh,yr,mo,da,mu_gsim,sd,cv] = get_GSIM_discharge('BR_0002835',1);
[lon_gsim,lat_gsim,Sh,yr,mo,da,mu_gsim,sd,cv] = get_GSIM_discharge('US_0000919',1);%
area_gsim = 4*pi*areaint(Sh.Y,Sh.X).*(re^2);

lon_lora = ncread('/Users/xudo627/DATA/Runoff/LORA.nc','lon');
lat_lora = ncread('/Users/xudo627/DATA/Runoff/LORA.nc','lat');
[lon_lora,lat_lora] = meshgrid(lon_lora,lat_lora);
lon_lora = lon_lora'; lat_lora = lat_lora';
lon_grun = ncread('/Users/xudo627/DATA/Runoff/GRUN_v1_GSWP3_WGS84_05_1902_2014.nc','X');
lat_grun = ncread('/Users/xudo627/DATA/Runoff/GRUN_v1_GSWP3_WGS84_05_1902_2014.nc','Y');
[lon_grun,lat_grun] = meshgrid(lon_grun,lat_grun);
lon_grun = lon_grun'; lat_grun = lat_grun';
load('/Users/xudo627/DATA/Runoff/GRFR.mat');
lon_grfr = lon;
lat_grfr = lat; clear lon lat;

in_grun = inpoly2([lon_grun(:) lat_grun(:)],[Sh.X' Sh.Y']);
in_lora = inpoly2([lon_lora(:) lat_lora(:)],[Sh.X' Sh.Y']);
in_grfr = inpoly2([lon_grfr(:) lat_grfr(:)],[Sh.X' Sh.Y']);

Runoff = ncread('/Users/xudo627/DATA/Runoff/GRUN_v1_GSWP3_WGS84_05_1902_2014.nc','Runoff'); % [mm/day]
mrro   = ncread('/Users/xudo627/DATA/Runoff/LORA.nc','mrro');% [mm/s]
yr_grun = 1902;
yr_lora = 1980;

Runoff = Runoff(:,:,(yr1 - yr_grun)*12+1 : (yr2 - yr_grun + 1) * 12);
mrro   = mrro(:,:,(yr1 - yr_lora)*12+1 : (yr2 - yr_lora + 1) * 12);
GRFR   = GRFR(:,:,(yr1 - yr_grfr)*12+1 : (yr2 - yr_grfr + 1) * 12);

mu_grun = NaN(nansum(in_grun),(yr2-yr1+1)*12);
mu_lora = NaN(nansum(in_lora),(yr2-yr1+1)*12);
mu_grfr = NaN(nansum(in_grfr),(yr2-yr1+1)*12);

for i = 1 : (yr2-yr1+1)*12
    tmp = Runoff(:,:,i);
    mu_grun(:,i) = tmp(in_grun);
    tmp = mrro(:,:,i);
    mu_lora(:,i) = tmp(in_lora);
%     if i <= (yr2-yr1)*12
%     tmp = mrro(:,:,i);
%     mu_lora(:,i) = tmp(in_lora);
%     end
    
    tmp = GRFR(:,:,i);
    mu_grfr(:,i) = tmp(in_grfr);
end

mu_grun = nanmean(mu_grun,1)'./1000.*area_gsim./86400;
mu_lora = nanmean(mu_lora,1)'./1000.*area_gsim;
mu_grfr = nanmean(mu_grfr,1)'./1000.*area_gsim./30.45/86400;

annual_gsim = nanmean(reshape(mu_gsim(yr >= yr1 & yr <= yr2),[12 yr2-yr1+1]),1);
annual_grun = nanmean(reshape(mu_grun,[12 yr2-yr1+1]),1);
annual_lora = nanmean(reshape(mu_lora,[12 yr2-yr1+1]),1);
annual_grfr = nanmean(reshape(mu_grfr,[12 yr2-yr1+1]),1);

%cols = {'k','r','b','g'};
cols = getPyPlot_cMap('Set1',9);
cols = [0 0 0; cols];
models = {'GAGE','GRUN','LORA','GFRF'};
figure; set(gcf,'Position',[10 10 1200 600]);
a1 = (annual_gsim-nanmean(annual_gsim))./nanmean(annual_gsim);
a2 = (annual_grun-nanmean(annual_grun))./nanmean(annual_grun);
a3 = (annual_lora-nanmean(annual_lora))./nanmean(annual_lora);
a4 = (annual_grfr-nanmean(annual_grfr))./nanmean(annual_grfr);

plot(yr1:yr2,annual_gsim,'-','Color',cols(1,:),'LineWidth',3); hold on; grid on;
plot(yr1:yr2,annual_grun,'-','Color',cols(2,:),'LineWidth',3); 
plot(yr1:yr2,annual_lora,'-','Color',cols(3,:),'LineWidth',3); 
plot(yr1:yr2,annual_grfr,'-','Color',cols(4,:),'LineWidth',3); 
set(gca,'FontSize',13);
% leg = legend('Observation','GRUN','LORA','GRFR');
% leg.FontSize = 15;
% leg.FontWeight = 'bold';
ylabel('Annual streamflow [m^{3}/s]','FontSize',16,'FontWeight','bold');

[slope(1), ~] = sens_slope(annual_gsim'); slope(1) = slope(1)/nanmean(annual_gsim)*100;
[z(1), trend(1), ~, pvalue(1)] = mk_test(annual_gsim');
[slope(2), ~] = sens_slope(annual_grun'); slope(2) = slope(2)/nanmean(annual_grun)*100;
[z(2), trend(2), ~, pvalue(2)] = mk_test(annual_grun');
[slope(3), ~] = sens_slope(annual_lora'); slope(3) = slope(3)/nanmean(annual_lora)*100;
[z(3), trend(3), ~, pvalue(3)] = mk_test(annual_lora');
[slope(4), ~] = sens_slope(annual_grfr'); slope(4) = slope(4)/nanmean(annual_grfr)*100;
[z(4), trend(4), ~, pvalue(4)] = mk_test(annual_grfr');
save('data/runoff_uncertainty.mat','annual_grun','annual_lora','annual_grfr','slope','pvalue');
for i = 1 : 4
    sens_slope = num2str(num2str(slope(i)));
    pv         = num2str(pvalue(i));
    if pvalue(i) < 0.001
        pv = ' < 0.05';
    else
        pv = [' = ' pv(1:4)];
    end
    if slope(i) < 0
        sens_slope = sens_slope(1:5);
    else
        sens_slope = ['+' sens_slope(1:4)];
    end
    strs{i} = [models{i} ': Sen slope = ' sens_slope ' [%/yr], p-value ' pv];
    t(i) = annotation('textbox',[0.525 0.850 - i*0.035 0.3 0.1],'String',strs{i},       ... 
                    'FitBoxToText','on','EdgeColor','none','Color',cols(i,:),     ...
                    'FontSize',14,'FontWeight','bold','fontname','Menlo'); 
end
            
% [slope(1), ~] = sens_slope(a1');
% [z(1), trend(1)] = mk_test(a1');
% [slope(2), ~] = sens_slope(a2');
% [z(2), trend(2)] = mk_test(a2');