Elaboradar  0.1
 Tutto Classi Namespace File Funzioni Variabili Tipi enumerati (enum) Gruppi
classifier.h
Vai alla documentazione di questo file.
1 
7 #include <radarelab/volume.h>
8 
9 #include<string>
10 #include<iostream>
11 #include<fstream>
12 #include<sstream>
13 
14 // TODO: prima o poi arriviamo a far senza di questi define
15 #define NUM_AZ_X_PPI 400
16 
17 namespace radarelab {
18 namespace volume {
19 
20 /*
21  * ==================================================================
22  * Enum: EchoClass
23  * Description: classes of radar echoes
24  * ==================================================================
25  */
34 enum EchoClass {
35  GC_AP, // ground clutter or anomalous propagation
36  BS, // biological scatterers
37  DS, // dry aggregated snow
38  WS, // wet snow
39  CR, // crystals of various orientation
40  GR, // graupel
41  BD, // big drops
42  RA, // light and moderate rain
43  HR, // heavy rain
44  RH, // mixture of rain and hail
45  NC // not classified
46 };
47 
48 inline std::ostream& operator<<(std::ostream& oss, EchoClass& ECl)
49 {
50  switch(ECl)
51  {
52  case GC_AP:
53  oss<<"GC_AP";
54  break;
55  case BS:
56  oss<<"BS";
57  break;
58  case DS:
59  oss<<"DS";
60  break;
61  case WS:
62  oss<<"WS";
63  break;
64  case CR:
65  oss<<"CR";
66  break;
67  case GR:
68  oss<<"GR";
69  break;
70  case BD:
71  oss<<"BD";
72  break;
73  case RA:
74  oss<<"RA";
75  break;
76  case HR:
77  oss<<"HR";
78  break;
79  case RH:
80  oss<<"RH";
81  break;
82  case NC:
83  oss<<"NC";
84  break;
85  default:
86  oss<<"unknown echo type "<<ECl;
87  }
88  return oss;
89 }
90 
91 /*
92  * ==================================================================
93  * Class: PROB
94  * Description: Matrix of probability (trapezoidal)
95  * ==================================================================
96  */
100 class PROB : public Matrix2D<double>
101 {
102 public:
106  PROB(double z, double zdr, double rhohv, double lkdp, double sdz, double sdphidp, double vrad);
107 
108 private:
112  double f_1(double Z) {return -0.05+2.5e-3*Z+7.5e-4*Z*Z;}
116  double f_2(double Z) {return 0.68-4.81e-2*Z+2.92e-3*Z*Z;}
120  double f_3(double Z) {return 1.42+6.67e-2*Z+4.85e-4*Z*Z;}
124  double g_1(double Z) {return -44.0+0.8*Z;}
128  double g_2(double Z) {return -22.0+0.5*Z;}
132  double trap(double x1, double x2, double x3, double x4, double val);
136  Matrix2D<double> prob_class(EchoClass classe,double z, double zdr, double rhohv, double lkdp, double sdz, double sdphidp, double vrad);
137 };
138 
139 /*
140  * ==================================================================
141  * Class: MLpoints
142  * Description: store ML points in a azimuth height matrix
143  * ==================================================================
144  */
149 class MLpoints : public Matrix2D<unsigned>
150 {
151 public:
152  double Hmin;
153  double Hmax;
154  unsigned count;
155 
156  MLpoints(double minHeight,double maxHeight,unsigned az_count,unsigned height_count)
157  : Matrix2D<unsigned>(Matrix2D::Constant(height_count,az_count,0)),
158  Hmin(minHeight),Hmax(maxHeight),count(0) {}
159 
160  double azimuth_deg(unsigned az_idx){return (double)az_idx*360./(double)this->cols();}
161  double azimuth_rad(unsigned az_idx){return (double)az_idx*2.*M_PI/(double)this->cols();}
162  unsigned rad2idx(double rad){return (unsigned)(rad*(double)this->cols()/(2.*M_PI));}
163  unsigned deg2idx(double deg){return (unsigned)(deg*(double)this->cols()/360.);}
164 
165  double height(unsigned h_idx){return Hmin+(double)h_idx*(Hmax-Hmin)/(double)this->rows();}
166  unsigned h_idx(double height){return (unsigned)((height-Hmin)*(double)this->rows()/(Hmax-Hmin));}
167 
168  void box_top_bottom(double box_width_deg, double bot_th, double top_th, std::vector<double>& ML_b, std::vector<double>& ML_t);
169 };
170 
171 
172 /*
173  * ==================================================================
174  * Class: CONF
175  * Description: Confidence vector
176  * ==================================================================
177  */
181 class CONF : public Matrix2D<double>
182 {
183 public:
184  CONF()
185  {
186  this->resize(6,1);
187  *this<<1.,1.,1.,1.,1.,1.;
188  }
189 
190  CONF(double phidp, double rhohv, double snr,
191  double gradphitheta, double gradphiphi, double gradZtheta, double gradZphi, double gradZdrtheta,double gradZdrphi,
192  double omega=0.9, double alpha=0.) // omega e alpha li predefinisco in attesa di implementare metodi specifici
193  {
194  double phidpZ=250.;
195  double phidpZdr=250.;
196  double delphidpt=10.;
197  double delrhohv1=0.2;
198  double delrhohv2=0.1;
199  double delZdrt=std::pow(10.,0.005);
200  double snrZ=1.; // 10^0.;
201  double snrKdp=1.;
202  double snrZdr=std::pow(10.,0.05);
203  double snrrhohv=std::pow(10.,0.05);
204 
205  double delZdr,csi,chi;
206  if(rhohv<0.8)
207  {
208  delZdr=0.;
209  csi=1.;
210  chi=0.;
211  }
212  else
213  {
214  chi=std::exp(-0.0000137*omega*omega*(gradphitheta*gradphitheta+gradphiphi*gradphiphi));
215  delZdr=0.02*omega*omega*(gradZtheta*gradZdrtheta+gradZphi*gradZdrphi);
216  chi=(1.-rhohv)/delrhohv1;
217  chi=chi*chi;
218  }
219  double delphi=0.02*omega*omega*(gradphitheta*gradZtheta + gradphiphi*gradZphi); // ho messo gradZh == gradZ
220  this->resize(6,1);
221 
222  *this<<std::exp(-0.69*( phidp*phidp/(phidpZ*phidpZ) + snrZ*snrZ/(snr*snr) + alpha*alpha/(50.*50.))),
223  std::exp(-0.69*(phidp*phidp/(phidpZdr*phidpZdr) + delZdr*delZdr/(delZdrt*delZdrt)
224  + (1.-rhohv)*(1.-rhohv)/(delrhohv1*delrhohv1) + snrZdr*snrZdr/(snr*snr) + alpha*alpha/(50.*50.))),
225  std::exp(-0.69*((1.-chi)*(1.-chi)/(delrhohv2*delrhohv2) + (1.-rhohv)*(1.-rhohv)/(delrhohv1*delrhohv1) + snrrhohv*snrrhohv/(snr*snr))),
226  std::exp(-0.69*(delphi*delphi/(delphidpt*delphidpt) + (1.-rhohv)*(1.-rhohv)/(delrhohv1*delrhohv1) + snrKdp*snrKdp/(snr*snr))),
227  std::exp(-0.69*snrZ*snrZ/(snr*snr)), // definire snrSDZ perfavore, io per adesso ci metto snrZ in analogia con snrSDphi che è snrKdp
228  std::exp(-0.69*snrKdp*snrKdp/(snr*snr));
229 
230  }
231 
232 };
233 
234 
235 /*
236  * ==================================================================
237  * Class: HCA_Park
238  * Description: compute aggregation values A_i and return class of echo
239  * ==================================================================
240  */
244 class HCA_Park
245 {
246 public:
247  /*================== MEMBERS ====================*/
251  double z,zdr,rhohv,lkdp,sdz,sdphidp,vrad;
252  double phidp,snr,gradphitheta,gradphiphi,gradZtheta,gradZphi,gradZdrtheta,gradZdrphi;
256  Eigen::VectorXd Ai;
257 
258  /*================== METHODS ====================*/
266  HCA_Park(double Z, double ZDR, double RHOHV, double LKDP, double SDZ, double SDPHIDP, double VRAD,
267  double PHIDP, double SNR, double GPHITH, double GPHIPHI, double GZTH, double GZPHI, double GZDRTH, double GZDRPHI);
271  inline bool non_meteo_echo()
272  {
273  unsigned idx;
274  Ai.maxCoeff(&idx);
275  if(idx==0||idx==1) return true;
276  else return false;
277  }
281  inline bool meteo_echo()
282  {
283  unsigned idx;
284  Ai.maxCoeff(&idx);
285  if(idx==0||idx==1) return false;
286  else return true;
287  }
291  EchoClass echo(double minimum=0.)
292  {
293  unsigned idx;
294  Ai.maxCoeff(&idx);
295  EchoClass hca=NC;
296  if(Ai(idx)>minimum) hca=(EchoClass)idx;
297  return hca;
298  }
302  void clearAi(){ Ai=Eigen::VectorXd::Zero(Ai.size()); }
303 };
304 
305 /*
306  * ==================================================================
307  * Class: MeltingLayer
308  * Description: compute melting layer boundaries from polarimetry
309  * ==================================================================
310  */
317 {
318 public:
319  std::vector<double> top;
320  std::vector<double> bot;
321 
322  Volume<double> vol_z_0_5km;
323  Volume<double> vol_zdr_1km;
324  Volume<double> vol_rhohv_1km;
325 
326  MeltingLayer(Volume<double>& vol_z,Volume<double>& vol_zdr,Volume<double>& vol_rhohv,
327  std::vector< std::vector< std::vector< HCA_Park> > >& HCA);
328 
329  void seek4mlfile(time_t now, MLpoints&);
330  void fill_empty_azimuths();
331 
332 };
333 
334 /*
335  * ==================================================================
336  * Class: Classifier
337  * Description: compute hydrometeor classification Park et al. (2009)
338  * ==================================================================
339  */
350 {
351 public:
352  /*================== MEMBERS ====================*/
360  std::vector< std::vector< std::vector<HCA_Park> > > vol_Ai;
364  std::string pathname;
449 
450 
451  /*================== METHODS ====================*/
459  classifier(const std::string& file);
467  void correct_for_snr();
476  void compute_lkdp();
481  void HCA_Park_2009();
490  void class_designation(unsigned win_rg=1, unsigned win_az=1);
494  void print_ppi_class(int elev=-1);
495 };
496 
497 } // namespace volume
498 } // namespace radarelab
void compute_lkdp()
Initialize vol_lkdp 10log10 of the moving average slope of phidp along beam path. ...
Definition: classifier.cpp:221
Volume< double > vol_lkdp_6km
Definition: classifier.h:416
void correct_for_snr()
correct rhohv and zdr for noise (Schuur et al. 2003)
Definition: classifier.cpp:315
Volume< EchoClass > vol_hca
Definition: classifier.h:356
Volume< double > vol_phidp
Definition: classifier.h:380
Volume< double > vol_snr
Definition: classifier.h:388
MeltingLayer Melting Layer Detection Algorithm MLDA Giangrande et al. 2008.
Definition: classifier.h:316
Definisce le principali strutture che contengono i dati.
void class_designation(unsigned win_rg=1, unsigned win_az=1)
Designate class echo Find the maximum of aggregation values.
Definition: classifier.cpp:440
double Hmax
height min [km]
Definition: classifier.h:153
unsigned count
height max [km]
Definition: classifier.h:154
Volume< double > vol_phidp_2km
Definition: classifier.h:404
EchoClass
List classes of radar echoes Classes defined in Park et al. (2009) Overload &lt;&lt; operator to print clas...
Definition: classifier.h:34
Volume< double > vol_rhohv_2km
Definition: classifier.h:400
Volume< double > vol_sdphidp
Definition: classifier.h:424
Base for all matrices we use, since we rely on row-major data.
Definition: matrix.h:36
PROB(double z, double zdr, double rhohv, double lkdp, double sdz, double sdphidp, double vrad)
Definition: classifier.cpp:32
Volume< double > vol_grad_phi_phi
Definition: classifier.h:444
MLpoints Melting Layer Points matrix AzH.
Definition: classifier.h:149
Volume< double > vol_zdr_2km
Definition: classifier.h:396
classifier(const std::string &file)
Constructor from odim file.
Definition: classifier.cpp:171
void correct_for_attenuation()
correct Z and Zdr for path attenuation
Definition: classifier.cpp:302
Given radar variables compute matrix of probability.
Definition: classifier.h:100
void melting_layer_classification(MeltingLayer &ML)
Check consistency respect to Melting Layer height.
Definition: classifier.cpp:363
Volume< double > vol_phidp_6km
Definition: classifier.h:408
EchoClass echo(double minimum=0.)
Definition: classifier.h:291
Volume< double > vol_grad_zdr_theta
Definition: classifier.h:440
compute hydrometeor classification
Definition: classifier.h:349
Volume< double > vol_grad_z_phi
Definition: classifier.h:428
Volume< double > vol_sdz
Definition: classifier.h:420
std::vector< std::vector< std::vector< HCA_Park > > > vol_Ai
Definition: classifier.h:360
void print_ppi_class(int elev=-1)
print PPI of EchoClass
Definition: classifier.cpp:600
compute confidence vector of radar variables
Definition: classifier.h:181
Volume< double > vol_grad_zdr_phi
Definition: classifier.h:436
Volume< double > vol_grad_phi_theta
Definition: classifier.h:448
MLpoints(double minHeight, double maxHeight, unsigned az_count, unsigned height_count)
counter
Definition: classifier.h:156
Volume< double > vol_zdr
Definition: classifier.h:372
Volume< double > vol_z_1km
Definition: classifier.h:392
void compute_derived_volumes()
Initialize derived input data.
Definition: classifier.cpp:334
Volume< double > vol_lkdp_2km
Definition: classifier.h:412
Volume< double > vol_grad_z_theta
Definition: classifier.h:432
Volume< double > vol_vrad
Definition: classifier.h:384
Volume< double > vol_rhohv
Definition: classifier.h:376
void HCA_Park_2009()
Compute Echo Classes Park et al. (2009) HCA algorithm.
Definition: classifier.cpp:539