clear;close all;clc;

addpath('/Users/xudo627/donghui/CODE/dengwirda-inpoly-355c57c/');

datadir = '/Volumes/LaCie/PROJECTS/Pantanal/';

region = 'Susquehanna';

if strcmp(region,'Pantanal')
    numc = 191;
    [lon_gsim,lat_gsim,Sh,yr,mo,da,mu_gsim,sd,cv] = get_GSIM_discharge('BR_0002835',1);
    
    xv=ncread('inputdata_run/domain_lnd_Pantanal_all_template.nc','xv');
    yv=ncread('inputdata_run/domain_lnd_Pantanal_all_template.nc','yv');
    xc=ncread('inputdata_run/domain_lnd_Pantanal_all_template.nc','xc');
    yc=ncread('inputdata_run/domain_lnd_Pantanal_all_template.nc','yc');
else
    numc = 32;
    [lon_gsim,lat_gsim,Sh,yr,mo,da,mu_gsim,sd,cv] = get_GSIM_discharge('US_0000919',1);
    
    xv=ncread('inputdata_run/domain_lnd_Susquehanna_run3_template.nc','xv');
    yv=ncread('inputdata_run/domain_lnd_Susquehanna_run3_template.nc','yv');
    xc=ncread('inputdata_run/domain_lnd_Susquehanna_run3_template.nc','xc');
    yc=ncread('inputdata_run/domain_lnd_Susquehanna_run3_template.nc','yc');
end

%
long = ncread('/Volumes/LaCie/DATA/GSCD_v2/GSCD_v2.0.nc','longitude');
lati = ncread('/Volumes/LaCie/DATA/GSCD_v2/GSCD_v2.0.nc','latitude');
[long,lati] = meshgrid(long,lati);
long = long'; lati = lati';
BFI1 = ncread('/Volumes/LaCie/DATA/GSCD_v2/GSCD_v2.0.nc','BFI1');
BFI2 = ncread('/Volumes/LaCie/DATA/GSCD_v2/GSCD_v2.0.nc','BFI2');
BFI3 = ncread('/Volumes/LaCie/DATA/GSCD_v2/GSCD_v2.0.nc','BFI3');
BFI4 = ncread('/Volumes/LaCie/DATA/GSCD_v2/GSCD_v2.0.nc','BFI4');

BFI1 = BFI1'; BFI2 = BFI2'; BFI3 = BFI3'; BFI4 = BFI4';
in = inpoly2([long(:) lati(:)],[Sh.X' Sh.Y']);
xc2 = long(in); yc2 = lati(in);
xv2 = NaN(4,length(xc2)); yv2 = NaN(4,length(yc2));
dx = 0.05; dy = 0.05;
xv2(1,:) = xc2 - dx/2; xv2(2,:) = xc2 + dx/2; xv2(3,:) = xc2 + dx/2; xv2(4,:) = xc2 - dx/2;
yv2(1,:) = yc2 - dy/2; yv2(2,:) = yc2 - dy/2; yv2(3,:) = yc2 + dy/2; yv2(4,:) = yc2 + dy/2;


figure;
subplot(2,2,1);
patch(xv2,yv2,BFI1(in),'LineStyle','none'); colormap(jet(100)); caxis([0 1]); hold on;
plot(Sh.X,Sh.Y,'k-','LineWidth',3);
title('BFI1','FontSize',15,'FontWeight','bold');
subplot(2,2,2);
patch(xv2,yv2,BFI2(in),'LineStyle','none');colormap(jet(100)); caxis([0 1]); hold on;
plot(Sh.X,Sh.Y,'k-','LineWidth',3);
title('BFI2','FontSize',15,'FontWeight','bold');
subplot(2,2,3);
patch(xv2,yv2,BFI3(in),'LineStyle','none'); colormap(jet(100)); caxis([0 1]); hold on;
plot(Sh.X,Sh.Y,'k-','LineWidth',3);
title('BFI3','FontSize',15,'FontWeight','bold');
subplot(2,2,4);
patch(xv2,yv2,BFI4(in),'LineStyle','none'); colormap(jet(100)); caxis([0 1]); hold on;
plot(Sh.X,Sh.Y,'k-','LineWidth',3);
title('BFI4','FontSize',15,'FontWeight','bold');

BFI3_half = convert_res(double(BFI3),1,10)./100;
lon_half  = convert_res(double(long),1,10)./100;
lat_half  = convert_res(double(lati),1,10)./100;

BFI = NaN(length(xc),1);
for i = 1 : length(xc)
    ind = find(lon_half == xc(i) & lat_half == yc(i));
    disp(ind);
    BFI(i) = BFI3_half(ind);
end

if strcmp(region,'Susquehanna')
    save('SUS_BFI3.mat','BFI');
else
    save('BFI3.mat','BFI');
end

% tag = 'outputs_geo_CLMCRUNCEPv7';
% rng(1);
% if strcmp(tag,'outputs_sur_CLMCRUNCEPv7')
%     ntot = 500;
%     xtrain(:,1)= ones(ntot,1).*2.5;
%     xtrain(:,2)= 0.1 + (5-0.1)*rand(ntot,1);
%     xtrain(:,3)= 0.001 + (1-0.001)*rand(ntot,1); % micro_sigma
%     xtrain(:,4)= 0.01 + (1-0.01)*rand(ntot,1); % fc  
% else
%     ntot = 100;
%     xtrain(:,1) = 0.1 + (5-0.1)*rand(100,1);
% end
% % Get runoff from default parameters
% % Read MOSART
% files = dir([datadir tag '/Pantanal_half*elm.h0.*.nc']);
% 
% qover   = NaN(numc,length(files));
% qdrai   = NaN(numc,length(files));
% qrunoff = NaN(numc,length(files));
% qdraiall   = NaN(ntot,length(files));
% qrunoffall = NaN(ntot,length(files));
% 
% ind = 91;
% for i = 1 : length(files)
%     filename = fullfile(files(i).folder,files(i).name);
%     disp(filename);
%     tmp1 = ncread(filename,'FH2OSFC');
%     tmp2 = ncread(filename,'QOVER');
%     tmp3 = ncread(filename,'QDRAI');
%     tmp4 = ncread(filename,'QRUNOFF');
%     
%     qover(:,i)   = tmp2((ind-1)*numc+1:ind*numc);
%     qdrai(:,i)   = tmp3((ind-1)*numc+1:ind*numc);
%     qrunoff(:,i) = tmp4((ind-1)*numc+1:ind*numc);
%     
%     for j = 1 : ntot
%         qdraiall(j,i) = nanmean(tmp3((j-1)*numc+1:j*numc));
%         qrunoffall(j,i) = nanmean(tmp4((j-1)*numc+1:j*numc));
%     end
% end
% 
% figure;
% patch(xv,yv,nanmean(qdrai./qrunoff,2),'LineStyle','none'); colorbar;
% colormap(jet(100)); caxis([0 1]); hold on;
% plot(Sh.X,Sh.Y,'k-','LineWidth',3);
% 
% BFI = nanmean(qdraiall,2)./nanmean(qrunoffall,2);
% for i = 1 : size(xtrain,2)
%     subplot(size(xtrain,2),1,i);
%     plot(xtrain(:,i),BFI,'bo');
% end