How to calculate Probability of detection (hit rate) (POD),False alarm ratio (FAR), Accuracy (ACC) between two precipitation products.
10 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Antworten (1)
shankar sharma
am 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 Kommentare
Siehe auch
Kategorien
Mehr zu Surface and Mesh Plots finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!