%Michael Angus, 2021
%Verification based on hit/miss rates. 
%For this we need the observed 90% percentile value at each location. 
%Use IMD for step 1. Add others/change up in future. 

addpath('C:\Users\micha\Documents\Work_HEPPI\Data\IMD\IMD_conv')
load IMD_processed
load IMD_dates
%First and last dates in 6955 to 7289 (see IMD_dates)
%for reasons unknown the 24th of June is missing, as in 31st Oct 2019.
%in the forecast. So need to skip in the obs to line up correctly. 
%also starting one day later, as forecast date not actual date. 
IMD_mapped_frdates=IMD_map_rorder(:,:,[6956:6978,6980:7290]);
%Points from forecast - need to convert to IMD
%fcst_loc_lon_points=[76 80 94 93 97 101 140 145 177];
%fcst_loc_lat_points=[152 207 83 266 236 145 219 170 212];
%Mumbai 18.3N, 73.3E
%Rajasthan %25N, 74E
%Kerala %10N, 76.6E
%Shimla %32.04N, 76.4E
%Delhi %28.44N, 77.12E
%Hyderabad %17.52N, 77.8E
%Patna %26.4N, 84.8E
%Bhub %20.5N, 85.7E
%Megh %25.56N, 91.52E
%just do this manually, then you'll have them 
IMD_loc_lon_points=[29 32 42 41 44 47 75 78 101];
IMD_loc_lat_points=[48 75 15 103 89 45 81 57 77];
loc_opts={'Mumbai','Rajasthan','Kerala','Shimla','Delhi',...
    'Hyderabad','Patna','Bhubaneswar','Meghalaya'};


%% Baseline 1 - forecast 
%Find 90th percentile at all points. 
%For this, we need to first find IMD_map_rorder with zeros removed.
IMD_map_rorder_drynan=IMD_map_rorder;
IMD_map_rorder_drynan(IMD_map_rorder==0)=NaN;
IMD_90=squeeze(prctile(IMD_map_rorder_drynan,90,3));
%For each location, find a hit rate/false alarm based on how many ensemble
%members get the 90% observed threshold correct. 
%IMERG 90th percentile, in case of calculating IMERG
IMERG_drynan=IMERG_to_IMD_frdates;
IMERG_drynan(IMERG_to_IMD_frdates==0)=NaN;
IMERG_90=squeeze(prctile(IMERG_drynan,90,3));

obs_exceed=IMD_mapped_frdates-IMERG_90;
fcst_exceed=fr12km_d1_ensall-IMERG_90;
UQM_exceed=IMD_QM_DGQM_ensall-IMERG_90;
EMOS_exceed=fr12km_EMOS_ensall-IMERG_90;

date_list=1:334;
%%
%preallocate for speed
%number of obs
ob_map=NaN(length(IMD_lat_points),length(IMD_lon_points));
neg_fa=NaN(length(IMD_lat_points),length(IMD_lon_points));
fcst_mem_hits=cell(129,135);
UQM_mem_hits=cell(129,135);
EMOS_mem_hits=cell(129,135);
fcst_mem_Fas=cell(129,135);
UQM_mem_Fas=cell(129,135);
EMOS_mem_Fas=cell(129,135);
for lati = 1:length(IMD_lat_points)
    for lonj = 1:length(IMD_lon_points)
    %question 1 - does the IMD dataset exceed the 90th percentile?
        if nansum(IMD_mapped_frdates(lati,lonj,:))>0
        %find hits. 
        iob=find(squeeze(obs_exceed(lati,lonj,:))>0);
        %find dry dates. 
        idry=find(squeeze(IMERG_to_IMD_frdates(lati,lonj,:))==0);
        %find not dry but not hit.
        iobn=date_list(~ismember(date_list,[iob;idry]));
        %save to map
        ob_map(lati,lonj)=length(iob);
        neg_fa(lati,lonj)=length(iobn);
        %question 2 - does each ensemble member exceed the same threshold?
            for ens=1:23
            ifc{ens}=find(squeeze(fcst_exceed(lati,lonj,ens,:))>=0);
            ufc{ens}=find(squeeze(UQM_exceed(lati,lonj,ens,:))>=0);
            efc{ens}=find(squeeze(EMOS_exceed(lati,lonj,ens,:))>=0);
            hits_fc{ens}=iob(ismember(iob,ifc{ens}));
            FAs_fc{ens}=iobn(ismember(iobn,ifc{ens}));
            hits_UQM{ens}=iob(ismember(iob,ufc{ens}));
            FAs_UQM{ens}=iobn(ismember(iobn,ufc{ens}));           
            hits_EMOS{ens}=iob(ismember(iob,efc{ens}));
            FAs_EMOS{ens}=iobn(ismember(iobn,efc{ens}));
            end
            for obs=1:length(iob)
            fcst_mem_hits{lati,lonj}(obs)=length(find(cat(1,hits_fc{:})==iob(obs)));
            UQM_mem_hits{lati,lonj}(obs)=length(find(cat(1,hits_UQM{:})==iob(obs)));
            EMOS_mem_hits{lati,lonj}(obs)=length(find(cat(1,hits_EMOS{:})==iob(obs)));
            end
            fcst_uniq_alarms=unique(cat(2,FAs_fc{:}));
            UQM_uniq_alarms=unique(cat(2,FAs_UQM{:}));
            EMOS_uniq_alarms=unique(cat(2,FAs_EMOS{:}));
            for ffas=1:length(fcst_uniq_alarms)
            fcst_mem_Fas{lati,lonj}(ffas)=length(find(cat(2,FAs_fc{:})==fcst_uniq_alarms(ffas)));
            end
            for ufas=1:length(UQM_uniq_alarms)
            UQM_mem_Fas{lati,lonj}(ufas)=length(find(cat(2,FAs_UQM{:})==UQM_uniq_alarms(ufas)));
            end
            for efas=1:length(EMOS_uniq_alarms)
            EMOS_mem_Fas{lati,lonj}(efas)=length(find(cat(2,FAs_EMOS{:})==EMOS_uniq_alarms(efas)));
            end
        end
    end
end
%% Create x and y (Hit rate, False alarm rate) at every grid cell for ROC skill score. 
hr_mat_fcst=zeros(129,135,25); %two longer than the ensemble as need 1 and 0
%at each end for ROC graph - no forecasting and forecasting constantly.
hr_mat_fcst(:,:,1)=1;
fa_mat_fcst=zeros(129,135,25);
fa_mat_fcst(:,:,1)=1;
hr_mat_UQM=zeros(129,135,25);
hr_mat_UQM(:,:,1)=1;
fa_mat_UQM=zeros(129,135,25);
fa_mat_UQM(:,:,1)=1;
hr_mat_EMOS=zeros(129,135,25);
hr_mat_EMOS(:,:,1)=1;
fa_mat_EMOS=zeros(129,135,25);
fa_mat_EMOS(:,:,1)=1;

for li=1:129
    for lj=1:135
        for ensk=1:23
           hr_mat_fcst(li,lj,ensk+1)=length(find(fcst_mem_hits{li,lj}(:)>=ensk))/ob_map(li,lj);
           fa_mat_fcst(li,lj,ensk+1)=length(find(fcst_mem_Fas{li,lj}(:)>=ensk))/neg_fa(li,lj);
           hr_mat_UQM(li,lj,ensk+1)=length(find(UQM_mem_hits{li,lj}(:)>=ensk))/ob_map(li,lj);
           fa_mat_UQM(li,lj,ensk+1)=length(find(UQM_mem_Fas{li,lj}(:)>=ensk))/neg_fa(li,lj);
           hr_mat_EMOS(li,lj,ensk+1)=length(find(EMOS_mem_hits{li,lj}(:)>=ensk))/ob_map(li,lj);
           fa_mat_EMOS(li,lj,ensk+1)=length(find(EMOS_mem_Fas{li,lj}(:)>=ensk))/neg_fa(li,lj);
        end
    end
end

%% Turn this into both a map of roc skill score and the ROC graphs for the
%9 key points. 
figure('Renderer', 'painters', 'Position', [40 40 1200 600])
for locs_ROC=1:9
    subplot(3,3,locs_ROC)
    plot(squeeze(fa_mat_fcst(IMD_loc_lat_points(locs_ROC),...
        IMD_loc_lon_points(locs_ROC),:)),...
        squeeze(hr_mat_fcst(IMD_loc_lat_points(locs_ROC),...
        IMD_loc_lon_points(locs_ROC),:)),...
        'linewidth',1.2,'color',[0 66 37 132]/255);
    hold on
    grid on
    plot(squeeze(fa_mat_UQM(IMD_loc_lat_points(locs_ROC),...
        IMD_loc_lon_points(locs_ROC),:)),...
        squeeze(hr_mat_UQM(IMD_loc_lat_points(locs_ROC),...
        IMD_loc_lon_points(locs_ROC),:)),...
        'linewidth',1.2,'color',[242 9 242 132]/255);
    plot(squeeze(fa_mat_EMOS(IMD_loc_lat_points(locs_ROC),...
        IMD_loc_lon_points(locs_ROC),:)),...
        squeeze(hr_mat_EMOS(IMD_loc_lat_points(locs_ROC),...
        IMD_loc_lon_points(locs_ROC),:)),...
        'linewidth',1.2,'color',[0 204 204 132]/255);
    plot(0:.1:1,0:.1:1,'k:')
    ylim([0 1]);
    xlim([0 1]);
    set(gca,'xtick',0:.2:1);
    set(gca,'ytick',0:.2:1);
    axis square
    title(loc_opts(locs_ROC));
    if locs_ROC==1
    axP = get(gca,'Position');
    legend('Forecast','UQM','EMOS',...
        'location','northwestoutside')
    legend boxoff   
    set(gca, 'Position', axP)
    end
end
    % add some general axes
    ht = text(-5.8,1.65,...
        'Hit Rate');
    set(ht,'Rotation',90)
    set(ht,'FontSize',14)
    ht2=text(-2.7,-0.3,'False Alarm Rate');
    set(ht2,'Fontsize',14)
saveas(gcf,'ROCs_IMERG','png');

%% Skill score map, and also find the best performing in each grid cell.
%For skill score we need Area under curve. Use abs(trapz).
%Need a loop, only works on vectors. 
Area_fcst=NaN(129,135);
Area_UQM=NaN(129,135);
Area_EMOS=NaN(129,135);
for lia=1:129
    for ljb=1:135
        if nansum(IMD_mapped_frdates(lia,ljb,:))>0
Area_fcst(lia,ljb)=abs(trapz(squeeze(fa_mat_fcst(lia,ljb,:)),...
    squeeze(hr_mat_fcst(lia,ljb,:))));
Area_UQM(lia,ljb)=abs(trapz(squeeze(fa_mat_UQM(lia,ljb,:)),...
    squeeze(hr_mat_UQM(lia,ljb,:))));
Area_EMOS(lia,ljb)=abs(trapz(squeeze(fa_mat_EMOS(lia,ljb,:)),...
    squeeze(hr_mat_EMOS(lia,ljb,:))));
        end
    end
end

%% Create 4 panel final figure:

ROCSS_UQM=(Area_UQM-Area_fcst)./(1-Area_fcst);
ROCSS_EMOS=(Area_EMOS-Area_fcst)./(1-Area_fcst);

% Which one is best, forecast UQM or EMOS?
[~,max_mat]=max(cat(3,Area_fcst,Area_UQM,Area_EMOS),[],3);
max_mat(isnan(Area_fcst))=NaN;

%combine matricies together for the figure plotting
ROC_score_mat=cat(3,ROCSS_UQM,ROCSS_EMOS,max_mat,IMERG_90);
%%
cmap_diff=cbrewer2('div','PiYG',10);
cmap_rain=cbrewer2('seq','Blues',10);
cmap_choice=cat(2,[0 66 37;242 9 242; 0 204 204]/255);
climits_diff=[-.5 .5];
climits_rain=[0 150];
climits_choice=[1 3];
unit_diff='BS skill score';
unit_rain='Rainfall (mm/day)';
unit_choice='Best Performance';
NAME_ME='CS_IMERG_perform';
%grid spacing on figure. 
lon_spacing=10; 
lat_spacing=10;
%------------------------------------------------%
%if m_grid is giving errors, make sure lat/lon/map are in double. 
lon_choice=[66.5 100];
lat_choice=[6.5 38];
projection='Equidistant';
figure;
axes('position',[-0.05 0.08 .85 .85]);
for i=1:4
subplot(2,2,i)
m_proj(projection,'lat',lat_choice,'lon',lon_choice);
m_pcolor(IMD_lon_points,IMD_lat_points,CSS_score_mat(:,:,i));
m_grid('xtick',ceil(lon_choice(1)):lon_spacing:ceil(lon_choice(2)),...
    'ytick',floor(lat_choice(1)):lat_spacing:ceil(lat_choice(2)),...
    'ticklength',0.01,...
    'linest',':',...
    'color',[100 100 100]/255,...
    'FontName','AvantGarde',...
    'FontSize',5,...
    'box','off');
m_coast('line','color','k','linewidth',.2);
if ismember(i,[1 2])
colormap(gca,cmap_diff);
caxis(climits_diff)
elseif i==3
colormap(gca,cmap_choice);
caxis(climits_choice);
else
colormap(gca,cmap_rain);
caxis(climits_rain);    
end
%control location and settings of the colorbar
if i==1
h1=colorbar;
h1.Ticks=climits_diff(1):(climits_diff(2)-climits_diff(1))/size(cmap_diff,1):climits_diff(2);
h1.Orientation='vertical';
h1.FontName = 'AvantGarde';    
h1.FontSize = 6;
h1.Box = 'off';
h1.Label.String = unit_diff;
h1.Label.FontSize = 10;
h1.Label.FontName = 'Arial';
h1.Position=[.45 .6 .03 .3]; %left,bottom,width,height
end
if i==3
h2=colorbar;
h2.Ticks=0:1:3;
h2.Orientation='vertical';
h2.FontName = 'AvantGarde';    
h2.FontSize = 6;
h2.Box = 'off';
h2.Label.String = unit_choice;
h2.Label.FontSize = 10;
h2.Label.FontName = 'Arial';
h2.Position=[.45 .12 .03 .3]; %left,bottom,width,height
end
if i==4
h2=colorbar;
h2.Ticks=climits_rain(1):(climits_rain(2)-climits_rain(1))/size(cmap_rain,1):climits_rain(2);;
h2.Orientation='vertical';
h2.FontName = 'AvantGarde';    
h2.FontSize = 6;
h2.Box = 'off';
h2.Label.String = unit_rain;
h2.Label.FontSize = 10;
h2.Label.FontName = 'Arial';
h2.Position=[.88 .12 .03 .3]; %left,bottom,width,height
end
end
saveas(gcf,NAME_ME,'png');

%% Brier skill score
%observations
BS_on=obs_exceed;
BS_on(obs_exceed>0)=1;
BS_on(obs_exceed<0)=0;
%forecast fraction
Bs_fcst_fctn=fcst_exceed;
Bs_fcst_fctn(fcst_exceed>0)=1;
Bs_fcst_fctn(fcst_exceed<0)=0;
Bs_fcst_fctn_tot=squeeze(nansum(Bs_fcst_fctn,3))./23;
%UQM forecast fraction
Bs_UQM_fctn=UQM_exceed;
Bs_UQM_fctn(UQM_exceed>0)=1;
Bs_UQM_fctn(UQM_exceed<0)=0;
Bs_UQM_fctn_tot=squeeze(nansum(Bs_UQM_fctn,3))./23;
%EMOS fraction
Bs_EMOS_fctn=EMOS_exceed;
Bs_EMOS_fctn(EMOS_exceed>0)=1;
Bs_EMOS_fctn(EMOS_exceed<0)=0;
Bs_EMOS_fctn_tot=squeeze(nansum(Bs_EMOS_fctn,3))./23;
%find difference at each grid cell, only for wet points. 
BS_fcst=NaN(129,135);
BS_UQM=NaN(129,135);
BS_EMOS=NaN(129,135);
for lati = 1:length(IMD_lat_points)
    for lonj = 1:length(IMD_lon_points)
    %question 1 - does the IMD dataset exceed the 90th percentile?
        if nansum(IMD_mapped_frdates(lati,lonj,:))>0
        loc_pts=find(squeeze(IMERG_to_IMD_frdates(lati,lonj,:))>0);
        BS_fcst(lati,lonj)=nanmean((Bs_fcst_fctn_tot(lati,lonj,loc_pts)-...
            BS_on(lati,lonj,loc_pts)).^2);
        BS_UQM(lati,lonj)=nanmean((Bs_UQM_fctn_tot(lati,lonj,loc_pts)-...
            BS_on(lati,lonj,loc_pts)).^2);
        BS_EMOS(lati,lonj)=nanmean((Bs_EMOS_fctn_tot(lati,lonj,loc_pts)-...
            BS_on(lati,lonj,loc_pts)).^2);
        end
    end
end
% Brier skill scores
BSS_UQM=1-(BS_UQM./BS_fcst);
BSS_EMOS=1-(BS_EMOS./BS_fcst);
% Which one is best, forecast UQM or EMOS?
[~,min_mat]=min(cat(3,BS_fcst,BS_UQM,BS_EMOS),[],3);
min_mat(isnan(IMD_90))=NaN;
BSS_score_mat=cat(3,BSS_UQM,BSS_EMOS,min_mat,IMERG_90);

%% CRPS
IMERG_fcst_CRPS=NaN(length(IMD_lat_points),length(IMD_lon_points));
IMERG_UQM_CRPS=NaN(length(IMD_lat_points),length(IMD_lon_points));
IMERG_EMOS_CRPS=NaN(length(IMD_lat_points),length(IMD_lon_points));

for i=1:length(IMD_lat_points)
    for j=1:length(IMD_lon_points)
        if isfinite(nanmean(IMD_mapped_frdates(i,j,:)))
        loc_pts_CRPS=find(squeeze(IMERG_to_IMD_frdates(i,j,:))>0);
        IMERG_fcst_CRPS(i,j)=crps(squeeze(fr12km_d1_ensall(i,j,:,loc_pts_CRPS))',...
            squeeze(IMERG_to_IMD_frdates(i,j,loc_pts_CRPS)),'ecdf');
        IMERG_UQM_CRPS(i,j)=crps(squeeze(IMD_QM_DGQM_ensall(i,j,:,loc_pts_CRPS))',...
            squeeze(IMERG_to_IMD_frdates(i,j,loc_pts_CRPS)),'ecdf');
        IMERG_EMOS_CRPS(i,j)=crps(squeeze(fr12km_EMOS_ensall(i,j,:,loc_pts_CRPS))',...
            squeeze(IMERG_to_IMD_frdates(i,j,loc_pts_CRPS)),'ecdf');
        end
    end
    i
end

%% CRPS skill scores
CSS_UQM=1-(IMERG_UQM_CRPS./IMERG_fcst_CRPS);
CSS_EMOS=1-(IMERG_EMOS_CRPS./IMERG_fcst_CRPS);
% Which one is best, forecast UQM or EMOS?
[~,CSS_min_mat]=min(cat(3,IMERG_fcst_CRPS,IMERG_UQM_CRPS,IMERG_EMOS_CRPS),[],3);
CSS_min_mat(isnan(IMD_90))=NaN;
CSS_score_mat=cat(3,CSS_UQM,CSS_EMOS,CSS_min_mat,IMERG_90);