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

10 次查看(过去 30 天)
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 个评论

请先登录,再进行评论。

回答(1 个)

shankar sharma
shankar sharma 2022-8-18
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)

类别

Help CenterFile Exchange 中查找有关 Surface and Mesh Plots 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by