21 #include <radarlib/radar.hpp>
27 using namespace radarelab;
28 using namespace volume;
30 namespace odim = OdimH5v21;
32 PROB::PROB(
double z,
double zdr,
double rhohv,
double lkdp,
double sdz,
double sdphidp,
double vrad)
35 this->row(0)=prob_class(GC_AP, z, zdr, rhohv, lkdp, sdz, sdphidp, vrad);
36 this->row(1)=prob_class( BS, z, zdr, rhohv, lkdp, sdz, sdphidp, vrad);
37 this->row(2)=prob_class( DS, z, zdr, rhohv, lkdp, sdz, sdphidp, vrad);
38 this->row(3)=prob_class( WS, z, zdr, rhohv, lkdp, sdz, sdphidp, vrad);
39 this->row(4)=prob_class( CR, z, zdr, rhohv, lkdp, sdz, sdphidp, vrad);
40 this->row(5)=prob_class( GR, z, zdr, rhohv, lkdp, sdz, sdphidp, vrad);
41 this->row(6)=prob_class( BD, z, zdr, rhohv, lkdp, sdz, sdphidp, vrad);
42 this->row(7)=prob_class( RA, z, zdr, rhohv, lkdp, sdz, sdphidp, vrad);
43 this->row(8)=prob_class( HR, z, zdr, rhohv, lkdp, sdz, sdphidp, vrad);
44 this->row(9)=prob_class( RH, z, zdr, rhohv, lkdp, sdz, sdphidp, vrad);
47 Matrix2D<double> PROB::prob_class(
EchoClass classe,
double z,
double zdr,
double rhohv,
double lkdp,
double sdz,
double sdphidp,
double vrad)
50 probability=Matrix2D::Constant(1,6,0.);
61 probability<<trap(15.,20.,70.,80.,z),trap(-4.,-2.,1.,2.,zdr),trap(0.5,0.6,0.9,0.95,rhohv),
62 trap(-30.,-25.,10.,20.,lkdp),trap(2.,4.,35.,50.,sdz),trap(20.,30.,50.,60.,sdphidp);
70 probability<<trap(5.,10.,20.,30.,z),trap(0.,2.,10.,12.,zdr),trap(0.3,0.5,0.8,0.83,rhohv),
71 trap(-30.,-25.,10.,10.,lkdp),trap(1.,2.,4.,7.,sdz),trap(8.,10.,40.,60.,sdphidp);
77 probability<<trap(5.,10.,35.,40.,z),trap(-0.3,0.,0.3,0.6,zdr),trap(0.95,0.98,1.,1.01,rhohv),
78 trap(-30.,-25.,10.,20.,lkdp),trap(0.,0.5,3.,6.,sdz),trap(0.,1.,15.,30.,sdphidp);
84 probability<<trap(25.,30.,50.,60.,z),trap(0.5,1.,2.,3.0,zdr),trap(0.88,0.92,0.95,0.985,rhohv),
85 trap(-30.,-25.,10.,20.,lkdp),trap(0.,0.5,3.,6.,sdz),trap(0.,1.,15.,30.,sdphidp);
92 probability<<trap(0.,5.,20.,25.,z),trap(0.1,0.4,3.,3.3,zdr),trap(0.95,0.98,1.,1.01,rhohv),
93 trap(-5.,0.,10.,20.,lkdp),trap(0.,0.5,3.,6.,sdz),trap(0.,1.,15.,30.,sdphidp);
99 probability<<trap(25.,35.,50.,55.,z),trap(-0.3,0,f1,f1+0.3,zdr),trap(0.9,0.97,1.,1.01,rhohv),
100 trap(-30.,-25.,10.,20.,lkdp),trap(0.,0.5,3.,6.,sdz),trap(0.,1.,15.,30.,sdphidp);
106 probability<<trap(20.,25.,45.,50.,z),trap(f2-0.3,f2,f3,f3+1.,zdr),trap(0.92,0.95,1.,1.01,rhohv),
107 trap(g1-1.,g1,g2,g2+1.,lkdp),trap(0.,0.5,3.,6.,sdz),trap(0.,1.,15.,30.,sdphidp);
113 probability<<trap(5.,10.,45.,50.,z),trap(f1-0.3,f1,f2,f2+0.5,zdr),trap(0.95,0.97,1.,1.01,rhohv),
114 trap(g1-1.,g1,g2,g2+1.,lkdp),trap(0.,0.5,3.,6.,sdz),trap(0.,1.,15.,30.,sdphidp);
120 probability<<trap(40.,45.,55.,60.,z),trap(f1-0.3,f1,f2,f2+0.5,zdr),trap(0.92,0.95,1.,1.01,rhohv),
121 trap(g1-1.,g1,g2,g2+1.,lkdp),trap(0.,0.5,3.,6.,sdz),trap(0.,1.,15.,30.,sdphidp);
127 probability<<trap(45.,50.,75.,80.,z),trap(-0.3,0.,f1,f1+0.5,zdr),trap(0.85,0.9,1.,1.01,rhohv),
128 trap(-10.,-4.,g1,g1+1.,lkdp),trap(0.,0.5,3.,6.,sdz),trap(0.,1.,15.,30.,sdphidp);
132 cout<<
"ERROR!!! PROBability calculator does not known the requested echo type "<<classe<<endl;
137 double PROB::trap(
double x1,
double x2,
double x3,
double x4,
double val)
139 if(val<=x3&&val>=x2)
return 1.;
140 else if(val<x2&&val>x1)
return val/(x2-x1)-x1/(x2-x1);
141 else if (val<x4&&val>x3)
return val/(x3-x4)-x4/(x3-x4);
145 HCA_Park::HCA_Park(
double Z,
double ZDR,
double RHOHV,
double LKDP,
double SDZ,
double SDPHIDP,
double VRAD,
146 double PHIDP,
double SNR,
double GPHITH,
double GPHIPHI,
double GZTH,
double GZPHI,
double GZDRTH,
double GZDRPHI)
147 : z(Z),zdr(ZDR),rhohv(RHOHV),lkdp(LKDP),sdz(SDZ),sdphidp(SDPHIDP), vrad(VRAD),phidp(PHIDP),snr(SNR),gradphitheta(GPHITH),
148 gradphiphi(GPHIPHI),gradZtheta(GZTH),gradZphi(GZPHI),gradZdrtheta(GZDRTH),gradZdrphi(GZDRPHI)
150 PROB Pij(z,zdr,rhohv,lkdp,sdz,sdphidp,vrad);
152 CONF Qi(phidp, rhohv, snr, gradphitheta, gradphiphi, gradZtheta, gradZphi, gradZdrtheta, gradZdrphi);
156 Wij << 0.2, 0.4, 1.0, 0.0, 0.6, 0.8,
157 0.4, 0.6, 1.0, 0.0, 0.8, 0.8,
158 1.0, 0.8, 0.6, 0.0, 0.2, 0.2,
159 0.6, 0.8, 1.0, 0.0, 0.2, 0.2,
160 1.0, 0.6, 0.4, 0.5, 0.2, 0.2,
161 0.8, 1.0, 0.4, 0.0, 0.2, 0.2,
162 0.8, 1.0, 0.6, 0.0, 0.2, 0.2,
163 1.0, 0.8, 0.6, 0.0, 0.2, 0.2,
164 1.0, 0.8, 0.6, 1.0, 0.2, 0.2,
165 1.0, 0.8, 0.6, 1.0, 0.2, 0.2;
168 Ai=((Wij.array()*Pij.array()).matrix()*Qi).array()/(Wij*Qi).array();
172 vol_hca(NUM_AZ_X_PPI), pathname(file), vol_z(NUM_AZ_X_PPI),
173 vol_zdr(NUM_AZ_X_PPI), vol_rhohv(NUM_AZ_X_PPI), vol_phidp(NUM_AZ_X_PPI),
174 vol_vrad(NUM_AZ_X_PPI), vol_snr(NUM_AZ_X_PPI), vol_z_1km(NUM_AZ_X_PPI),
175 vol_zdr_2km(NUM_AZ_X_PPI), vol_rhohv_2km(NUM_AZ_X_PPI),
176 vol_phidp_2km(NUM_AZ_X_PPI), vol_phidp_6km(NUM_AZ_X_PPI),
177 vol_lkdp_2km(NUM_AZ_X_PPI), vol_lkdp_6km(NUM_AZ_X_PPI),
178 vol_sdz(NUM_AZ_X_PPI), vol_sdphidp(NUM_AZ_X_PPI),
179 vol_grad_z_phi(NUM_AZ_X_PPI), vol_grad_z_theta(NUM_AZ_X_PPI),
180 vol_grad_zdr_phi(NUM_AZ_X_PPI), vol_grad_zdr_theta(NUM_AZ_X_PPI),
181 vol_grad_phi_phi(NUM_AZ_X_PPI), vol_grad_phi_theta(NUM_AZ_X_PPI)
183 printf(
"il nome del mio file è %s\n",
pathname.c_str());
194 loader_all.request_quantity(odim::PRODUCT_QUANTITY_DBZH,&full_volume_z);
195 loader_all.request_quantity(odim::PRODUCT_QUANTITY_ZDR,&full_volume_zdr);
196 loader_all.request_quantity(odim::PRODUCT_QUANTITY_RHOHV,&full_volume_rhohv);
197 loader_all.request_quantity(odim::PRODUCT_QUANTITY_PHIDP,&full_volume_phidp);
198 loader_all.request_quantity(odim::PRODUCT_QUANTITY_VRAD,&full_volume_vrad);
199 loader_all.request_quantity(odim::PRODUCT_QUANTITY_SNR,&full_volume_snr);
203 auto elev_array = loader_all.get_nominal_elevations();
204 for (
auto i: loader_all.to_load)
205 i.second->normalize_elevations(elev_array);
207 printf(
"Non so se è andato tutto bene, ma almeno sono arrivato in fondo\n");
209 algo::azimuthresample::Closest<double> resampler;
210 resampler.resample_volume(full_volume_z,
vol_z, 1);
211 resampler.resample_volume(full_volume_zdr,
vol_zdr, 1);
212 resampler.resample_volume(full_volume_rhohv,
vol_rhohv, 1);
213 resampler.resample_volume(full_volume_phidp,
vol_phidp, 1);
214 resampler.resample_volume(full_volume_vrad,
vol_vrad, 1);
215 resampler.resample_volume(full_volume_snr,
vol_snr, 1);
224 unsigned half_win6km=12;
227 double phidp1, phidp2;
230 for(
unsigned el=0;el<
vol_phidp.size();el++)
239 for(
unsigned rg=0;rg<
vol_phidp[el].beam_size;rg++)
241 if(rg<half_win6km||rg>(
vol_phidp[el].beam_size-half_win6km-1)){kdp=0.;}
244 phidp1=
vol_phidp[el].get(az,rg-half_win6km);
245 phidp2=
vol_phidp[el].get(az,rg+half_win6km);
246 if(phidp1==undet||phidp1==nodat||phidp2==undet||phidp2==nodat){kdp=0.;}
249 kdp=0.5*(phidp2-phidp1)/6.;
253 if((kdp<-2.)||(kdp>20.))
257 kdp=0.5*(phidp2-phidp1+180.)/6.;
304 for(
unsigned el=0;el<
vol_z.size();el++)
306 for(
unsigned rg=0;rg<
vol_z[el].beam_size;rg++)
317 cout<<
"inizio snr"<<endl;
319 dB2lin(
vol_snr,vol_snr_linear);
321 dB2lin(
vol_zdr,vol_zdr_linear);
323 for(
unsigned el=0;el<
vol_rhohv.size();el++)
327 vol_snr_linear[el].array()*=alpha;
328 vol_zdr_linear[el].array()=vol_snr_linear[el].array()*vol_zdr_linear[el].array()/(vol_snr_linear[el].array()+alpha-vol_zdr_linear[el].array());
330 lin2dB(vol_zdr_linear,
vol_zdr);
331 cout<<
"finito snr"<<endl;
343 printf(
"filtro Z 1 km\n");
345 printf(
"filtro Zdr 2 km\n");
347 printf(
"filtro rhohv 2 km\n");
354 printf(
"calcolo i gradienti\n");
367 for(
unsigned el=0;el<
vol_z.size();el++)
370 cout<<
"El "<<el<<
"\t"<<elev<<endl;
371 for(
unsigned az=0;az<
vol_z.
scan(el).beam_count;az++)
381 vol_Ai[el][az][rg].Ai[DS]=0.;
382 vol_Ai[el][az][rg].Ai[WS]=0.;
383 vol_Ai[el][az][rg].Ai[CR]=0.;
384 vol_Ai[el][az][rg].Ai[GR]=0.;
394 vol_Ai[el][az][rg].Ai[DS]=0.;
395 vol_Ai[el][az][rg].Ai[CR]=0.;
406 vol_Ai[el][az][rg].Ai[CR]=0.;
407 vol_Ai[el][az][rg].Ai[RA]=0.;
408 vol_Ai[el][az][rg].Ai[HR]=0.;
418 vol_Ai[el][az][rg].Ai[RA]=0.;
419 vol_Ai[el][az][rg].Ai[HR]=0.;
427 vol_Ai[el][az][rg].Ai[GC_AP]=0.;
428 vol_Ai[el][az][rg].Ai[BS]=0.;
429 vol_Ai[el][az][rg].Ai[WS]=0.;
430 vol_Ai[el][az][rg].Ai[BD]=0.;
431 vol_Ai[el][az][rg].Ai[RA]=0.;
432 vol_Ai[el][az][rg].Ai[HR]=0.;
444 std::vector< std::vector< std::vector<HCA_Park> > > vol_Ai_filtered;
446 unsigned half_rg=0.5*(win_rg-1);
447 unsigned half_az=0.5*(win_az-1);
448 int pre_rg,post_rg,pre_az,post_az;
451 for(
unsigned el=0;el<
vol_Ai.size();el++)
459 for(
unsigned az=0;az<
vol_Ai[el].size();az++)
460 for(
unsigned rg=0;rg<
vol_Ai[el][az].size();rg++)
463 vol_Ai_filtered[el][az][rg].clearAi();
464 vol_Ai_filtered[el][az][rg].Ai+=
vol_Ai[el][az][rg].Ai;
467 for(
unsigned j=1;j<half_az+1;j++)
472 if(pre_az<0)pre_az+=
vol_Ai[el].size();
473 if(post_az>=
vol_Ai[el].size())post_az-=
vol_Ai[el].size();
475 vol_Ai_filtered[el][az][rg].Ai+=
vol_Ai[el][pre_az][rg].Ai;
476 vol_Ai_filtered[el][az][rg].Ai+=
vol_Ai[el][post_az][rg].Ai;
480 for(
unsigned i=1;i<half_rg+1;i++)
486 vol_Ai_filtered[el][az][rg].Ai+=
vol_Ai[el][az][pre_rg].Ai;
489 if(post_rg<
vol_Ai[el][az].size())
491 vol_Ai_filtered[el][az][rg].Ai+=
vol_Ai[el][az][post_rg].Ai;
495 for(
unsigned j=1;j<half_az+1;j++)
501 if(pre_az<0)pre_az+=
vol_Ai[el].size();
502 if(post_az>=
vol_Ai[el].size())post_az-=
vol_Ai[el].size();
505 vol_Ai_filtered[el][az][rg].Ai+=
vol_Ai[el][pre_az][pre_rg].Ai;
506 vol_Ai_filtered[el][az][rg].Ai+=
vol_Ai[el][pre_az][pre_rg].Ai;
509 if(post_rg<
vol_Ai[el][az].size())
511 vol_Ai_filtered[el][az][rg].Ai+=
vol_Ai[el][post_az][post_rg].Ai;
512 vol_Ai_filtered[el][az][rg].Ai+=
vol_Ai[el][post_az][post_rg].Ai;
518 vol_Ai_filtered[el][az][rg].Ai.array()/=(double)count;
520 vol_hca[el](az,rg)=vol_Ai_filtered[el][az][rg].echo(0.00001);
528 for(
unsigned el=0;el<
vol_z.size();el++)
532 for(
unsigned az=0;az<
vol_z.
scan(el).beam_count;az++)
533 for(
unsigned rg=0;rg<
vol_z.
scan(el).beam_size;rg++)
541 vector< vector< HCA_Park > > SCAN;
542 vector< HCA_Park > BEAM;
543 printf(
"inizio HCA\n");
545 double Z,Zdr,rhohv,lkdp,sdz,sdphidp,vrad;
546 double phidp,snr,gradphitheta,gradphiphi,gradZtheta,gradZphi,gradZdrtheta,gradZdrphi;
547 for(
unsigned el=0;el<
vol_z.size();el++)
549 cout<<
"\tHCA el "<<el<<endl;
551 for(
unsigned az=0;az<
vol_z.
scan(el).beam_count;az++)
554 for(
unsigned rg=0;rg<
vol_z.
scan(el).beam_size;rg++)
572 HCA_Park hca(Z,Zdr,rhohv,lkdp,sdz,sdphidp,vrad,phidp,snr,gradphitheta,gradphiphi,gradZtheta,gradZphi,gradZdrtheta,gradZdrphi);
581 cout<<
"applico ML criteria ad HCA"<<endl;
605 ostringstream filename;
606 filename<<
"class_"<<elev<<
".png";
608 img=
vol_hca[elev].cast<
unsigned short>();
610 write_image(img, filename.str(),
"PNG");
614 for(
unsigned el=0;el<
vol_hca.size();el++)
616 ostringstream filename;
617 filename<<
"class_"<<el<<
".png";
619 img=
vol_hca[el].cast<
unsigned short>();
621 write_image(img, filename.str(),
"PNG");
void compute_lkdp()
Initialize vol_lkdp 10log10 of the moving average slope of phidp along beam path. ...
Volume< double > vol_lkdp_6km
PolarScan< T > & scan(unsigned idx)
Access a polar scan.
void correct_for_snr()
correct rhohv and zdr for noise (Schuur et al. 2003)
Codice per il caricamento di volumi ODIM in radarelab.
Struttura che eredita da Loader e definisce i metodi per accedere ai dati ODIM.
Volume< EchoClass > vol_hca
Volume< double > vol_phidp
MeltingLayer Melting Layer Detection Algorithm MLDA Giangrande et al. 2008.
PolarScan - structure to describe a polarScan containing a matrix of data and conversion factors...
void class_designation(unsigned win_rg=1, unsigned win_az=1)
Designate class echo Find the maximum of aggregation values.
EchoClass
List classes of radar echoes Classes defined in Park et al. (2009) Overload << operator to print clas...
Volume< double > vol_rhohv_2km
Volume< double > vol_sdphidp
Classes to compute hydrometeor classification.
const unsigned beam_count
Number of beam_count used ast each elevations.
PROB(double z, double zdr, double rhohv, double lkdp, double sdz, double sdphidp, double vrad)
std::shared_ptr< LoadInfo > load_info
Polar volume information.
Volume< double > vol_grad_phi_phi
Volume< double > vol_zdr_2km
classifier(const std::string &file)
Constructor from odim file.
void correct_for_attenuation()
correct Z and Zdr for path attenuation
Given radar variables compute matrix of probability.
void melting_layer_classification(MeltingLayer &ML)
Check consistency respect to Melting Layer height.
std::string units
Data units according to ODIM documentation.
Volume< double > vol_phidp_6km
Volume< double > vol_grad_zdr_theta
Volume< double > vol_grad_z_phi
std::string quantity
Odim quantity name.
std::vector< std::vector< std::vector< HCA_Park > > > vol_Ai
void print_ppi_class(int elev=-1)
print PPI of EchoClass
compute confidence vector of radar variables
Volume< double > vol_grad_zdr_phi
void gdal_init_once()
Initialize the GDAL library when called for the first time; does nothing all other times...
Volume< double > vol_grad_phi_theta
Volume< double > vol_z_1km
void compute_derived_volumes()
Initialize derived input data.
Volume< double > vol_lkdp_2km
Volume< double > vol_grad_z_theta
Volume< double > vol_vrad
Volume< double > vol_rhohv
void HCA_Park_2009()
Compute Echo Classes Park et al. (2009) HCA algorithm.