Elaboradar 0.1
Caricamento in corso...
Ricerca in corso...
Nessun risultato
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
17namespace radarelab {
18namespace volume {
19
20/*
21 * ==================================================================
22 * Enum: EchoClass
23 * Description: classes of radar echoes
24 * ==================================================================
25 */
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
48inline 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 */
100class PROB : public Matrix2D<double>
101{
102public:
106 PROB(double z, double zdr, double rhohv, double lkdp, double sdz, double sdphidp, double vrad);
107
108private:
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 */
149class MLpoints : public Matrix2D<unsigned>
150{
151public:
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 */
181class CONF : public Matrix2D<double>
182{
183public:
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 */
245{
246public:
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{
318public:
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{
351public:
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(float Ht, float Hb);
485 void melting_layer_classification(MeltingLayer& ML, float Ht, float Hb);
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
Homogeneous volume with a common beam count for all PolarScans.
Definition: volume.h:431
compute confidence vector of radar variables
Definition: classifier.h:182
EchoClass echo(double minimum=0.)
Definition: classifier.h:291
MLpoints(double minHeight, double maxHeight, unsigned az_count, unsigned height_count)
counter
Definition: classifier.h:156
unsigned count
height max [km]
Definition: classifier.h:154
double Hmax
height min [km]
Definition: classifier.h:153
MLpoints Melting Layer Points matrix AzH.
Definition: classifier.h:150
MeltingLayer Melting Layer Detection Algorithm MLDA Giangrande et al. 2008.
Definition: classifier.h:317
Given radar variables compute matrix of probability.
Definition: classifier.h:101
void compute_derived_volumes()
Initialize derived input data.
Definition: classifier.cpp:334
Volume< double > vol_rhohv
Definition: classifier.h:376
void class_designation(unsigned win_rg=1, unsigned win_az=1)
Designate class echo Find the maximum of aggregation values.
Definition: classifier.cpp:446
Volume< double > vol_lkdp_6km
Definition: classifier.h:416
Volume< double > vol_vrad
Definition: classifier.h:384
Volume< double > vol_zdr_2km
Definition: classifier.h:396
Volume< double > vol_grad_z_phi
Definition: classifier.h:428
Volume< double > vol_phidp_6km
Definition: classifier.h:408
void compute_lkdp()
Initialize vol_lkdp 10log10 of the moving average slope of phidp along beam path.
Definition: classifier.cpp:221
Volume< double > vol_zdr
Definition: classifier.h:372
Volume< double > vol_grad_z_theta
Definition: classifier.h:432
Volume< double > vol_sdz
Definition: classifier.h:420
std::vector< std::vector< std::vector< HCA_Park > > > vol_Ai
Definition: classifier.h:360
Volume< double > vol_grad_zdr_phi
Definition: classifier.h:436
Volume< double > vol_lkdp_2km
Definition: classifier.h:412
Volume< double > vol_grad_phi_phi
Definition: classifier.h:444
Volume< EchoClass > vol_hca
Definition: classifier.h:356
Volume< double > vol_rhohv_2km
Definition: classifier.h:400
void HCA_Park_2009(float Ht, float Hb)
Compute Echo Classes Park et al. (2009) HCA algorithm.
Definition: classifier.cpp:545
Volume< double > vol_grad_zdr_theta
Definition: classifier.h:440
void correct_for_attenuation()
correct Z and Zdr for path attenuation
Definition: classifier.cpp:302
Volume< double > vol_z_1km
Definition: classifier.h:392
Volume< double > vol_phidp
Definition: classifier.h:380
Volume< double > vol_sdphidp
Definition: classifier.h:424
void correct_for_snr()
correct rhohv and zdr for noise (Schuur et al. 2003)
Definition: classifier.cpp:315
Volume< double > vol_grad_phi_theta
Definition: classifier.h:448
void print_ppi_class(int elev=-1)
print PPI of EchoClass
Definition: classifier.cpp:606
Volume< double > vol_phidp_2km
Definition: classifier.h:404
void melting_layer_classification(MeltingLayer &ML, float Ht, float Hb)
Check consistency respect to Melting Layer height.
Definition: classifier.cpp:363
Volume< double > vol_snr
Definition: classifier.h:388
compute hydrometeor classification
Definition: classifier.h:350
EchoClass
List classes of radar echoes Classes defined in Park et al. (2009) Overload << operator to print clas...
Definition: classifier.h:34
String functions.
Definition: cart.cpp:4
Base for all matrices we use, since we rely on row-major data.
Definition: matrix.h:37
Definisce le principali strutture che contengono i dati.