clc
clear all
close all

cn51_read_path = '/discover/nobackup/projects/geoscm/fzeng/Catchment-CN5.1/GEOSldas_CN51/output/SMAP_EASEv2_M36/cat/ens0000/';
clm45_read_path = '/discover/nobackup/projects/geoscm/fzeng/GEOSldas45_M36_rst_ldas06_15/output/SMAP_EASEv2_M36/cat/ens0000/';
cn40_read_path = '/discover/nobackup/projects/geoscm/fzeng/Catchment-CN40_9km/GEOSldas_1981_present/GEOSldas_CNCLM40_9km/output/SMAP_EASEv2_M09_GLOBAL/cat/ens0000/';
fluxcom_read_path = '/discover/nobackup/projects/geoscm/fzeng/Fluxcom/';
alexi_read_path = '/discover/nobackup/projects/geoscm/fzeng/ALEXI/monthly/';

start_year = 2001;
stop_year = 2019; 

leapYears = 1980:4:2022;
monthLength_nly = [31 28 31 30 31 30 31 31 30 31 30 31];
monthLength_ly = [31 29 31 30 31 30 31 31 30 31 30 31];

% compute mapping indices

cn51_read_file = [cn51_read_path 'Y2001/M01/GEOSldas_CN51.tavg24_1d_lnd_Nt.monthly.200101.nc4'];
clm45_read_file = [clm45_read_path 'Y2001/M01/GEOSldas45_M36_rst_ldas06_15.tavg24_1d_lnd_Nt.monthly.200101.nc4'];
cn40_read_file = [cn40_read_path 'Y2001/M01/GEOSldas_CNCLM40_9km.tavg24_1d_lnd_Nt.monthly.200101.nc4'];

fluxcom_et_read_file = [fluxcom_read_path 'LE.RS.EBC-ALL.MLM-ALL.METEO-NONE.720_360.monthly.2001.nc'];
alexi_et_read_file = [alexi_read_path 'ALEXI_ET_201501.nc'];

lat_vec = double(ncread(cn51_read_file,'lat'));
lon_vec = double(ncread(cn51_read_file,'lon'));

clm45_lat_vec = double(ncread(clm45_read_file,'lat'));
clm45_lon_vec = double(ncread(clm45_read_file,'lon'));

cn40_lat_vec = double(ncread(clm45_read_file,'lat'));
cn40_lon_vec = double(ncread(clm45_read_file,'lon'));

fc_lat_vec = ncread(fluxcom_et_read_file,'lat');
fc_lon_vec = ncread(fluxcom_et_read_file,'lon');

alexi_lat_vec = ncread(alexi_et_read_file,'latitude');
alexi_lon_vec = ncread(alexi_et_read_file,'longitude');

for i = 1:length(lat_vec)
    lat_dist = sqrt((fc_lat_vec-lat_vec(i)).^2);
    fc_lat_index(i) = find(lat_dist==min(lat_dist));

    lon_dist = sqrt((fc_lon_vec-lon_vec(i)).^2);
    tmp = find(lon_dist==min(lon_dist));
    fc_lon_index(i) = tmp(1);

    lat_dist = sqrt((alexi_lat_vec-lat_vec(i)).^2);
    lon_dist = sqrt((alexi_lon_vec-lon_vec(i)).^2);
    if (min(lat_dist)<=0.25) & (min(lon_dist)<=0.25)
       alexi_lat_index(i) = find(lat_dist==min(lat_dist));
       alexi_lon_index(i) = find(lon_dist==min(lon_dist));
    else
       alexi_lat_index(i) = 0;
       alexi_lon_index(i) = 0;
    end


    dist_vec = sqrt((clm45_lat_vec-lat_vec(i)).^2 + (clm45_lon_vec-lon_vec(i)).^2);
    min_index = find(dist_vec==min(dist_vec));
    data_index(i) = min_index(1);

    dist_vec = sqrt((cn40_lat_vec-lat_vec(i)).^2 + (cn40_lon_vec-lon_vec(i)).^2);
    min_index = find(dist_vec==min(dist_vec));
    cn40_data_index(i) = min_index(1);

end

% read data

mat_index = 0;
y_index = 0;

for y = start_year:stop_year

    y_index = y_index + 1;
    y_str = num2str(y,'%04i');

    % define month length based on leap year status
    if sum(y==leapYears)>0
        monthLength = monthLength_ly;
    else
        monthLength = monthLength_nly;
    end

    % loop over months
    for m = 1:12

        mat_index = mat_index + 1;

        m_str = num2str(m,'%02i');

        cn51_read_file = [cn51_read_path 'Y' y_str '/M' m_str '/GEOSldas_CN51.tavg24_1d_lnd_Nt.monthly.' y_str m_str '.nc4'];
        clm45_read_file = [clm45_read_path 'Y' y_str '/M' m_str '/GEOSldas45_M36_rst_ldas06_15.tavg24_1d_lnd_Nt.monthly.' y_str m_str '.nc4'];
        cn40_read_file = [cn40_read_path 'Y' y_str '/M' m_str '/GEOSldas_CNCLM40_9km.tavg24_1d_lnd_Nt.monthly.' y_str m_str '.nc4'];
        fluxcom_et_read_file = [fluxcom_read_path 'LE.RS.EBC-ALL.MLM-ALL.METEO-NONE.720_360.monthly.' y_str '.nc'];
        fluxcom_rn_read_file = [fluxcom_read_path 'Rn.RS.EBC-NONE.MLM-ALL.METEO-NONE.720_360.monthly.' y_str '.nc'];

        alexi_et_read_file = [alexi_read_path 'ALEXI_ET_' y_str m_str '.nc'];

        cn51_forc_file = [cn51_read_path 'Y' y_str '/M' m_str '/GEOSldas_CN51.tavg24_1d_lfs_Nt.monthly.' y_str m_str '.nc4'];

        if (mat_index == 1)

           cn51_et_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           cn45_et_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           cn40_et_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           fc_et_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           alexi_et_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           cn51_precip_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);

        end

        evpsoil_tmp_mat = double(ncread(cn51_read_file,'EVPSOIL'));
        evptrns_tmp_mat = double(ncread(cn51_read_file,'EVPTRNS'));
        evpintr_tmp_mat = double(ncread(cn51_read_file,'EVPINTR'));
        evpsbln_tmp_mat = double(ncread(cn51_read_file,'EVPSBLN'));
        cn51_et_mat(:,mat_index) = evpsoil_tmp_mat + evptrns_tmp_mat + ...
                                   evpintr_tmp_mat + evpsbln_tmp_mat;

        cn51_precip_mat(:,mat_index) = double(ncread(cn51_forc_file,'RainfSnowf'));

        if (y<=2016)
           evpsoil_tmp_mat = double(ncread(clm45_read_file,'EVPSOIL'));
           evptrns_tmp_mat = double(ncread(clm45_read_file,'EVPTRNS'));
           evpintr_tmp_mat = double(ncread(clm45_read_file,'EVPINTR'));
           evpsbln_tmp_mat = double(ncread(clm45_read_file,'EVPSBLN'));
           cn45_et_mat(:,mat_index) = evpsoil_tmp_mat(data_index) + evptrns_tmp_mat(data_index) + ... 
                                   evpintr_tmp_mat(data_index) + evpsbln_tmp_mat(data_index);
        end

        evpsoil_tmp_mat = double(ncread(cn40_read_file,'EVPSOIL'));
        evptrns_tmp_mat = double(ncread(cn40_read_file,'EVPTRNS'));
        evpintr_tmp_mat = double(ncread(cn40_read_file,'EVPINTR'));
        evpsbln_tmp_mat = double(ncread(cn40_read_file,'EVPSBLN'));
        cn40_et_mat(:,mat_index) = evpsoil_tmp_mat(cn40_data_index) + evptrns_tmp_mat(cn40_data_index) + ... 
                                   evpintr_tmp_mat(cn40_data_index) + evpsbln_tmp_mat(cn40_data_index);

        if (y>=2015)
           alexi_et_tmp_mat = ncread(alexi_et_read_file,'ET');

           for i = 1:length(lat_vec)
               if (alexi_lat_index(i)~=0)
                  alexi_et_mat(i,mat_index) = alexi_et_tmp_mat(alexi_lat_index(i),alexi_lon_index(i));
               end
           end
        end

        datenum_vec(mat_index) = datenum(y,m,1);
        year_vec(mat_index) = y;
        month_vec(mat_index) = m;

    end

    if (y<=2015)

       et_tmp_mat = ncread(fluxcom_et_read_file,'LE');
       et_tmp_mat(et_tmp_mat==-9999) = NaN;

       rn_tmp_mat = ncread(fluxcom_rn_read_file,'Rn');
       rn_tmp_mat(rn_tmp_mat==-9999) = NaN;

       for i = 1:length(lat_vec)

           start_index = (y_index-1)*12 + 1;
           stop_index = y_index*12;

           fc_et_mat(i,start_index:stop_index) = et_tmp_mat(fc_lon_index(i),fc_lat_index(i),:);
           fc_rn_mat(i,start_index:stop_index) = rn_tmp_mat(fc_lon_index(i),fc_lat_index(i),:);

       end
    end
end
       

% remove extra entries

cn51_et_mat(:,(mat_index+1):end) = [];
cn45_et_mat(:,(mat_index+1):end) = [];
cn40_et_mat(:,(mat_index+1):end) = [];
fc_et_mat(:,(mat_index+1):end) = [];
alexi_et_mat(:,(mat_index+1):end) = [];
cn51_precip_mat(:,(mat_index+1):end) = [];

% convert Fluxcom variables from MJ m^-2 day^-1 to W m^-2

fc_et_mat = fc_et_mat./86400;
fc_et_mat = fc_et_mat*1000000;
fc_rn_mat = fc_rn_mat./86400;
fc_rn_mat = fc_rn_mat*1000000;

% convert ALEXI ET data from mm/day to W m^-2

alexi_et_mm_mat = alexi_et_mat;
alexi_et_mat = alexi_et_mat./(0.0864*0.408);

% convert CN51 precipitation data to mm/day and create ET matrix in mm/day

cn51_precip_mm_mat = cn51_precip_mat * 86400;
cn51_et_mm_mat = cn51_et_mat*0.0864*0.408;
cn45_et_mm_mat = cn45_et_mat*0.0864*0.408;
fc_et_mm_mat = fc_et_mat*0.0864*0.408;

ann_index = 0;

for y = start_year:stop_year

    ann_index = ann_index + 1;
    year_index = (year_vec==y);

    tmp = cn51_et_mat(:,year_index);
    mean_global_ann_cn51_et(ann_index) = nanmean(tmp(:));
    mean_local_ann_cn51_et(:,ann_index) = nanmean(tmp,2);

    tmp = cn40_et_mat(:,year_index);
    mean_global_ann_cn40_et(ann_index) = nanmean(tmp(:));
    mean_local_ann_cn40_et(:,ann_index) = nanmean(tmp,2);

    tmp = fc_et_mat(:,year_index);
    mean_global_ann_fc_et(ann_index) = nanmean(tmp(:));
    mean_local_ann_fc_et(:,ann_index) = nanmean(tmp,2);

    tmp = alexi_et_mat(:,year_index);
    mean_global_ann_alexi_et(ann_index) = nanmean(tmp(:));
    mean_local_ann_fs_alexi_et(:,ann_index) = nanmean(tmp,2);

    tmp = cn45_et_mat(:,year_index);
    mean_global_ann_cn45_et(ann_index) = nanmean(tmp(:));
    mean_local_ann_cn45_et(:,ann_index) = nanmean(tmp,2);

end

figure
plot(start_year:stop_year,mean_global_ann_cn51_et,'LineWidth',2);
hold on
plot(start_year:stop_year,mean_global_ann_cn45_et,'LineWidth',2);
hold on
%plot(start_year:stop_year,mean_global_ann_cn40_et,'LineWidth',2);
%hold on
plot(start_year:stop_year,mean_global_ann_fc_et,'LineWidth',2,'LineStyle','--');
hold on
plot(start_year:stop_year,mean_global_ann_alexi_et,'LineWidth',2,'LineStyle','--');
hold off
grid on
xlabel('year')
ylabel('ET [W m^{-2}]')
title('annual average ET [W m^{-2}]')
legend('CN51','CN45','FluxCom','ALEXI','Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/CN51_eval_Fluxcom_ET'],'png')

% compute and plot global difference maps

diff_et_cn51_fc = cn51_et_mat - fc_et_mat;
mean_diff_et_cn51_fc = nanmean(diff_et_cn51_fc,2);
tmp = (fc_et_mat-cn51_et_mat).^2;
mse_et_cn51_fc = nanmean(tmp,2);
rmse_et_cn51_fc = sqrt(mse_et_cn51_fc);
emp_cn51_fc = fc_et_mm_mat - cn51_precip_mm_mat;
emp_cn51_fc = nanmean(emp_cn51_fc,2);

diff_et_cn51_alexi = cn51_et_mat - alexi_et_mat;
mean_diff_et_cn51_alexi = nanmean(diff_et_cn51_alexi,2);
tmp = (alexi_et_mat-cn51_et_mat).^2;
mse_et_cn51_alexi = nanmean(tmp,2);
rmse_et_cn51_alexi = sqrt(mse_et_cn51_alexi);
emp_cn51_alexi = alexi_et_mm_mat - cn51_precip_mm_mat;
emp_cn51_alexi = nanmean(emp_cn51_alexi,2);

diff_et_cn45_fc = cn45_et_mat - fc_et_mat;
mean_diff_et_cn45_fc = nanmean(diff_et_cn45_fc,2);
tmp = (fc_et_mat-cn45_et_mat).^2;
mse_et_cn45_fc = nanmean(tmp,2);
rmse_et_cn45_fc = sqrt(mse_et_cn45_fc);

diff_et_cn45_alexi = cn45_et_mat - alexi_et_mat;
mean_diff_et_cn45_alexi = nanmean(diff_et_cn45_alexi,2);
tmp = (alexi_et_mat-cn45_et_mat).^2;
mse_et_cn45_alexi = nanmean(tmp,2);
rmse_et_cn45_alexi = sqrt(mse_et_cn45_alexi);

closeness_et_fc = rmse_et_cn51_fc - rmse_et_cn45_fc;
closeness_et_alexi = rmse_et_cn51_alexi - rmse_et_cn45_alexi;

test_lat = unique(lat_vec);
test_lon = unique(lon_vec);

[lon_mat,lat_mat] = meshgrid(test_lon,test_lat);

mean_diff_et_cn51_fc_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
mean_diff_et_cn45_fc_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));

rmse_et_cn51_fc_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
rmse_et_cn45_fc_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
closeness_et_fc_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
closeness_et_r_fc_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
r_et_cn51_fc_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
r_et_cn45_fc_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));

rmse_et_cn51_alexi_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
rmse_et_cn45_alexi_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
closeness_et_alexi_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
closeness_et_r_alexi_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
r_et_cn51_alexi_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
r_et_cn45_alexi_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));

emp_cn51_fc_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
emp_cn51_alexi_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));

for i = 1:length(lon_vec)

    lat_index = test_lat == lat_vec(i);
    lon_index = test_lon == lon_vec(i);

    mean_diff_et_cn51_fc_mat(lat_index,lon_index) = mean_diff_et_cn51_fc(i);
    mean_diff_et_cn45_fc_mat(lat_index,lon_index) = mean_diff_et_cn45_fc(i);

    rmse_et_cn51_fc_mat(lat_index,lon_index) = rmse_et_cn51_fc(i);
    rmse_et_cn45_fc_mat(lat_index,lon_index) = rmse_et_cn45_fc(i);
    closeness_et_fc_mat(lat_index,lon_index) = closeness_et_fc(i);

    rmse_et_cn51_alexi_mat(lat_index,lon_index) = rmse_et_cn51_alexi(i);
    rmse_et_cn45_alexi_mat(lat_index,lon_index) = rmse_et_cn45_alexi(i);
    closeness_et_alexi_mat(lat_index,lon_index) = closeness_et_alexi(i);

    emp_cn51_fc_mat(lat_index,lon_index) = emp_cn51_fc(i);
    emp_cn51_alexi_mat(lat_index,lon_index) = emp_cn51_alexi(i);

    % compute correlations

    data_index = isnan(cn51_et_mat(i,:)) + isnan(fc_et_mat(i,:));
    data_index = data_index == 0;
    if sum(data_index)>10
       tmp = corrcoef(cn51_et_mat(i,data_index),fc_et_mat(i,data_index));
       r_et_cn51_fc_mat(lat_index,lon_index) = tmp(1,2);
    end

    data_index = isnan(cn45_et_mat(i,:)) + isnan(fc_et_mat(i,:));
    data_index = data_index == 0;
    if sum(data_index)>10
       tmp = corrcoef(cn45_et_mat(i,data_index),fc_et_mat(i,data_index));
       r_et_cn45_fc_mat(lat_index,lon_index) = tmp(1,2);
    end

    data_index = isnan(cn51_et_mat(i,:)) + isnan(alexi_et_mat(i,:));
    data_index = data_index == 0;
    if sum(data_index)>10
       tmp = corrcoef(cn51_et_mat(i,data_index),alexi_et_mat(i,data_index));
       r_et_cn51_alexi_mat(lat_index,lon_index) = tmp(1,2);
    end 

    data_index = isnan(cn45_et_mat(i,:)) + isnan(alexi_et_mat(i,:));
    data_index = data_index == 0;
    if sum(data_index)>10
       tmp = corrcoef(cn45_et_mat(i,data_index),alexi_et_mat(i,data_index));
       r_et_cn45_alexi_mat(lat_index,lon_index) = tmp(1,2);
    end 

end

% compute correlation-based closeness

closeness_et_r_fc_mat = r_et_cn51_fc_mat - r_et_cn45_fc_mat;
closeness_et_r_alexi_mat = r_et_cn51_alexi_mat - r_et_cn45_alexi_mat;

% plots
 
load coastlines;
[latcells, loncells] = polysplit(coastlat, coastlon);

latlim = [-60 90];
lonlim = [-180 180];

cmap_mat = cmocean('balance',12);
%cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
cmap = cmap_mat; 
s1 = quantile(mean_diff_et_cn51_fc_mat(:),0.01); 
s2 = quantile(mean_diff_et_cn51_fc_mat(:),0.99);
caxis_min = -1*max(abs(s1),abs(s2));
caxis_max = max(abs(s1),abs(s2));

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,mean_diff_et_cn51_fc_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['mean diff. ET (CN51 - FluxCom) [W m^{-2}]'])
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/diff_ET_CN51_FC'],'png')

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on; 
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,mean_diff_et_cn45_fc_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['mean diff. ET (CN45 - FluxCom) [W m^{-2}]'])
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/diff_ET_CN45_FC'],'png')

%------------------
% RMSE plots
%------------------

cmap_mat = cbrewer2('YlOrRd',12); 
%cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
cmap = cmap_mat; 
s1 = quantile(rmse_et_cn51_fc_mat(:),0.01); 
s2 = quantile(rmse_et_cn51_fc_mat(:),0.99);
caxis_min = s1;
caxis_max = s2;

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,rmse_et_cn51_fc_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['CN51 vs. FluxCom RMSE [W m^{-2}]; mean: ' num2str(nanmean(rmse_et_cn51_fc_mat(:)))])
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/RMSE_ET_CN51_FC'],'png')

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,rmse_et_cn45_fc_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['CN45 vs. FluxCom RMSE [W m^{-2}]; mean: ' num2str(nanmean(rmse_et_cn45_fc_mat(:)))])
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/RMSE_ET_CN45_FC'],'png')

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,rmse_et_cn51_alexi_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['CN51 vs. ALEXI RMSE [W m^{-2}]; mean: ' num2str(nanmean(rmse_et_cn51_alexi_mat(:)))])
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/RMSE_ET_CN51_ALEXI'],'png')

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,rmse_et_cn45_alexi_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['CN45 vs. ALEXI RMSE [W m^{-2}]; mean: ' num2str(nanmean(rmse_et_cn45_alexi_mat(:)))])
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/RMSE_ET_CN45_ALEXI'],'png')


%-----------------------------
% closeness plots (RMSE-based)
%-----------------------------

cmap_mat = cmocean('balance',12);
%cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
cmap = cmap_mat;
s1 = quantile(closeness_et_fc_mat(:),0.01);
s2 = quantile(closeness_et_fc_mat(:),0.99);
caxis_min = -1*max(abs(s1),abs(s2));
caxis_max = max(abs(s1),abs(s2));

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,closeness_et_fc_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['closeness to FluxCom (RMSE-based) [W m^{-2}] '])
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/closeness_ET_CN51_CN45_FC'],'png')

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,closeness_et_alexi_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['closeness to ALEXI (RMSE-based) [W m^{-2}] '])
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/closeness_ET_CN51_CN45_ALEXI'],'png')

%-----------------------------
% closeness plots (R-based)
%-----------------------------

cmap_mat = cmocean('balance',12);
cmap_mat = flipud(cmap_mat);
%cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
cmap = cmap_mat;
s1 = quantile(closeness_et_r_fc_mat(:),0.01);
s2 = quantile(closeness_et_r_fc_mat(:),0.99);
caxis_min = -1*max(abs(s1),abs(s2));
caxis_max = max(abs(s1),abs(s2));

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,closeness_et_r_fc_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['closeness to FluxCom (R-based) [-] '])
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/closeness_ET_R_CN51_CN45_FC'],'png')

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,closeness_et_r_alexi_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')      
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);    
set(gca,'YLim',latlim);
tightmap
title(['closeness to ALEXI (R-based) [-] '])
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/closeness_ET_R_CN51_CN45_ALEXI'],'png')

%-----------------------------
% EmP plots
%-----------------------------

cmap_mat = cmocean('deep',12);
cmap_mat = flipud(cmap_mat);
%cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
cmap = cmap_mat;
s1 = quantile(emp_cn51_fc_mat(:),0.01);
s2 = quantile(emp_cn51_fc_mat(:),0.99);
caxis_min = s1;
caxis_max = s2;

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])

set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['Fluxcom ET - CN51 P [mm/day] '])
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/EmP_CN51_FC'],'png')

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,emp_cn51_alexi_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['ALEXI ET - CN51 P [mm/day] '])
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/EmP_CN51_ALEXI'],'png')


