How to calculate Probability of detection (hit rate) (POD),False alarm ratio (FAR), Accuracy (ACC) between two precipitation products.
10 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Réponses (1)
shankar sharma
le 18 Août 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 commentaires
Voir également
Catégories
En savoir plus sur Surface and Mesh Plots dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!