%% k_stats by E Peterson 8-Mar-2021
% DOI 
sourcepath='D:\ARC317Jan2021\';
% WorldPop, The spatial distribution of population in 2020, 
% in Global High Resolution Population Denominators Project - Funded by The Bill and Melinda Gates Foundation (OPP1134076). 
% U.o.S.D.o.G.a.G. WorldPop (www.worldpop.org - School of Geography and Environmental Science
% University of Louisville; Departement de Geographie
% Universite de Namur) and Center for International Earth Science Information Network (CIESIN)
% Columbia University (2018). , Editor. 2019.
geographyfile='geog.csv';
geography=strcat(sourcepath,geographyfile);
% mapfile='D:\WB\from_arc3\RESULTS_map.csv'; % n=14,711 records (>=3yrs) with GIS intersects
version='1987-2020';
stats=strcat(sourcepath,'RESULTSstats_',version,'.csv');
rates=strcat(sourcepath,'RESULTSrates_',version,'.csv');
out=strcat(sourcepath,'RESULTSout_',version,'.csv');
%% population density included in geography file
opts = delimitedTextImportOptions("NumVariables", 17);
% Specify range and delimiter
opts.DataLines = [2, Inf];
opts.Delimiter = ",";
% Specify column names and types
opts.VariableNames = ["WMOstn", "WBANo", "GRIDCODE", "Kottek2006", "Peel2007", "Peel", "Becketal", "Beck", "consensus", "latitude", "longitude", "elev", "Placename", "CC", "ST", "province", "gpw_v4_pop"];
opts.VariableTypes = ["double", "double", "double", "categorical", "double", "categorical", "double", "categorical", "categorical", "double", "double", "double", "string", "categorical", "string", "categorical", "double"];
opts = setvaropts(opts, [13, 15], "WhitespaceRule", "preserve");
opts = setvaropts(opts, [4, 6, 8, 9, 13, 14, 15, 16], "EmptyFieldRule", "auto");
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";
% Import the data
geog = readtable("D:\ARC317Jan2021\geog.csv", opts);

% Clear temporary variables
clear opts
%% stats
opts = delimitedTextImportOptions("NumVariables", 42);
% Specify range and delimiter
opts.DataLines = [2, Inf];
opts.Delimiter = ",";
% Specify column names and types
opts.VariableNames = ["WMOstn", "WBANo", "Bracket3", "adjR2_H", "ob_H", "pVal_H", "DuWa_H", "pVDWH", "Bracket9", "adjR2_C", "obC", "pVal_C", "DuWaC", "pVDWC", "Bracket15", "adjR2_E", "ob_E", "pVal_E", "DuWaE", "pVDWE", "Bracket21", "adjR2_D", "obD", "pValD", "DuWaD", "pVDWD", "Bracket27", "adjR2_W", "obW", "pValW", "DuWaW", "pVDWW", "Bracket33", "adjR2_Y", "obY", "pValY", "DuWaY", "pVDWY", "Bracket39", "base", "stop", "Placename"];
opts.VariableTypes = ["double", "double", "categorical", "double", "double", "double", "double", "double", "categorical", "double", "double", "double", "double", "double", "categorical", "double", "double", "double", "double", "double", "categorical", "double", "double", "double", "double", "double", "categorical", "double", "double", "double", "double", "double", "categorical", "double", "double", "double", "double", "double", "categorical", "double", "double", "string"];
opts = setvaropts(opts, 42, "WhitespaceRule", "preserve");
opts = setvaropts(opts, [3, 9, 15, 21, 27, 33, 39, 42], "EmptyFieldRule", "auto");
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";
% Import the data
stats = readtable(stats, opts);
% Clear temporary variables
clear opts
%% rates
opts = delimitedTextImportOptions("NumVariables", 42);
% Specify range and delimiter
opts.DataLines = [2, Inf];
opts.Delimiter = ",";
% Specify column names and types
opts.VariableNames = ["WMOstn", "WBANo", "Bracket3", "HdD18", "std_H", "slope_H", "SEE_H", "HDD18", "Bracket9", "CdD24", "stdC", "slope_C", "SEEC", "CDD24", "Bracket15", "EnthDys", "stdE", "slope_E", "SEEE", "EnthD", "Bracket21", "DXd", "std_D", "slope_D", "SEE_D", "DXD", "Bracket27", "Wsd", "std_W", "slope_W", "SEE_W", "WsD", "Bracket33", "Ysd", "std_Y", "slope_Y", "SEE_Y", "YsD", "Bracket39", "base", "stop", "Placename"];
opts.VariableTypes = ["double", "double", "categorical", "double", "double", "double", "double", "double", "categorical", "double", "double", "double", "double", "double", "categorical", "double", "double", "double", "double", "double", "categorical", "double", "double", "double", "double", "double", "categorical", "double", "double", "double", "double", "double", "categorical", "double", "double", "double", "double", "double", "categorical", "double", "double", "string"];
opts = setvaropts(opts, 42, "WhitespaceRule", "preserve");
opts = setvaropts(opts, [3, 9, 15, 21, 27, 33, 39, 42], "EmptyFieldRule", "auto");
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";
% Import the data
rates = readtable(rates, opts);
% Clear temporary variables
clear opts
%% out
opts = delimitedTextImportOptions("NumVariables", 43);
% Specify range and delimiter
opts.DataLines = [2, Inf];
opts.Delimiter = ",";
% Specify column names and types
opts.VariableNames = ["WMOstn", "WBANo", "Acz", "HVAC", "meanPa", "WBref", "hBASE", "dpy", "Bracket9", "avHDD", "avCDD", "avChD", "Bracket13", "lowshed", "minshed", "aveshed", "maxshed", "bestshed", "Bracket19", "WsD", "YsD", "Bracket22", "yrs", "dys", "Bracket25", "Icy", "Hea", "Hum", "Pas", "Vtn", "EvC", "DX", "Bracket33", "Cd", "Fd", "latitude", "longitude", "elev", "base", "stop", "CC", "ST", "Placename"];
opts.VariableTypes = ["double", "double", "double", "categorical", "double", "double", "double", "double", "categorical", "double", "double", "double", "categorical", "double", "double", "double", "double", "double", "categorical", "double", "double", "categorical", "double", "double", "categorical", "double", "double", "double", "double", "double", "double", "double", "categorical", "double", "double", "double", "double", "double", "double", "double", "categorical", "string", "string"];
opts = setvaropts(opts, [42, 43], "WhitespaceRule", "preserve");
opts = setvaropts(opts, [4, 9, 13, 19, 22, 25, 33, 41, 42, 43], "EmptyFieldRule", "auto");
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";
% Import the data
out = readtable(out, opts);
% Clear temporary variables
clear opts
%% check that WMO and WBAN are aligned in all 4 files
if isequal(stats(:,1:2), rates(:,1:2)) && isequal(stats(:,1:2), out(:,1:2)) && isequal(stats(:,1:2), geog(:,1:2)) 
%% start of working block
% geog [17] popdensity per square km /// province [16] Koppen-Geiger [3:9]
geovec=[3:9 16 17];
for geodex=geovec
    eval(strcat(geog.Properties.VariableNames{geodex},'=geog{:,',num2str(geodex),'};'));
end
% stats[3   9     15     21     27     33     ]
statsvect=[1 2 4:8 10:14  16:20  22:26  28:32  34:38];
for statsdex=statsvect % 
    eval(strcat(stats.Properties.VariableNames{statsdex},'=stats{:,',num2str(statsdex),'};'));
end
% rates  [3   9     15     21     27     33     ]
ratesvect=[4:8 10:14  16:20  22:26  28:32  34:38];
for ratesdex=ratesvect % 
    eval(strcat(rates.Properties.VariableNames{ratesdex},'=rates{:,',num2str(ratesdex),'};'));
end
% out  [    9,    13,    19,    22,    25,    33]
outvect=[3:8 10:12  14:18  20 21  23 24  26:32  34:43];
for outdex=outvect % 
    eval(strcat(out.Properties.VariableNames{outdex},'=out{:,',num2str(outdex),'};'));
end
%% now find delta_H and delta_C
delta_H(1:length(WMOstn))=NaN; delta_H=delta_H';
delta_C(1:length(WMOstn))=NaN; delta_C=delta_C';
pos_H(1:length(WMOstn))=NaN; pos_H=pos_H';
q1_H(1:length(WMOstn))=NaN; q1_H=q1_H';
q1_C(1:length(WMOstn))=NaN; q1_C=q1_C';
pos_C(1:length(WMOstn))=NaN; pos_C=pos_C';
q2_H(1:length(WMOstn))=NaN; q2_H=q2_H';
q2_C(1:length(WMOstn))=NaN; q2_C=q2_C';
neg_H(1:length(WMOstn))=NaN; neg_H=neg_H';
q3_H(1:length(WMOstn))=NaN; q3_H=q3_H';
q3_C(1:length(WMOstn))=NaN; q3_C=q3_C';
neg_C(1:length(WMOstn))=NaN; neg_C=neg_C';
q4_H(1:length(WMOstn))=NaN; q4_H=q4_H';
q4_C(1:length(WMOstn))=NaN; q4_C=q4_C';
for dd=1:length(WMOstn)
    if ~isnan(ob_H(dd)) && ~isnan(pVal_H(dd)) && ~isnan(slope_H(dd)) && ob_H(dd)>=14 && pVal_H(dd)<=0.05 && slope_H(dd)~=0
        delta_H(dd)=slope_H(dd)./HdD18(dd);
    end
end
for dd=1:length(WMOstn)
    if ~isnan(obC(dd)) && ~isnan(pVal_C(dd)) && ~isnan(slope_C(dd)) && obC(dd)>=14 && pVal_C(dd)<=0.05 && slope_C(dd)~=0
        delta_C(dd)=slope_C(dd)./CdD24(dd);
    end
end
%% plotting axis and quadrant without population scaller
eightmax(1:8)=0; % absolutely
eightcount(1:8)=0; % pos_H q1 pos_C q2 neg_H q3 neg_C q4
eightsum(1:8)=0; % pos_H q1 pos_C q2 neg_H q3 neg_C q4
%                    1    2   3    4   5    6   7    8
for dd=1:length(WMOstn)
    if ~isnan(delta_H(dd)) && ~isnan(delta_C(dd)) 
        if delta_C(dd) == 0
            if delta_H(dd) > 0 
                pos_H(dd)=delta_H(dd);
                if ~isinf(pos_H(dd))
                    eightmax(1)=max(eightmax(1),abs(pos_H(dd)));
                    eightsum(1)=eightsum(1)+    abs(pos_H(dd)); % axis 1 / eight 1
                    eightcount(1)=eightcount(1)+1;
                end
            else
                if delta_H(dd) < 0
                    neg_H(dd)=-delta_H(dd);
                    if ~isinf(neg_H(dd))
                        eightmax(5)=max(eightmax(5),abs(neg_H(dd)));
                        eightsum(5)=eightsum(5)+    abs(neg_H(dd)); % axis 3 / eight 5
                        eightcount(5)=eightcount(5)+1;
                    end
                end
            end
        else
            if delta_H(dd) > 0 && delta_C(dd) > 0
                q1_H(dd)=delta_H(dd);
                q1_C(dd)=delta_C(dd);
                if ~isinf(q1_H(dd)) && ~isinf(q1_C(dd))
                    eightmax(2)=max(eightmax(2),sqrt(q1_H(dd)^2+q1_C(dd)^2));
                    eightsum(2)=eightsum(2)+    sqrt(q1_H(dd)^2+q1_C(dd)^2); % quadrant 1 / eight 2\
                    eightcount(2)=eightcount(2)+1;
                end
            else
                if delta_H(dd) < 0 && delta_C(dd) > 0
                    q2_H(dd)=delta_H(dd);
                    q2_C(dd)=delta_C(dd);
                    if ~isinf(q2_H(dd)) && ~isinf(q2_C(dd))
                        eightmax(4)=max(eightmax(4),sqrt(q2_H(dd)^2+q2_C(dd)^2));
                        eightsum(4)=eightsum(4)+    sqrt(q2_H(dd)^2+q2_C(dd)^2); % quadrant 2 / eight 4
                        eightcount(4)=eightcount(4)+1;
                    end
                end
            end
        end
    end
    if ~isnan(delta_H(dd)) && isnan(delta_C(dd))
        if delta_H(dd) > 0 
            pos_H(dd)=delta_H(dd);
            if ~isinf(pos_H(dd))
                eightmax(1)=max(eightmax(1),abs(pos_H(dd)));
                eightsum(1)=eightsum(1)+    abs(pos_H(dd)); % axis 1 / eight 1
                eightcount(1)=eightcount(1)+1;
            end
        else
            if delta_H(dd) < 0
                neg_H(dd)=-delta_H(dd);
                if ~isinf(neg_H(dd))
                    eightmax(5)=max(eightmax(5),abs(neg_H(dd)));
                    eightsum(5)=eightsum(5)+    abs(neg_H(dd)); % axis 3 / eight 5
                    eightcount(5)=eightcount(5)+1;
                end
            end
        end
    end
end
for dd=1:length(WMOstn)
    if ~isnan(delta_H(dd)) && ~isnan(delta_C(dd)) 
        if delta_H(dd) == 0
            if delta_C(dd) > 0 
                pos_C(dd)=delta_C(dd);
                if ~isinf(pos_C(dd))
                    eightmax(3)=max(eightmax(3),abs(pos_C(dd)));
                    eightsum(3)=eightsum(3)+    abs(pos_C(dd)); % axis 2 / eight 3
                    eightcount(3)=eightcount(3)+1;
                end
            else
                if delta_C(dd) < 0
                    neg_C(dd)=-delta_C(dd);
                    if ~isinf(neg_C(dd))
                        eightmax(7)=max(eightmax(7),abs(neg_C(dd)));
                        eightsum(7)=eightsum(7)+    abs(neg_C(dd)); % axis 4 / eight 7
                        eightcount(7)=eightcount(7)+1;
                    end
                end
            end
        else
            if delta_H(dd) < 0 && delta_C(dd) < 0
                q3_H(dd)=delta_H(dd);
                q3_C(dd)=delta_C(dd);
                if  ~isinf(q3_H(dd)) && ~isinf(q3_C(dd))
                    eightmax(6)=max(eightmax(6),sqrt(q3_H(dd)^2+q3_C(dd)^2));
                    eightsum(6)=eightsum(6)+    sqrt(q3_H(dd)^2+q3_C(dd)^2); % quadrant 3 / eight 6
                    eightcount(6)=eightcount(6)+1;
                end
            else
                if delta_H(dd) > 0 && delta_C(dd) < 0
                    q4_H(dd)=delta_H(dd);
                    q4_C(dd)=delta_C(dd);
                    if  ~isinf(q4_H(dd)) && ~isinf(q4_C(dd))
                        eightmax(8)=max(eightmax(8),sqrt(q4_H(dd)^2+q4_C(dd)^2));
                        eightsum(8)=eightsum(8)+    sqrt(q4_H(dd)^2+q4_C(dd)^2); % quadrant 4 / eight 8
                        eightcount(8)=eightcount(8)+1;
                    end
                end
            end
        end
    end
    if ~isnan(delta_C(dd)) && isnan(delta_H(dd))
        if delta_C(dd) > 0 
            pos_C(dd)=delta_C(dd);
            if ~isinf(pos_C(dd))
                eightmax(3)=max(eightmax(3),abs(pos_C(dd)));
                eightsum(3)=eightsum(3)+    abs(pos_C(dd)); % axis 2 / eight 3
                eightcount(3)=eightcount(3)+1;
            end
        else
            if delta_C(dd) < 0
                neg_C(dd)=-delta_C(dd);
                if ~isinf(neg_C(dd))
                    eightmax(7)=max(eightmax(7),abs(neg_C(dd)));
                    eightsum(7)=eightsum(7)+    abs(neg_C(dd)); % axis 4 / eight 7
                    eightcount(7)=eightcount(7)+1;
                end
            end
        end
    end
end
sp=100*eightsum/sum(eightsum);
for step=1:8; spp{step}=strcat(num2str(round(sp(step))),'%'); end
s_str = {{'\DeltaH>0',spp{1}},{'H&C',spp{2}},{'\DeltaC>0',spp{3}},{'H~>C',spp{4}},{'\DeltaH<0',spp{5}},{'',spp{6}},{'\DeltaC<0',spp{7}},{'C~>H',spp{8}}};

eightx(1) = eightmax(1);   eighty(1) = 0;
eightx(3) = 0;             eighty(3) = eightmax(3); 
eightx(5) =-eightmax(5);   eighty(5) = 0;
eightx(7) = 0;             eighty(7) =-eightmax(7); 
eightx(2) = eightmax(1)/2; eighty(2) = eightmax(3)/2;
eightx(4) =-eightmax(5)/2; eighty(4) = eightmax(3)/2;
eightx(6) =-eightmax(5)/2; eighty(6) =-eightmax(7)/2;
eightx(8) = eightmax(1)/2; eighty(8) =-eightmax(7)/2;

%% how many dp loc ChD NaNs and yrs<3
dpcount=0;
loccount=0;
ChDcount=0;
yrscount=0;
for gd=1:length(WMOstn)
    if dpy(gd)<1  
        dpcount=dpcount+1;
    else
        if isnan(latitude(gd))
            loccount=loccount+1;
        else
            if isnan(avChD(gd))
                ChDcount=ChDcount+1;
            else
                if yrs(gd)<3
                    yrscount=yrscount+1;
                end
            end
        end
    end
end
throughout=sum([dpcount loccount ChDcount yrscount]);
validsample=length(WMOstn)-throughout; % 14,711
sprintf('n=%g valid sites after throwing out those without dew point data, enthalpy, geography, or less than 3 years',validsample)
%% pop density gpw_v4_pop with latitude, longitude, 
nanmean(gpw_v4_pop)
sumdens=0; popcount=0;
for gd=1:length(WMOstn)
    if dpy(gd)>=1 && ~isnan(latitude(gd)) && ~isnan(avChD(gd)) && yrs(gd)>=3 && ~isnan(gpw_v4_pop(gd)) && gpw_v4_pop(gd)~=0 
        sumdens=sumdens+gpw_v4_pop(gd);
        popcount=popcount+1;
    end
end
qapd=sumdens/popcount; % qualified average pop density
sprintf('n=%g validated sites intersected by non-zero population data, with %5.1f average density per km^2',popcount, qapd)
% multiply delta_H and delta_C each by (sitepopdens/qapd)
p_delta_H=delta_H.*(gpw_v4_pop/qapd);
p_delta_C=delta_C.*(gpw_v4_pop/qapd);
%
%% plotting axis and quadrant WITH population weighting
p_eightmax(1:8)=0; % absolutely
p_eightcount(1:8)=0;
p_eightsum(1:8)=0; % pos_H q1 pos_C q2 neg_H q3 neg_C q4
%                     1    2   3    4   5    6   7    8
for dd=1:length(WMOstn)
    if ~isnan(p_delta_H(dd)) && ~isnan(p_delta_C(dd)) 
        if p_delta_C(dd) == 0
            if p_delta_H(dd) > 0 
                p_pos_H(dd)=p_delta_H(dd);
                if ~isinf(p_pos_H(dd))
                    p_eightmax(1)=max(p_eightmax(1),abs(p_pos_H(dd)));
                    p_eightsum(1)=p_eightsum(1)+    abs(p_pos_H(dd)); % axis 1 / eight 1
                    p_eightcount(1)=p_eightcount(1)+1;
                end
            else
                if p_delta_H(dd) < 0
                    p_neg_H(dd)=-p_delta_H(dd);
                    if ~isinf(p_neg_H(dd))
                        p_eightmax(5)=max(p_eightmax(5),abs(p_neg_H(dd)));
                        p_eightsum(5)=p_eightsum(5)+    abs(p_neg_H(dd)); % axis 3 / eight 5
                        p_eightcount(5)=p_eightcount(5)+1;
                    end
                end
            end
        else
            if p_delta_H(dd) > 0 && p_delta_C(dd) > 0
                p_q1_H(dd)=p_delta_H(dd);
                p_q1_C(dd)=p_delta_C(dd);
                if ~isinf(p_q1_H(dd)) && ~isinf(p_q1_C(dd))
                    p_eightmax(2)=max(p_eightmax(2),sqrt(p_q1_H(dd)^2+p_q1_C(dd)^2));
                    p_eightsum(2)=p_eightsum(2)+    sqrt(p_q1_H(dd)^2+p_q1_C(dd)^2); % quadrant 1 / eight 2
                    p_eightcount(2)=p_eightcount(2)+1;
                end
            else
                if p_delta_H(dd) < 0 && p_delta_C(dd) > 0
                    p_q2_H(dd)=p_delta_H(dd);
                    p_q2_C(dd)=p_delta_C(dd);
                    if ~isinf(p_q2_H(dd)) && ~isinf(p_q2_C(dd))
                        p_eightmax(4)=max(p_eightmax(4),sqrt(p_q2_H(dd)^2+p_q2_C(dd)^2));
                        p_eightsum(4)=p_eightsum(4)+    sqrt(p_q2_H(dd)^2+p_q2_C(dd)^2); % quadrant 2 / eight 4
                        p_eightcount(4)=p_eightcount(4)+1;
                    end
                end
            end
        end
    end
    if ~isnan(p_delta_H(dd)) && isnan(p_delta_C(dd))
        if p_delta_H(dd) > 0 
            p_pos_H(dd)=p_delta_H(dd);
            if ~isinf(p_pos_H(dd))
                p_eightmax(1)=max(p_eightmax(1),abs(p_pos_H(dd)));
                p_eightsum(1)=p_eightsum(1)+    abs(p_pos_H(dd)); % axis 1 / eight 1
                p_eightcount(1)=p_eightcount(1)+1;
            end
        else
            if p_delta_H(dd) < 0
                p_neg_H(dd)=-p_delta_H(dd);
                if ~isinf(p_neg_H(dd))
                    p_eightmax(5)=max(p_eightmax(5),abs(p_neg_H(dd)));
                    p_eightsum(5)=p_eightsum(5)+    abs(p_neg_H(dd)); % axis 3 / eight 5
                    p_eightcount(5)=p_eightcount(5)+1;
                end
            end
        end
    end
end
for dd=1:length(WMOstn)
    if ~isnan(p_delta_H(dd)) && ~isnan(p_delta_C(dd)) 
        if p_delta_H(dd) == 0
            if p_delta_C(dd) > 0 
                p_pos_C(dd)=p_delta_C(dd);
                if ~isinf(p_pos_C(dd))
                    p_eightmax(3)=max(p_eightmax(3),abs(p_pos_C(dd)));
                    p_eightsum(3)=p_eightsum(3)+    abs(p_pos_C(dd)); % axis 2 / eight 3
                    p_eightcount(3)=p_eightcount(3)+1;
                end
            else
                if p_delta_C(dd) < 0
                    p_neg_C(dd)=-p_delta_C(dd);
                    if ~isinf(p_neg_C(dd))
                        p_eightmax(4)=max(p_eightmax(4),abs(p_neg_C(dd)));
                        p_eightsum(4)=p_eightsum(4)+    abs(p_neg_C(dd)); % axis 4 / eight 7
                        p_eightcount(4)=p_eightcount(4)+1;
                    end
                end
            end
        else
            if p_delta_H(dd) < 0 && p_delta_C(dd) < 0
                p_q3_H(dd)=p_delta_H(dd);
                p_q3_C(dd)=p_delta_C(dd);
                if ~isinf(p_q3_H(dd)) && ~isinf(p_q3_C(dd))
                    p_eightmax(6)=max(p_eightmax(6),sqrt(p_q3_H(dd)^2+p_q3_C(dd)^2));
                    p_eightsum(6)=p_eightsum(6)+    sqrt(p_q3_H(dd)^2+p_q3_C(dd)^2); % quadrant 3 / eight 6
                    p_eightcount(6)=p_eightcount(6)+1;
                end
            else
                if p_delta_H(dd) > 0 && p_delta_C(dd) < 0
                    p_q4_H(dd)=p_delta_H(dd);
                    p_q4_C(dd)=p_delta_C(dd);
                    if ~isinf(p_q4_H(dd)) && ~isinf(p_q4_C(dd))
                        p_eightmax(8)=max(p_eightmax(8),sqrt(p_q4_H(dd)^2+p_q4_C(dd)^2));
                        p_eightsum(8)=p_eightsum(8)+    sqrt(p_q4_H(dd)^2+p_q4_C(dd)^2); % quadrant 4 / eight 8
                        p_eightcount(8)=p_eightcount(8)+1;
                    end
                end
            end
        end
    end
    if ~isnan(p_delta_C(dd)) && isnan(p_delta_H(dd))
        if p_delta_C(dd) > 0 
            p_pos_C(dd)=p_delta_C(dd);
            if ~isinf(p_pos_C(dd))
                p_eightmax(3)=max(p_eightmax(3),abs(p_pos_C(dd)));
                p_eightsum(3)=p_eightsum(3)+    abs(p_pos_C(dd)); % axis 2 / eight 3
                p_eightcount(3)=p_eightcount(3)+1;
            end
        else
            if p_delta_C(dd) < 0
                p_neg_C(dd)=-p_delta_C(dd);
                if ~isinf(p_neg_C(dd))
                    p_eightmax(7)=max(p_eightmax(7),abs(p_neg_C(dd)));
                    p_eightsum(7)=p_eightsum(7)+    abs(p_neg_C(dd)); % axis 4 / eight 7
                    p_eightcount(7)=p_eightcount(7)+1;
                end
            end
        end
    end
end
pp=100*p_eightsum/sum(p_eightsum);
for step=1:8; ppp{step}=strcat(num2str(round(pp(step))),'%'); end
p_str = {{'\DeltaH>0',ppp{1}},{'H&C',ppp{2}},{'\DeltaC>0',ppp{3}},{'H~>C',ppp{4}},{'\DeltaH<0',ppp{5}},{'',ppp{6}},{'\DeltaC<0',ppp{7}},{'C~>H',ppp{8}}};

p_eightx(1) = p_eightmax(1);   p_eighty(1) = 0;
p_eightx(3) = 0;               p_eighty(3) = p_eightmax(3); 
p_eightx(5) =-p_eightmax(5);   p_eighty(5) = 0;
p_eightx(7) = 0;               p_eighty(7) =-p_eightmax(7); 
p_eightx(2) = p_eightmax(1)/2; p_eighty(2) = p_eightmax(3)/2;
p_eightx(4) =-p_eightmax(5)/2; p_eighty(4) = p_eightmax(3)/2;
p_eightx(6) =-p_eightmax(5)/2; p_eighty(6) =-p_eightmax(7)/2;
p_eightx(8) = p_eightmax(1)/2; p_eighty(8) =-p_eightmax(7)/2;

%% graphics
pointsize =5; 
cyan        = [0.2  0.8  0.8];
malachite   = [0.04 0.85 0.32];
orange      = [1    0.5  0];
red         = [1    0.2  0.2];
maroon      = [0.78 0.13 0.28];
green       = [0    0.6  0.3];
navy        = [0    0    1];
blue        = [0    0.5  1];
%
figure('Position', [10 10 1000 600]);
subplot(1,2,1)
% scatter(delta_H,delta_C, '.')
scatter(pos_H,0*pos_H,pointsize,cyan,'filled')
hold on
scatter(q1_H,q1_C,pointsize,malachite,'filled')
scatter(0*pos_C,pos_C,pointsize,orange,'filled')
scatter(q2_H,q2_C,pointsize,red,'filled')
scatter(-neg_H,0*neg_H,pointsize,maroon,'filled')
scatter(q3_H,q3_C,pointsize,green,'filled')
scatter(0*neg_C,-neg_C,pointsize,navy,'filled')
scatter(q4_H,q4_C,pointsize,blue,'filled')
ytickformat('percentage')
xtickformat('percentage')
text(eightx,eighty,s_str)
title('annual HDD and CDD % change per annum')
%%

subplot(1,2,2)
scatter(p_pos_H,0*p_pos_H,pointsize,cyan,'filled')
hold on
scatter(p_q1_H,p_q1_C,pointsize,malachite,'filled')
scatter(0*p_pos_C,p_pos_C,pointsize,orange,'filled')
scatter(p_q2_H,p_q2_C,pointsize,red,'filled')
scatter(-p_neg_H,0*p_neg_H,pointsize,maroon,'filled')
scatter(p_q3_H,p_q3_C,pointsize,green,'filled')
scatter(0*p_neg_C,-p_neg_C,pointsize,navy,'filled')
scatter(p_q4_H,p_q4_C,pointsize,blue,'filled')
ytickformat('percentage')
xtickformat('percentage')
text(p_eightx,p_eighty,p_str)
title('pop weighted H&C DD % change per annum')
%
print('Prediliction.svg','-dsvg');
%% end of working block
else
    sprintf('these files do NOT have identical WMO and WBAN columns')
    stats
    rates
    out
end
%%
% of of script