How to calculate Probability of detection (hit rate) (POD),False alarm ratio (FAR), Accuracy (ACC) between two precipitation products.

i want to calculate POD, FAR, and ACC, between 387 station daily data sets with GSMAP precipitaion with similar dimension data sets, code and data sets are attached; data also contain some NaN values.

clc
 clear all
 load('pre_product.mat')
 load('pre_gauge.mat')
for m=1:length(pre_gauge);%%(number of staitons) 
    aa=0;
    bb=0;
    cc=0;
    dd=0;
   
    pre_g_hit=[];
    pre_p_hit=[];
    pre_g_false=[];
    pre_p_false=[];
    pre_g_miss=[];
    pre_p_miss=[];
    
    for k=1:365;%%(day)
        if (pre_gauge(k)>=0.1)&&(pre_product(k)>=0.1)
            aa=aa+1;
            pre_g_hit=[pre_g_hit;pre_gauge(k)];
            pre_p_hit=[pre_p_hit;pre_product(k)];
        elseif (pre_gauge(k)<0.1)&&(pre_product(k)>=0.1)
            bb=bb+1;
            pre_g_false=[pre_g_false;pre_gauge(k)];
            pre_p_false=[pre_p_false;pre_product(k)];
        elseif (pre_gauge(k)>=0.1)&&(pre_product(k)<0.1)
            cc=cc+1;
            pre_g_miss=[pre_g_miss;pre_gauge(k)];
            pre_p_miss=[pre_p_miss;pre_product(k)];
        elseif (pre_gauge(k)<0.1)&&(pre_product(k)<0.1)
            dd=dd+1;
        end
    end
    
    r0=corrcoef(pre_product,pre_gauge);
    R(m)=r0(1,2);
    bias(m)=sum(pre_product-pre_gauge)/sum(pre_gauge);
    bias_hit(m)=sum(pre_p_hit-pre_g_hit)/sum(pre_gauge);
    bias_false(m)=sum(pre_p_false-pre_g_false)/sum(pre_gauge);
    bias_miss(m)=sum(pre_p_miss-pre_g_miss)/sum(pre_gauge);
    RMSE(m)=(sum((pre_product-pre_gauge).^2)/length(pre_gauge))^0.5;
    POD(m)=aa/(aa+cc);
    FAR(m)=bb/(aa+bb);
    ACC(m)=(aa+dd)/(aa+bb+cc+dd);
end

[Mohammed ALsiddig]

Hello Shanka

am interesting to do the same of your work, could you please explane how the data inside these files

Thanks in advance

[shankar sharma]

can you please contact me at [e-mail]

[Mohammed ALsiddig]

Thanks for your reply, I have sent an email, please check it

Answer 1

[shankar sharma]

Recently, i have solved the question,

here are the codes

    pre_gauge=reshape(obs,387,1095); %% making station data timeseries (387 stations and 1095 days (3year))
    pre_product=reshape(gpm,387,1095); %% making satellite data timeseries (387 stations and 1095 days (3year))
    
    %%loop for each stations (it calculates for each station, a total of 387
    %%staitons)
    for m=1:387;%%(number of staitons)
        
        aa=0;
        bb=0;
        cc=0;
        dd=0;
        pre_g_hit=[];
        pre_p_hit=[];
        pre_g_false=[];
        pre_p_false=[];
        pre_g_miss=[];
        pre_p_miss=[];
        
        
        for k=1:1095;  %%(number of day)
            if (pre_gauge(m,k)>=1)&&(pre_product(m,k)>=1)
                aa=aa+1;
                pre_g_hit=[pre_g_hit;pre_gauge(m,k)];
                pre_p_hit=[pre_p_hit;pre_product(m,k)];
            elseif (pre_gauge(m,k)<1)&&(pre_product(m,k)>=1)
                bb=bb+1;
                pre_g_false=[pre_g_false;pre_gauge(m,k)];
                pre_p_false=[pre_p_false;pre_product(m,k)];
            elseif (pre_gauge(m,k)>=1)&&(pre_product(m,k)<1)
                cc=cc+1;
                pre_g_miss=[pre_g_miss;pre_gauge(m,k)];
                pre_p_miss=[pre_p_miss;pre_product(m,k)];
            elseif (pre_gauge(m,k)<1)&&(pre_product(m,k)<1)
                dd=dd+1;
            end
        end
        %% formula for pod, far, acc
        
        POD(m)=aa/(aa+cc);
        FAR(m)=bb/(aa+bb);
        ACC(m)=(aa+dd)/(aa+bb+cc+dd);
        CSI(m)=aa/(aa+bb+cc);
        FBI(m)=(bb+aa)/(aa+cc);
    
    end
    %% average for whole stations
    
    ALL_POD=AA/(AA+CC)
    ALL_FAR=BB/(AA+BB)
    ALL_ACC=(AA+DD)/(AA+BB+CC+DD)
    ALL_FBI=(AA+BB)/(AA+CC)
    ALL_CSI=AA/(AA+BB+CC)