clc
clear all 
close all

usgs_sf_read_path = '/gpfsm/dnb05/projects/p15/iau/merra_land/DATA/USGS_StreamFlow/';
cn51_read_path = '/gpfsm/dnb05/projects/p51/SMAP_Nature/SMAP_Nature_v10/SMAP_Nature_v10.0_M36/output/SMAP_EASEv2_M36_GLOBAL/cat/ens0000/';

start_year = 2015;
stop_year = 2019;
start_month = 1;
stop_month = 4;

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];

% read station information

station_info_file = [usgs_sf_read_path 'Riverflow_Station_Information_NA-M36.nc'];

station_lat = double(ncread(station_info_file,'sta_lat'));
station_lon = double(ncread(station_info_file,'sta_lon'));
comp_basin_area = double(ncread(station_info_file,'CompBasinArea'));  % km^2
%comp_basin_area = comp_basin_area*1000000;         % m^2
station_name = double(ncread(station_info_file,'Station_Name'));
station_ID = double(ncread(station_info_file,'Station_ID'));
n_catch_basin = double(ncread(station_info_file,'N_catch_basin'));
global_id = double(ncread(station_info_file,'GlobalID'));
local_id = double(ncread(station_info_file,'LocalID'));
n_cells_basin = double(ncread(station_info_file,'N_cells_basin'));
smap_id = double(ncread(station_info_file,'SMAPID'));
smap_frac = double(ncread(station_info_file,'SMAPFrac'));
smap_i = double(ncread(station_info_file,'I_SMAP'));
smap_j = double(ncread(station_info_file,'J_SMAP'));

for i = 1:length(station_lat)

    tmp_str = [num2str((station_ID(1,i)-48)) ...
                      num2str((station_ID(2,i)-48)) ...
                      num2str((station_ID(3,i)-48)) ...
                      num2str((station_ID(4,i)-48)) ...
                      num2str((station_ID(5,i)-48)) ...
                      num2str((station_ID(6,i)-48)) ...
                      num2str((station_ID(7,i)-48)) ...
                      num2str((station_ID(8,i)-48))];
    station_ID_new{i} = tmp_str;
end

% read model data

mat_index = 0;
streamflow_station_mat_30d = [];
streamflow_station_mat_15d = [];
streamflow_station_mat_10d = [];

for y = start_year:stop_year

    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

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

        for d = 1:monthLength(m)

            mat_index = mat_index + 1;

            d_str = num2str(d,'%02i');
   
            if (y<2015)
            cn51_read_file = [cn51_read_path 'Y' y_str '/M' m_str '/SMAP_Nature_v10.0_M36.SMAP_L4_SM_gph.' y_str m_str d_str '_1200z.nc4'];
            else
            cn51_read_file = [cn51_read_path 'Y' y_str '/M' m_str '/SMAP_Nature_v10.0_M36.SMAP_L4_SM_gph.' y_str m_str d_str '.nc4'];
            end

            if (mat_index==1)

               cn51_lat_vec = double(ncread(cn51_read_file,'lat'));
               cn51_lon_vec = double(ncread(cn51_read_file,'lon'));

               cn51_streamflow_30d_mat = NaN*ones(length(station_lat),(stop_year-start_year+1)*12*31);
               cn51_streamflow_15d_mat = NaN*ones(length(station_lat),(stop_year-start_year+1)*12*31);
               cn51_streamflow_10d_mat = NaN*ones(length(station_lat),(stop_year-start_year+1)*12*31);

            end
             
            runoff_tmp_mat = double(ncread(cn51_read_file,'overland_runoff_flux'));
            baseflow_tmp_mat = double(ncread(cn51_read_file,'baseflow_flux'));
 
            streamflow_station_vec = NaN*ones(size(smap_id,2),1);

            for i = 1:size(smap_id,2)
                data_index = smap_id(:,i);
                data_frac = smap_frac(:,i);
                data_frac(data_index<0) = [];
                data_index(data_index<0) = [];
                streamflow_tmp_vec = runoff_tmp_mat(data_index) + baseflow_tmp_mat(data_index); % streamflow in each tile
                streamflow_tmp_vec = streamflow_tmp_vec.*data_frac; % streamflow in each tile weighted by fraction of tile in basin
                streamflow_sum = nansum(streamflow_tmp_vec)./sum(data_frac); % basin-wide average streamflow in kg m^-2 s^-1
                %streamflow_sum = streamflow_sum*comp_basin_area(i); % basin-wide average streamflow in kg s^-1
                streamflow_station_vec(i) = streamflow_sum;
            end

            streamflow_station_mat_30d = [streamflow_station_mat_30d streamflow_station_vec];
            streamflow_station_mat_15d = [streamflow_station_mat_15d streamflow_station_vec];
            streamflow_station_mat_10d = [streamflow_station_mat_10d streamflow_station_vec];

            if (size(streamflow_station_mat_30d,2)>30) % only retain previous 30 days for temporal averaging
               streamflow_station_mat_30d(:,1) = [];
            end

            if (size(streamflow_station_mat_15d,2)>15) % only retain previous 15 days for temporal averaging
               streamflow_station_mat_15d(:,1) = [];
            end

            if (size(streamflow_station_mat_10d,2)>10) % only retain previous 10 days for temporal averaging
               streamflow_station_mat_10d(:,1) = [];
            end

            cn51_streamflow_30d_mat(:,mat_index) = nanmean(streamflow_station_mat_30d,2);
            cn51_streamflow_15d_mat(:,mat_index) = nanmean(streamflow_station_mat_15d,2);
            cn51_streamflow_10d_mat(:,mat_index) = nanmean(streamflow_station_mat_10d,2);
            cn51_datetime(:,mat_index) = datetime(y,m,d);

        end
     end
end

cn51_streamflow_30d_mat(:,(mat_index+1):end) = [];
cn51_streamflow_15d_mat(:,(mat_index+1):end) = [];
cn51_streamflow_10d_mat(:,(mat_index+1):end) = [];

% convert model streamflow from kg m^-2 s^-1 to mm/day

cn51_streamflow_30d_mat = cn51_streamflow_30d_mat*86400;
cn51_streamflow_15d_mat = cn51_streamflow_15d_mat*86400;
cn51_streamflow_10d_mat = cn51_streamflow_10d_mat*86400;

% read station data

usgs_data_dir = [usgs_sf_read_path '2022-11-18/'];
usgs_sf_files = dir(usgs_data_dir);
usgs_sf_files = usgs_sf_files(3:end);

cf2cm = 0.0283168466; % ft^3 to m^3

for i = 1:length(usgs_sf_files)

    read_file = [usgs_data_dir usgs_sf_files(i).name];
 
    tmp = detectImportOptions(read_file);
    header_lines = tmp.DataLines(1)-1;

    fid = fopen(read_file);
   % data = textscan(fid,'%4c %8c %{yyyy-mm-dd}D %f %s','HeaderLines',31);
    %data = textscan(fid,'%4c %8c %s %f %s','HeaderLines',header_lines);
    data = textscan(fid,'%4c %8c %s %f %s %f %s %f %s','HeaderLines',header_lines);
    fclose(fid);

    usgs_station_id = data{2};
    usgs_datetime = data{3};
    usgs_year = NaN*ones(length(usgs_datetime),1);
    usgs_month = NaN*ones(length(usgs_datetime),1);
    usgs_day = NaN*ones(length(usgs_datetime),1);
    for j = 1:length(usgs_datetime)
        usgs_year(j) = str2double(usgs_datetime{j}(1:4));
        usgs_month(j) = str2double(usgs_datetime{j}(6:7));
        usgs_day(j) = str2double(usgs_datetime{j}(9:10));
    end
    usgs_sf_data = data{4};

    for j = 1:length(station_ID_new)
        if strcmp(station_ID_new{j},usgs_station_id(1,:))
           cn51_mat_index = j;
           plot_station_id{i} = station_ID_new{j};
           plot_station_lat(i) = station_lat(j);
           plot_station_lon(i) = station_lon(j);
        end 
    end

    % convert USGS data from ft^3 s^-1 to mm/day

    usgs_sf_data = usgs_sf_data * cf2cm *86.4 / comp_basin_area(cn51_mat_index);

    start_index = 0;
    stop_index = 0;

    start_index =  (usgs_year==cn51_datetime(1).Year) + ...
                  (usgs_month==cn51_datetime(1).Month) + ...
                  (usgs_day==cn51_datetime(1).Day);
    start_index = find(start_index==3);
    stop_index =  (usgs_year==cn51_datetime(end).Year) + ...
                  (usgs_month==cn51_datetime(end).Month) + ...
                  (usgs_day==cn51_datetime(end).Day);
    stop_index = find(stop_index==3);

    if (isempty(start_index) | isempty(stop_index))
       continue
    end

    usgs_sf_tmp_mat = usgs_sf_data(start_index:stop_index);
    cn51_sf_30d_tmp_mat = cn51_streamflow_30d_mat(cn51_mat_index,:);
    cn51_sf_15d_tmp_mat = cn51_streamflow_15d_mat(cn51_mat_index,:);
    cn51_sf_10d_tmp_mat = cn51_streamflow_10d_mat(cn51_mat_index,:);

    % compute metrics based on daily data

    r_tmp = corrcoef(usgs_sf_tmp_mat,cn51_sf_30d_tmp_mat);
    r_30d(i) = r_tmp(1,2);
    r_tmp = corrcoef(usgs_sf_tmp_mat,cn51_sf_15d_tmp_mat);
    r_15d(i) = r_tmp(1,2);
    r_tmp = corrcoef(usgs_sf_tmp_mat,cn51_sf_10d_tmp_mat);
    r_10d(i) = r_tmp(1,2);

    mean_usgs_sf(i) = nanmean(usgs_sf_tmp_mat(:));
    mean_cn51_30d(i) = nanmean(cn51_sf_30d_tmp_mat(:));
    mean_cn51_15d(i) = nanmean(cn51_sf_15d_tmp_mat(:));
    mean_cn51_10d(i) = nanmean(cn51_sf_10d_tmp_mat(:));
 
    diff_30d = cn51_sf_30d_tmp_mat(:) - usgs_sf_tmp_mat(:);
    abs_bias_30d(i) = nanmean(abs(diff_30d)); 
    bias_30d(i) = nanmean(diff_30d);
    diff_15d = cn51_sf_15d_tmp_mat(:) - usgs_sf_tmp_mat(:);
    abs_bias_15d(i) = nanmean(abs(diff_15d));
    bias_15d(i) = nanmean(diff_15d);
    diff_10d = cn51_sf_10d_tmp_mat(:) - usgs_sf_tmp_mat(:);
    abs_bias_10d(i) = nanmean(abs(diff_10d));
    bias_10d(i) = nanmean(diff_10d);

    rmse_30d(i) = sqrt(nanmean(diff_30d.^2));
    rmse_15d(i) = sqrt(nanmean(diff_15d.^2));
    rmse_10d(i) = sqrt(nanmean(diff_10d.^2));

    r2_30d = 1 - (sum(diff_30d.^2)./sum((usgs_sf_tmp_mat - mean_usgs_sf(i)).^2));
    r2_15d = 1 - (sum(diff_15d.^2)./sum((usgs_sf_tmp_mat - mean_usgs_sf(i)).^2));
    r2_10d = 1 - (sum(diff_10d.^2)./sum((usgs_sf_tmp_mat - mean_usgs_sf(i)).^2));

end

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

latlim = [25 50];
lonlim = [-125 -60];

%cmap_mat = cmocean('solar',12);
%cmap_mat = colormap(hot(13));
cmap_mat = cbrewer2('YlOrRd',12);
%cmap_mat(end,:) = [];
cmap_mat = flipud(cmap_mat);
cmap = cmap_mat
caxis_min = 0;
caxis_max = 1;

cmap_incr = 1/12;
cmap_index = ceil(r_30d/cmap_incr);
cmap_index(cmap_index<0) = 1;

figure
ax = axesm('miller','MapLatLimit',latlim,'MapLonLimit',lonlim,... 
    'FLatLimit',latlim,'FLonLimit',lonlim);
framem off; 
gridm on;
setm(ax,'MLineLocation',10,'PLineLocation',10,'MeridianLabel','off',...      
    'ParallelLabel','on','MLabelLocation',10,'PLabelLocation',10,...       
'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [1 1 1])
%geoshow(states)
mean_vec = [];
for i = 1:length(r_30d)
    if (isnan(r_30d(i))==0) & (cmap_index(i)~=0)
        scatterm(plot_station_lat(i),plot_station_lon(i),30,r_30d(i),"filled",'MarkerFaceColor',cmap(cmap_index(i),:),...
        'MarkerEdgeColor','k')
        mean_vec = [mean_vec r_30d(i)];
    end
end
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(['30d Catch vs USGS; R mean: ' num2str(nanmean(mean_vec)) ' (N = ' num2str(length(mean_vec)) ')'])
h = colorbar;
set(h,'Location','southoutside','Orientation','horizontal','FontSize',14,'XTick',[0 0.25 0.5 0.75 1]);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,'/discover/nobackup/jkolassa/CN51_plots/global/Catch_streamflow_R_30d','png')

cmap_index = ceil(r_15d/cmap_incr);
cmap_index(cmap_index<0) = 1;

figure
ax = axesm('miller','MapLatLimit',latlim,'MapLonLimit',lonlim,... 
    'FLatLimit',latlim,'FLonLimit',lonlim);
framem off;
gridm on;
setm(ax,'MLineLocation',10,'PLineLocation',10,'MeridianLabel','off',...      
    'ParallelLabel','on','MLabelLocation',10,'PLabelLocation',10,...       
'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [1 1 1])
%geoshow(states)
mean_vec = [];
for i = 1:length(r_15d)
    if (isnan(r_15d(i))==0) & (cmap_index(i)~=0)
        scatterm(plot_station_lat(i),plot_station_lon(i),30,r_15d(i),"filled",'MarkerFaceColor',cmap(cmap_index(i),:),...
        'MarkerEdgeColor','k')
        mean_vec = [mean_vec r_15d(i)];
    end
end
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(['15d Catch vs USGS; R mean: ' num2str(nanmean(mean_vec)) ' (N = ' num2str(length(mean_vec)) ')'])
h = colorbar;
set(h,'Location','southoutside','Orientation','horizontal','FontSize',14,'XTick',[0 0.25 0.5 0.75 1]);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,'/discover/nobackup/jkolassa/CN51_plots/global/Catch_streamflow_R_15d','png')

cmap_index = ceil(r_10d/cmap_incr);
cmap_index(cmap_index<0) = 1;

figure
ax = axesm('miller','MapLatLimit',latlim,'MapLonLimit',lonlim,... 
    'FLatLimit',latlim,'FLonLimit',lonlim);
framem off;
gridm on;
setm(ax,'MLineLocation',10,'PLineLocation',10,'MeridianLabel','off',...      
    'ParallelLabel','on','MLabelLocation',10,'PLabelLocation',10,...       
'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [1 1 1])
%geoshow(states)
mean_vec = [];
for i = 1:length(r_10d)
    if (isnan(r_10d(i))==0) & (cmap_index(i)~=0)
        scatterm(plot_station_lat(i),plot_station_lon(i),30,r_10d(i),"filled",'MarkerFaceColor',cmap(cmap_index(i),:),...
        'MarkerEdgeColor','k')
        mean_vec = [mean_vec r_10d(i)];
    end
end
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(['10d Catch vs USGS; R mean: ' num2str(nanmean(mean_vec)) ' (N = ' num2str(length(mean_vec)) ')'])
h = colorbar;
set(h,'Location','southoutside','Orientation','horizontal','FontSize',14,'XTick',[0 0.25 0.5 0.75 1]);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,'/discover/nobackup/jkolassa/CN51_plots/global/Catch_streamflow_R_10d','png')

%cmap_mat = cmocean('solar',12);
%cmap_mat = colormap(hot(13));
cmap_mat = cbrewer2('BrBG',12);
%cmap_mat(end,:) = [];
cmap = cmap_mat
s1 = quantile(bias_30d,0.05);
s2 = quantile(bias_30d,0.95);
caxis_min = -1*max(abs(s1),abs(s2));
caxis_max = max(abs(s1),abs(s2));

cmap_incr = (caxis_max-caxis_min)/12;
cmap_index = ceil((bias_30d-caxis_min)/cmap_incr);
cmap_index(cmap_index<0) = 1;
cmap_index(cmap_index>length(cmap_mat)) = length(cmap_mat);

figure
ax = axesm('miller','MapLatLimit',latlim,'MapLonLimit',lonlim,... 
    'FLatLimit',latlim,'FLonLimit',lonlim);
framem off;
gridm on;
setm(ax,'MLineLocation',10,'PLineLocation',10,'MeridianLabel','off',...      
    'ParallelLabel','on','MLabelLocation',10,'PLabelLocation',10,...       
'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [1 1 1])
%geoshow(states)
mean_vec = [];
for i = 1:length(bias_30d)
    if (isnan(bias_30d(i))==0) & (cmap_index(i)~=0)
        scatterm(plot_station_lat(i),plot_station_lon(i),30,bias_30d(i),"filled",'MarkerFaceColor',cmap(cmap_index(i),:),...
        'MarkerEdgeColor','k')
        mean_vec = [mean_vec bias_30d(i)];
    end
end
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(['30d Catch vs USGS; bias mean: ' num2str(nanmean(mean_vec)) ' (N = ' num2str(length(mean_vec)) ')'])
h = colorbar;
set(h,'Location','southoutside','Orientation','horizontal','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,'/discover/nobackup/jkolassa/CN51_plots/global/Catch_streamflow_bias_30d','png')

%figure 
%geoscatter(plot_station_lat,plot_station_lon,30,bias_30d)
%caxis([caxis_min caxis_max]);
%set(gcf,'Colormap',cmap);
%tightmap
%title(['30d Catch vs USGS; bias mean: ' num2str(nanmean(mean_vec)) ' (N = ' num2str(length(mean_vec)) ')'])
%h = colorbar;
%set(h,'Location','southoutside','Orientation','horizontal','FontSize',14);
%set(gcf,'paperunits','centimeters')
%set(gcf,'paperposition',[1 1 17 15])
%saveas(gcf,'/discover/nobackup/jkolassa/CN51_plots/global/test','png')

%cmap_mat = cmocean('solar',12);
%cmap_mat = colormap(hot(13));
cmap_mat = cbrewer2('BrBG',12);
%cmap_mat(end,:) = [];
cmap = cmap_mat
s1 = quantile(bias_15d,0.05);
s2 = quantile(bias_15d,0.95);
caxis_min = -1*max(abs(s1),abs(s2));
caxis_max = max(abs(s1),abs(s2));

cmap_incr = (caxis_max-caxis_min)/12;

cmap_index = ceil((bias_15d-caxis_min)/cmap_incr);
cmap_index(cmap_index<0) = 1;
cmap_index(cmap_index>length(cmap_mat)) = length(cmap_mat);

figure
ax = axesm('miller','MapLatLimit',latlim,'MapLonLimit',lonlim,... 
    'FLatLimit',latlim,'FLonLimit',lonlim);
framem off;
gridm on;
setm(ax,'MLineLocation',10,'PLineLocation',10,'MeridianLabel','off',...      
    'ParallelLabel','on','MLabelLocation',10,'PLabelLocation',10,...       
'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [1 1 1])
%geoshow(states)
mean_vec = [];
for i = 1:length(bias_15d)
    if (isnan(bias_15d(i))==0) & (cmap_index(i)~=0)
        scatterm(plot_station_lat(i),plot_station_lon(i),30,bias_15d(i),"filled",'MarkerFaceColor',cmap(cmap_index(i),:),...
        'MarkerEdgeColor','k')
        mean_vec = [mean_vec bias_15d(i)];
    end
end 
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(['15d Catch vs USGS; bias mean: ' num2str(nanmean(mean_vec)) ' (N = ' num2str(length(mean_vec)) ')'])
h = colorbar;
set(h,'Location','southoutside','Orientation','horizontal','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,'/discover/nobackup/jkolassa/CN51_plots/global/Catch_streamflow_bias_15d','png')

%cmap_mat = cmocean('solar',12);
%cmap_mat = colormap(hot(13));
cmap_mat = cbrewer2('BrBG',12);
%cmap_mat(end,:) = [];
cmap = cmap_mat
s1 = quantile(bias_10d,0.05);
s2 = quantile(bias_10d,0.95);
caxis_min = -1*max(abs(s1),abs(s2));
caxis_max = max(abs(s1),abs(s2));

cmap_incr = (caxis_max-caxis_min)/12;

cmap_index = ceil((bias_10d-caxis_min)/cmap_incr);
cmap_index(cmap_index<0) = 1;
cmap_index(cmap_index>length(cmap_mat)) = length(cmap_mat);

figure
ax = axesm('miller','MapLatLimit',latlim,'MapLonLimit',lonlim,... 
    'FLatLimit',latlim,'FLonLimit',lonlim);
framem off;
gridm on;
setm(ax,'MLineLocation',10,'PLineLocation',10,'MeridianLabel','off',...      
    'ParallelLabel','on','MLabelLocation',10,'PLabelLocation',10,...       
'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [1 1 1])
%geoshow(states)
mean_vec = [];
for i = 1:length(bias_10d)
    if (isnan(bias_10d(i))==0) & (cmap_index(i)~=0)
        scatterm(plot_station_lat(i),plot_station_lon(i),30,bias_10d(i),"filled",'MarkerFaceColor',cmap(cmap_index(i),:),...
        'MarkerEdgeColor','k')
        mean_vec = [mean_vec bias_10d(i)];
    end
end 
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(['10d Catch vs USGS; bias mean: ' num2str(nanmean(mean_vec)) ' (N = ' num2str(length(mean_vec)) ')'])
h = colorbar;
set(h,'Location','southoutside','Orientation','horizontal','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,'/discover/nobackup/jkolassa/CN51_plots/global/Catch_streamflow_bias_10d','png')
