How to calculate Probability of detection (hit rate) (POD),False alarm ratio (FAR), Accuracy (ACC) between two precipitation products.
10 views (last 30 days)
Show older comments
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
3 Comments
Answers (1)
shankar sharma
on 18 Aug 2022
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)
0 Comments
See Also
Categories
Find more on Surface and Mesh Plots in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!