21#include <radarlib/radar.hpp>
28using namespace volume;
30namespace odim = OdimH5v21;
32PROB::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);
47Matrix2D<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;
137double 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);
145HCA_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());
196 loader_all.
request_quantity(odim::PRODUCT_QUANTITY_RHOHV,&full_volume_rhohv);
197 loader_all.
request_quantity(odim::PRODUCT_QUANTITY_PHIDP,&full_volume_phidp);
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++)
373 if ( Ht < 0.0 || Hb < 0.0 ){
385 vol_Ai[el][az][rg].Ai[DS]=0.;
386 vol_Ai[el][az][rg].Ai[WS]=0.;
387 vol_Ai[el][az][rg].Ai[CR]=0.;
388 vol_Ai[el][az][rg].Ai[GR]=0.;
398 vol_Ai[el][az][rg].Ai[DS]=0.;
399 vol_Ai[el][az][rg].Ai[CR]=0.;
410 vol_Ai[el][az][rg].Ai[CR]=0.;
411 vol_Ai[el][az][rg].Ai[RA]=0.;
412 vol_Ai[el][az][rg].Ai[HR]=0.;
422 vol_Ai[el][az][rg].Ai[RA]=0.;
423 vol_Ai[el][az][rg].Ai[HR]=0.;
431 vol_Ai[el][az][rg].Ai[GC_AP]=0.;
432 vol_Ai[el][az][rg].Ai[BS]=0.;
433 vol_Ai[el][az][rg].Ai[WS]=0.;
434 vol_Ai[el][az][rg].Ai[BD]=0.;
435 vol_Ai[el][az][rg].Ai[RA]=0.;
436 vol_Ai[el][az][rg].Ai[HR]=0.;
450 std::vector< std::vector< std::vector<HCA_Park> > > vol_Ai_filtered;
452 unsigned half_rg=0.5*(win_rg-1);
453 unsigned half_az=0.5*(win_az-1);
454 int pre_rg,post_rg,pre_az,post_az;
457 for(
unsigned el=0;el<
vol_Ai.size();el++)
465 for(
unsigned az=0;az<
vol_Ai[el].size();az++)
466 for(
unsigned rg=0;rg<
vol_Ai[el][az].size();rg++)
469 vol_Ai_filtered[el][az][rg].clearAi();
470 vol_Ai_filtered[el][az][rg].Ai+=
vol_Ai[el][az][rg].Ai;
473 for(
unsigned j=1;j<half_az+1;j++)
478 if(pre_az<0)pre_az+=
vol_Ai[el].size();
479 if(post_az>=
vol_Ai[el].size())post_az-=
vol_Ai[el].size();
481 vol_Ai_filtered[el][az][rg].Ai+=
vol_Ai[el][pre_az][rg].Ai;
482 vol_Ai_filtered[el][az][rg].Ai+=
vol_Ai[el][post_az][rg].Ai;
486 for(
unsigned i=1;i<half_rg+1;i++)
492 vol_Ai_filtered[el][az][rg].Ai+=
vol_Ai[el][az][pre_rg].Ai;
495 if(post_rg<
vol_Ai[el][az].size())
497 vol_Ai_filtered[el][az][rg].Ai+=
vol_Ai[el][az][post_rg].Ai;
501 for(
unsigned j=1;j<half_az+1;j++)
507 if(pre_az<0)pre_az+=
vol_Ai[el].size();
508 if(post_az>=
vol_Ai[el].size())post_az-=
vol_Ai[el].size();
511 vol_Ai_filtered[el][az][rg].Ai+=
vol_Ai[el][pre_az][pre_rg].Ai;
512 vol_Ai_filtered[el][az][rg].Ai+=
vol_Ai[el][pre_az][pre_rg].Ai;
515 if(post_rg<
vol_Ai[el][az].size())
517 vol_Ai_filtered[el][az][rg].Ai+=
vol_Ai[el][post_az][post_rg].Ai;
518 vol_Ai_filtered[el][az][rg].Ai+=
vol_Ai[el][post_az][post_rg].Ai;
524 vol_Ai_filtered[el][az][rg].Ai.array()/=(double)count;
526 vol_hca[el](az,rg)=vol_Ai_filtered[el][az][rg].echo(0.00001);
534 for(
unsigned el=0;el<
vol_z.size();el++)
538 for(
unsigned az=0;az<
vol_z.
scan(el).beam_count;az++)
539 for(
unsigned rg=0;rg<
vol_z.
scan(el).beam_size;rg++)
547 vector< vector< HCA_Park > > SCAN;
548 vector< HCA_Park > BEAM;
549 printf(
"inizio HCA\n");
551 double Z,Zdr,rhohv,lkdp,sdz,sdphidp,vrad;
552 double phidp,snr,gradphitheta,gradphiphi,gradZtheta,gradZphi,gradZdrtheta,gradZdrphi;
553 for(
unsigned el=0;el<
vol_z.size();el++)
555 cout<<
"\tHCA el "<<el<<endl;
557 for(
unsigned az=0;az<
vol_z.
scan(el).beam_count;az++)
560 for(
unsigned rg=0;rg<
vol_z.
scan(el).beam_size;rg++)
578 HCA_Park hca(Z,Zdr,rhohv,lkdp,sdz,sdphidp,vrad,phidp,snr,gradphitheta,gradphiphi,gradZtheta,gradZphi,gradZdrtheta,gradZdrphi);
587 cout<<
"applico ML criteria ad HCA"<<endl;
611 ostringstream filename;
612 filename<<
"class_"<<elev<<
".png";
614 img=
vol_hca[elev].cast<
unsigned short>();
616 write_image(img, filename.str(),
"PNG");
620 for(
unsigned el=0;el<
vol_hca.size();el++)
622 ostringstream filename;
623 filename<<
"class_"<<el<<
".png";
625 img=
vol_hca[el].cast<
unsigned short>();
627 write_image(img, filename.str(),
"PNG");
Classes to compute hydrometeor classification.
PolarScan - structure to describe a polarScan containing a matrix of data and conversion factors.
const unsigned beam_count
Number of beam_count used ast each elevations.
Homogeneous volume with a common beam count for all PolarScans.
compute confidence vector of radar variables
MeltingLayer Melting Layer Detection Algorithm MLDA Giangrande et al. 2008.
PROB(double z, double zdr, double rhohv, double lkdp, double sdz, double sdphidp, double vrad)
Given radar variables compute matrix of probability.
std::string units
Data units according to ODIM documentation.
std::string quantity
Odim quantity name.
std::shared_ptr< LoadInfo > load_info
Polar volume information.
PolarScan< T > & scan(unsigned idx)
Access a polar scan.
Sequence of PolarScans which can have a different beam count for each elevation.
void compute_derived_volumes()
Initialize derived input data.
Volume< double > vol_rhohv
void class_designation(unsigned win_rg=1, unsigned win_az=1)
Designate class echo Find the maximum of aggregation values.
Volume< double > vol_lkdp_6km
Volume< double > vol_vrad
Volume< double > vol_zdr_2km
Volume< double > vol_grad_z_phi
Volume< double > vol_phidp_6km
void compute_lkdp()
Initialize vol_lkdp 10log10 of the moving average slope of phidp along beam path.
Volume< double > vol_grad_z_theta
std::vector< std::vector< std::vector< HCA_Park > > > vol_Ai
Volume< double > vol_grad_zdr_phi
Volume< double > vol_lkdp_2km
Volume< double > vol_grad_phi_phi
Volume< EchoClass > vol_hca
Volume< double > vol_rhohv_2km
void HCA_Park_2009(float Ht, float Hb)
Compute Echo Classes Park et al. (2009) HCA algorithm.
Volume< double > vol_grad_zdr_theta
void correct_for_attenuation()
correct Z and Zdr for path attenuation
Volume< double > vol_z_1km
Volume< double > vol_phidp
classifier(const std::string &file)
Constructor from odim file.
Volume< double > vol_sdphidp
void correct_for_snr()
correct rhohv and zdr for noise (Schuur et al. 2003)
Volume< double > vol_grad_phi_theta
void print_ppi_class(int elev=-1)
print PPI of EchoClass
void melting_layer_classification(MeltingLayer &ML, float Ht, float Hb)
Check consistency respect to Melting Layer height.
EchoClass
List classes of radar echoes Classes defined in Park et al. (2009) Overload << operator to print clas...
void gdal_init_once()
Initialize the GDAL library when called for the first time; does nothing all other times.
Codice per il caricamento di volumi ODIM in radarelab.
Base for all matrices we use, since we rely on row-major data.
std::map< std::string, Scans< double > * > to_load
Map used to specify quantity to be loaded.
std::vector< double > get_nominal_elevations()
Get list of available elevation.
void request_quantity(const std::string &name, Scans< double > *volume)
Define a request - Fill to_load attribute
void load(const std::string &pathname)
Load method.
Struttura che eredita da Loader e definisce i metodi per accedere ai dati ODIM.