23using namespace volume;
28 unsigned m_h_idx=matrix.h_idx(scan.
height(rg_idx));
29 unsigned m_az_idx=matrix.deg2idx((
double)az_idx*360./scan.
beam_count);
30 matrix(m_h_idx,m_az_idx)++;
34void MLpoints::box_top_bottom(
double box_width_deg,
double bot_th,
double top_th, std::vector<double>& ML_b, std::vector<double>& ML_t)
36 if(bot_th<0.||bot_th>1.) cout<<
"ERROR bot_th must be 0<%<1 "<<endl;
37 if(top_th<0.||top_th>1.) cout<<
"ERROR top_th must be 0<%<1 "<<endl;
38 if(top_th<bot_th) cout<<
"ERROR top_th must be > than bot_th"<<endl;
39 int width=1+2*std::floor(0.5*box_width_deg*this->cols()/360.);
40 int half=0.5*(width-1);
44 for(
unsigned az=0;az<this->cols();az++)
52 for(
int bm=az-half;bm<width;bm++)
54 if(bm<0) round_bm=this->cols()+bm;
55 else if(bm>=this->cols()) round_bm=bm-this->cols();
57 for(
unsigned h=0;h<this->rows();h++)box_count+=(*
this)(h,round_bm);
59 bottom_lim=bot_th*box_count;
60 top_lim=top_th*box_count;
65 for(
unsigned h=0;h<this->rows();h++)
67 for(
int bm=az-half;bm<width;bm++)
69 if(bm<0) round_bm=this->cols()+bm;
70 else if(bm>=this->cols()) round_bm=bm-this->cols();
72 box_count+=(*this)(h,round_bm);
74 if(ML_b[az]<0 && box_count>bottom_lim)ML_b[az]=this->Hmin+h*(this->
Hmax-this->Hmin)/this->rows();
75 if(ML_t[az]<0 && box_count>top_lim) ML_t[az]=this->Hmin+h*(this->
Hmax-this->Hmin)/this->rows();
81void MeltingLayer::seek4mlfile(time_t now,
MLpoints& mlp)
83 unsigned max_time=16*60;
84 ostringstream filename;
94 file.open(filename.str());
97 cout<<endl<<
"Apro MLfile "<<filename.str()<<endl;
101 if(file.eof())
break;
103 m_h_idx=mlp.h_idx(h);
113void MeltingLayer::fill_empty_azimuths()
120 for(
unsigned i=0;i<bot.size();i++)
122 if(bot[i]!=-99) mean_bot.feed(bot[i]);
123 if(top[i]!=-99) mean_top.feed(top[i]);
133 for(
unsigned i=0;i<bot.size();i++)
135 if(bot[i]==-99) bot[i]=mean_bot.mean;
136 if(top[i]==-99) top[i]=mean_top.mean;
141 vector< vector< vector< HCA_Park> > >& HCA)
142 : vol_z_0_5km(NUM_AZ_X_PPI), vol_zdr_1km(NUM_AZ_X_PPI), vol_rhohv_1km(NUM_AZ_X_PPI)
144 cout<<
"\tInizio melting Layer"<<endl;
150 filter(vol_z,vol_z_0_5km,500.,3.,
false);
151 filter(vol_zdr,dummy,1000.,3.,
false);
152 filter(dummy,vol_zdr_1km,1000.,3.,
false);
153 filter(vol_rhohv,dummy,1000.,3.,
false);
154 filter(dummy,vol_rhohv_1km,1000.,3.,
false);
156 cout<<
"filtrati"<<endl;
167 unsigned curr_rg=0,rg_maxz=0,rg_maxzdr=0;
168 bool confirmed=
false;
174 OUT.open(DATE.str());
176 for(
unsigned el=0;el<vol_rhohv_1km.size();el++)
192 if(zdr(az,rg)>0.2) average.feed(zdr(az,rg));
199 if(rho(az,rg)>=0.9 && rho(az,rg)<=0.95 && HCA[el][az][rg].meteo_echo() && z.
height(rg)>MIN_ML_H && z.
height(rg)<MAX_ML_H)
206 if (z(az,curr_rg) > max_z)
213 if (zdr(az,curr_rg)>max_zdr)
215 max_zdr=zdr(az,curr_rg);
225 if(z(az,rg_maxz)>20. && max_z <47 && max_zdr >min_zdr &&
226 max_zdr < 2.5 && z.
height(rg)>melting_points.Hmin)
235 increment(melting_points,z,az,rg);
237 ML<<el<<
"\t"<<az<<
"\t"<<z.
height(rg)<<endl;
246 cout<<endl<<
"I punti ML trovati sono "<<melting_points.count<<endl;
247 cout<<
"Cerco altri file di ML"<<endl;
248 seek4mlfile(vol_z.
load_info->acq_date, melting_points);
249 cout<<
"Ora i punti ML sono "<<melting_points.count<<endl;
250 melting_points.box_top_bottom(20.,0.2,0.8,bot,top);
251 fill_empty_azimuths();
Classes to compute hydrometeor classification.
PolarScan - structure to describe a polarScan containing a matrix of data and conversion factors.
T compute_mean(unsigned minimum=1)
Compute mean of the distribution of x values.
Generic Class to perform statistical analysis Statistic object could be used as accumulator of data a...
const unsigned beam_count
Number of beam_count used ast each elevations.
Homogeneous volume with a common beam count for all PolarScans.
unsigned count
height max [km]
double Hmax
height min [km]
MLpoints Melting Layer Points matrix AzH.
std::shared_ptr< LoadInfo > load_info
Polar volume information.
double diff_height(unsigned rg_start, unsigned rg_end)
Height difference in kilometers (legacy) between two range gates.
unsigned beam_size
Number of samples in each beam.
double elevation
Nominal elevation of this PolarScan, which may be different from the effective elevation of each sing...
double height(unsigned rg, double beam_half_width=0.0)
Height in kilometers (legacy) at range gate for beam elevation + beam_half_width.
unsigned beam_count
Count of beams in this scan.