Elaboradar 0.1
Caricamento in corso...
Ricerca in corso...
Nessun risultato
classifier.h
Vai alla documentazione di questo file.
1
6
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 */
33
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 }
278
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 }
288
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 }
299
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 */
316class MeltingLayer
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 */
348
496
497} // namespace volume
498} // namespace radarelab
Homogeneous volume with a common beam count for all PolarScans.
Definition volume.h:431
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
PROB(double z, double zdr, double rhohv, double lkdp, double sdz, double sdphidp, double vrad)
void compute_derived_volumes()
Initialize derived input data.
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
Definition classifier.h:416
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.
Volume< double > vol_grad_z_theta
Definition classifier.h:432
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.
Volume< double > vol_grad_zdr_theta
Definition classifier.h:440
void correct_for_attenuation()
correct Z and Zdr for path attenuation
classifier(const std::string &file)
Constructor from odim file.
Volume< double > vol_sdphidp
Definition classifier.h:424
void correct_for_snr()
correct rhohv and zdr for noise (Schuur et al. 2003)
Volume< double > vol_grad_phi_theta
Definition classifier.h:448
void print_ppi_class(int elev=-1)
print PPI of EchoClass
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.
EchoClass
List classes of radar echoes Classes defined in Park et al. (2009) Overload << operator to print clas...
Definition classifier.h:34
Namespace per volume dati.
Definition elev_fin.h:12
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.