Elaboradar  0.1
 Tutto Classi Namespace File Funzioni Variabili Tipi enumerati (enum) Gruppi
classifier.cpp
1 /*
2  * =====================================================================================
3  *
4  * Filename: classifier.cpp
5  *
6  * Description: implementation of classifier.h methods
7  *
8  * Version: 1.0
9  * Created: 23/06/2014 15:30
10  * Revision:
11  * Compiler: gcc
12  *
13  * Author: Davide Ori
14  * Organization: Dept.Physics University of Bologna
15  *
16  * =====================================================================================
17  */
18 
19 #include "classifier.h"
20 #include <radarelab/odim.h>
21 #include <radarlib/radar.hpp>
22 #include <sstream>
24 #include <radarelab/image.h>
26 
27 using namespace radarelab;
28 using namespace volume;
29 using namespace std;
30 namespace odim = OdimH5v21;
31 
32 PROB::PROB(double z,double zdr,double rhohv, double lkdp, double sdz, double sdphidp, double vrad)
33 {
34  this->resize(10,6);
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);
45 }
46 
47 Matrix2D<double> PROB::prob_class(EchoClass classe,double z, double zdr, double rhohv, double lkdp, double sdz, double sdphidp, double vrad)
48 {
49  Matrix2D<double> probability(1,6);
50  probability=Matrix2D::Constant(1,6,0.);
51  double f1=f_1(z);
52  double f2=f_2(z);
53  double f3=f_3(z);
54  double g1=g_1(z);
55  double g2=g_2(z);
56  switch(classe)
57  {
58  case GC_AP:
59  if(abs(vrad)<2.)
60  {
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);
63  // TODO : il clutter appennino da sdz > 15, quello alpino > 20 cambio le soglie da 2,4,10,15
64  // TODO : Anche con la sdphi voglio essere più cattivo e abbasso la minima da 30,40,50,60
65  }
66  return probability;
67  case BS:
68  if(rhohv<0.97)
69  {
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);
72  }
73  return probability;
74  case DS:
75  if(zdr<2.)
76  {
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);
79  }
80  return probability;
81  case WS:
82  if(z>20.&&zdr>0.)
83  {
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);
86  //TODO alzo le soglie z di ws da 25,30,40,50
87  }
88  return probability;
89  case CR:
90  if(z<40.)
91  {
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);
94  }
95  return probability;
96  case GR:
97  if(z>10.&&z<60.)
98  {
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);
101  }
102  return probability;
103  case BD:
104  if(zdr>(f2-0.3))
105  {
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);
108  }
109  return probability;
110  case RA:
111  if(z<50.)
112  {
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);
115  }
116  return probability;
117  case HR:
118  if(z>30.)
119  {
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);
122  }
123  return probability;
124  case RH:
125  if(z>40.)
126  {
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);
129  }
130  return probability;
131  default:
132  cout<<"ERROR!!! PROBability calculator does not known the requested echo type "<<classe<<endl;
133  return probability; // without it produces compile warnings and runtime error if reached
134  }
135 }
136 
137 double PROB::trap(double x1, double x2, double x3, double x4, double val)
138 {
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);
142  else return 0.; // (val<=x1||val>=x4)
143 }
144 
145 HCA_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)
149 {
150  PROB Pij(z,zdr,rhohv,lkdp,sdz,sdphidp,vrad);
151 // CONF Qi; // TODO: confidence vector calculation could produce shit!!
152  CONF Qi(phidp, rhohv, snr, gradphitheta, gradphiphi, gradZtheta, gradZphi, gradZdrtheta, gradZdrphi); // gradients must be precomputed
153 
154  Matrix2D<double> Wij(10,6);
155 // Z Zdr rhohv lkdp SDZ SDphidp
156  Wij << 0.2, 0.4, 1.0, 0.0, 0.6, 0.8, // GC_AP 3.0
157  0.4, 0.6, 1.0, 0.0, 0.8, 0.8, // BS 3.6
158  1.0, 0.8, 0.6, 0.0, 0.2, 0.2, // DS 2.8
159  0.6, 0.8, 1.0, 0.0, 0.2, 0.2, // WS 2.8
160  1.0, 0.6, 0.4, 0.5, 0.2, 0.2, // CR 2.9
161  0.8, 1.0, 0.4, 0.0, 0.2, 0.2, // GR 2.6
162  0.8, 1.0, 0.6, 0.0, 0.2, 0.2, // BD 2.8
163  1.0, 0.8, 0.6, 0.0, 0.2, 0.2, // RA 2.8
164  1.0, 0.8, 0.6, 1.0, 0.2, 0.2, // HR 3.8
165  1.0, 0.8, 0.6, 1.0, 0.2, 0.2; // RH 3.8
166 // TOT = 8.8 7.6 6.8 2.5 3.0 3.2
167  Ai.resize(10);
168  Ai=((Wij.array()*Pij.array()).matrix()*Qi).array()/(Wij*Qi).array();
169 }
170 
171 classifier::classifier(const string& file):
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)
182 {
183  printf("il nome del mio file è %s\n", pathname.c_str());
184 
185  volume::ODIMLoader loader_all;
186 
187  volume::Scans<double> full_volume_z;
188  volume::Scans<double> full_volume_zdr;
189  volume::Scans<double> full_volume_rhohv;
190  volume::Scans<double> full_volume_phidp;
191  volume::Scans<double> full_volume_vrad;
192  volume::Scans<double> full_volume_snr;
193 
194  loader_all.request_quantity(odim::PRODUCT_QUANTITY_DBZH,&full_volume_z);
195  loader_all.request_quantity(odim::PRODUCT_QUANTITY_ZDR,&full_volume_zdr);
196  loader_all.request_quantity(odim::PRODUCT_QUANTITY_RHOHV,&full_volume_rhohv);
197  loader_all.request_quantity(odim::PRODUCT_QUANTITY_PHIDP,&full_volume_phidp);
198  loader_all.request_quantity(odim::PRODUCT_QUANTITY_VRAD,&full_volume_vrad);
199  loader_all.request_quantity(odim::PRODUCT_QUANTITY_SNR,&full_volume_snr);
200 
201  loader_all.load(pathname);
202 
203  auto elev_array = loader_all.get_nominal_elevations();
204  for (auto i: loader_all.to_load)
205  i.second->normalize_elevations(elev_array);
206 
207  printf("Non so se è andato tutto bene, ma almeno sono arrivato in fondo\n");
208 
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);
216  vol_hca.quantity="CLASS";
217 
218  cout<<vol_z.load_info->acq_date<<endl;
219 }
220 
222 {
224  unsigned half_win6km=12;
225  double kdp;
226  unsigned tries;
227  double phidp1, phidp2;
228  double undet,nodat;
229 
230  for(unsigned el=0;el<vol_phidp.size();el++)
231  {
232  vol_phidp_6km.push_back(vol_phidp[el]);
233  undet=vol_phidp[el].undetect;
234  nodat=vol_phidp[el].nodata;
235 
236  for(unsigned az=0;az<vol_phidp[el].beam_count;az++)
237  {
238  vol_phidp_6km[el].set(az,0,0.);
239  for(unsigned rg=0;rg<vol_phidp[el].beam_size;rg++)
240  {
241  if(rg<half_win6km||rg>(vol_phidp[el].beam_size-half_win6km-1)){kdp=0.;}
242  else
243  {
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.;}
247  else
248  {
249  kdp=0.5*(phidp2-phidp1)/6.;
250  tries=0;
251  while(tries<3)
252  {
253  if((kdp<-2.)||(kdp>20.))
254  {
255  if(kdp<-10.) // vulpiani diceva -20, ma considerava ricetrasmettitori simultanei (360°) e L=7km
256  {
257  kdp=0.5*(phidp2-phidp1+180.)/6.;
258  }
259  else
260  {
261  kdp=0.;
262  }
263  tries++;
264  }
265  else {tries=4;}
266  }
267  }
268  }
269  if(rg){vol_phidp_6km[el].set(az,rg,vol_phidp_6km[el].get(az,rg-1)+2.*kdp*vol_phidp[el].cell_size*0.001);}
270  //vol_lkdp_6km[el](az,rg)=kdp>0.001?10.*log10(kdp):-30;
271  //cout<<"fil 6km rg "<<el<<" "<<rg<<" "<<kdp<<" "<<vol_phidp_6km[el].get(az,rg)<<endl;
272  }
273  }
274  } // finita la ricostruzione di phidp secondo Vulpiani 2012
275  moving_average_slope(vol_phidp_6km,vol_lkdp_6km,6000.,false);
276  moving_average_slope(vol_phidp_6km,vol_lkdp_2km,2000.,false);
277 
278  vol_phidp_6km.quantity="PHIDP";
279  vol_phidp_6km.units="°";
280  vol_lkdp_2km.quantity="LKDP";
281  vol_lkdp_2km.units="°/km";
282  vol_lkdp_6km.quantity="LKDP";
283  vol_lkdp_6km.units="°/km";
284  vol_lkdp_2km*=1000.;
285  vol_lkdp_6km*=1000.;
286  for(unsigned el=0;el<vol_lkdp_6km.size();el++)
287  {
288  vol_lkdp_2km[el].nodata=-9999;
289  vol_lkdp_2km[el].undetect=-30;
290  vol_lkdp_6km[el].nodata=-9999;
291  vol_lkdp_6km[el].undetect=-30;
292  for(unsigned az=0;az<vol_lkdp_6km[el].beam_count;az++)
293  for(unsigned rg=0;rg<vol_lkdp_6km[el].beam_size;rg++)
294  {
295  vol_lkdp_6km[el](az,rg)=vol_lkdp_6km[el](az,rg)>0.001?10*log10(vol_lkdp_6km[el](az,rg)):vol_lkdp_6km[el].undetect;
296  vol_lkdp_2km[el](az,rg)=vol_lkdp_2km[el](az,rg)>0.001?10*log10(vol_lkdp_2km[el](az,rg)):vol_lkdp_2km[el].undetect;
297  }
298  }
299 
300 }
301 
303 {
304  for(unsigned el=0;el<vol_z.size();el++)
305  for(unsigned az=0;az<vol_z[el].beam_count;az++)
306  for(unsigned rg=0;rg<vol_z[el].beam_size;rg++)
307  {
308  if((vol_z[el](az,rg)!=vol_z[el].nodata)&&(vol_z[el](az,rg)!=vol_z[el].undetect))
309  vol_z[el](az,rg)+=0.06*vol_phidp_6km[el](az,rg);
310  if((vol_zdr[el](az,rg)!=vol_zdr[el].nodata)&&(vol_zdr[el](az,rg)!=vol_zdr[el].undetect))
311  vol_zdr[el](az,rg)+=0.01*vol_phidp_6km[el](az,rg);
312  }
313 }
314 
316 {
317  cout<<"inizio snr"<<endl;
318  Volume<double> vol_snr_linear(NUM_AZ_X_PPI);
319  dB2lin(vol_snr,vol_snr_linear);
320  Volume<double> vol_zdr_linear(NUM_AZ_X_PPI);
321  dB2lin(vol_zdr,vol_zdr_linear);
322  double alpha=1.48; // horizontal to vertical noise ratio. 1.48 S-band (Schuur 2003)
323  for(unsigned el=0;el<vol_rhohv.size();el++)
324  {
325  // vol_rhohv[el].array()*=(vol_snr_linear[el].array()+1)/vol_snr_linear[el].array(); // TODO: rhohv è già > 1 in molti casi
326  // questa correzione la aumenta ancora di più
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());
329  }
330  lin2dB(vol_zdr_linear,vol_zdr);
331  cout<<"finito snr"<<endl;
332 }
333 
335 {
336  //correct_phidp(); //TODO: probabilmente inutile se adottiamo il metodo di Vulpiani che stima direttamente la kdp
337  //correct_for_snr(); //TODO: fa danni indicibili questa correzione. Maledetto Schuur 2003 !!!
338 
339  compute_lkdp();
341 
342  // filtro i volumi
343  printf("filtro Z 1 km\n");
344  filter(vol_z,vol_z_1km,1000.,0.,false);
345  printf("filtro Zdr 2 km\n");
346  filter(vol_zdr,vol_zdr_2km,2000.,0.,false); //TODO: test new filter on azimuth to enhance melting characteristics
347  printf("filtro rhohv 2 km\n");
348  filter(vol_rhohv,vol_rhohv_2km,2000.,3.,false);
349 
350  // calcolo le texture
351  textureSD(vol_z,vol_sdz,1000.,0.,false);
352  textureSD(vol_phidp,vol_sdphidp,2000.,0.,true);
353 
354  printf("calcolo i gradienti\n");
355  gradient_azimuth(vol_z_1km, vol_grad_z_phi, false);
356  gradient_elevation(vol_z_1km, vol_grad_z_theta, false);
357  gradient_azimuth(vol_zdr_2km, vol_grad_zdr_phi,false);
358  gradient_elevation(vol_zdr_2km, vol_grad_zdr_theta,false);
359  gradient_azimuth(vol_phidp, vol_grad_phi_phi,true);
360  gradient_elevation(vol_phidp, vol_grad_phi_theta,true);
361 }
362 
364 {
365  double elev;
366  //double Hb,Ht;
367  for(unsigned el=0;el<vol_z.size();el++)
368  {
369  elev=vol_z.scan(el).elevation;
370  cout<<"El "<<el<<"\t"<<elev<<endl;
371  for(unsigned az=0;az<vol_z.scan(el).beam_count;az++)
372  {
373  if ( Ht < 0.0 || Hb < 0.0 ){ //cioè se Ht o Hb sono dati non validi li calcolo
374  Ht=ML.top[az];
375  Hb=ML.bot[az];
376  }
377 
378 
379  bool flag=true;
380  unsigned rg=0;
381  while(flag&&rg<vol_z.scan(el).beam_size)
382  {
383  if(vol_z.scan(el).height(rg,+0.45)<Hb)
384  {
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.;
389  rg++;
390  }
391  else flag=false;
392  }
393  flag=true;
394  while(flag&&rg<vol_z.scan(el).beam_size)
395  {
396  if(vol_z.scan(el).height(rg)<Hb)
397  {
398  vol_Ai[el][az][rg].Ai[DS]=0.;
399  vol_Ai[el][az][rg].Ai[CR]=0.;
400  rg++;
401  //vol_Ai[el][az][rg].Ai[GR]=0.; // TODO: aggiunta a posteriori per vedere l'effetto che fa
402  }
403  else flag=false;
404  }
405  flag=true;
406  while(flag&&rg<vol_z.scan(el).beam_size)
407  {
408  if(vol_z.scan(el).height(rg)<Ht)
409  {
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.;
413  rg++;
414  }
415  else flag=false;
416  }
417  flag=true;
418  while(flag&&rg<vol_z.scan(el).beam_size)
419  {
420  if(vol_z.scan(el).height(rg,-0.45)<Ht)
421  {
422  vol_Ai[el][az][rg].Ai[RA]=0.;
423  vol_Ai[el][az][rg].Ai[HR]=0.;
424  rg++;
425  }
426  else flag=false;
427  }
428  flag=true;
429  while(rg<vol_z.scan(el).beam_size)
430  {
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.;
437  rg++;
438  }
439 
440  }
441 
442  }
443  ;
444 }
445 
446 void classifier::class_designation(unsigned win_rg, unsigned win_az)
447 {
448  if(win_rg||win_az)
449  {
450  std::vector< std::vector< std::vector<HCA_Park> > > vol_Ai_filtered;
451  vol_Ai_filtered=vol_Ai;
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;
455  unsigned count=0;
456  //cout<<half_az<<"\t"<<half_rg<<endl;
457  for(unsigned el=0;el<vol_Ai.size();el++)
458  {
459  vol_hca.push_back(PolarScan<EchoClass>(vol_z.scan(el).beam_count,vol_z.scan(el).beam_size, NC));
460  vol_hca[el].elevation=vol_z[el].elevation;
461  vol_hca[el].gain=1;
462  vol_hca[el].offset=0;
463  vol_hca[el].undetect=NC;
464  vol_hca[el].nodata=255;
465  for(unsigned az=0;az<vol_Ai[el].size();az++)
466  for(unsigned rg=0;rg<vol_Ai[el][az].size();rg++)
467  {
468  //cout<<el<<"\t"<<az<<"\t"<<rg<<endl;
469  vol_Ai_filtered[el][az][rg].clearAi();
470  vol_Ai_filtered[el][az][rg].Ai+=vol_Ai[el][az][rg].Ai;
471  count++;
472  //cout<<"faccio gli az"<<endl;
473  for(unsigned j=1;j<half_az+1;j++)
474  {
475  pre_az=az-j;
476  post_az=az+j;
477  //cout<<pre_az<<"\t"<<post_az<<endl;
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();
480  //cout<<pre_az<<"\t"<<post_az<<endl;
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;
483  count+=2;
484  }
485  //cout<<"faccio gli rg"<<flush;
486  for(unsigned i=1;i<half_rg+1;i++)
487  {
488  pre_rg=rg-i;
489  post_rg=rg+i;
490  if(pre_rg>=0)
491  {
492  vol_Ai_filtered[el][az][rg].Ai+=vol_Ai[el][az][pre_rg].Ai;
493  count++;
494  }
495  if(post_rg<vol_Ai[el][az].size())
496  {
497  vol_Ai_filtered[el][az][rg].Ai+=vol_Ai[el][az][post_rg].Ai;
498  count++;
499  }
500  //cout<<"faccio gli az "<<i<<flush;
501  for(unsigned j=1;j<half_az+1;j++)
502  {
503  pre_az=az-j;
504  post_az=az+j;
505  pre_az=az-j;
506  post_az=az+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();
509  if(pre_rg>=0)
510  {
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;
513  count+=2;
514  }
515  if(post_rg<vol_Ai[el][az].size())
516  {
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;
519  count+=2;
520  }
521  }
522  }
523  //cout<<" fatto"<<endl;
524  vol_Ai_filtered[el][az][rg].Ai.array()/=(double)count;
525  //cout<<"normalizzato"<<endl;
526  vol_hca[el](az,rg)=vol_Ai_filtered[el][az][rg].echo(0.00001);
527  //cout<<"azzero"<<endl;
528  count=0;
529  }
530  }
531  }
532  else
533  {
534  for(unsigned el=0;el<vol_z.size();el++)
535  {
536  vol_hca.push_back(PolarScan<EchoClass>(vol_z.scan(el).beam_count,vol_z.scan(el).beam_size, NC));
537  vol_hca[el].elevation=vol_z[el].elevation;
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++)
540  vol_hca[el](az,rg)=vol_Ai[el][az][rg].echo(0.00001);
541  }
542  }
543 }
544 
545 void classifier::HCA_Park_2009(float Ht, float Hb)
546 {
547  vector< vector< HCA_Park > > SCAN;
548  vector< HCA_Park > BEAM;
549  printf("inizio HCA\n");
550  vol_Ai.resize(vol_z.size());
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++)
554  {
555  cout<<"\tHCA el "<<el<<endl;
556  SCAN.resize(vol_z.scan(el).beam_count);
557  for(unsigned az=0;az<vol_z.scan(el).beam_count;az++)
558  {
559  BEAM.resize(vol_z.scan(el).beam_size);
560  for(unsigned rg=0;rg<vol_z.scan(el).beam_size;rg++)
561  {
562  Z=vol_z_1km.scan(el).get(az,rg);
563  Zdr=vol_zdr_2km.scan(el).get(az,rg);
564  rhohv=vol_rhohv_2km.scan(el).get(az,rg);
565  lkdp=Z>40?vol_lkdp_2km[el].get(az,rg):vol_lkdp_6km[el].get(az,rg);
566  sdz=vol_sdz.scan(el).get(az,rg);
567  sdphidp=vol_sdphidp.scan(el).get(az,rg);
568  vrad=vol_vrad.scan(el).get(az,rg);
569  phidp=vol_phidp[el](az,rg);
570  snr=vol_snr[el](az,rg);
571  gradphitheta=vol_grad_phi_theta[el](az,rg);
572  gradphiphi=vol_grad_phi_phi[el](az,rg);
573  gradZtheta=vol_grad_z_theta[el](az,rg);
574  gradZphi=vol_grad_z_phi[el](az,rg);
575  gradZdrtheta=vol_grad_zdr_theta[el](az,rg);
576  gradZdrphi=vol_grad_zdr_phi[el](az,rg);
577 
578  HCA_Park hca(Z,Zdr,rhohv,lkdp,sdz,sdphidp,vrad,phidp,snr,gradphitheta,gradphiphi,gradZtheta,gradZphi,gradZdrtheta,gradZdrphi);
579  BEAM[rg]=hca;
580  }
581  SCAN[az]=BEAM;
582  }
583  vol_Ai[el]=SCAN;
584  }
585  // Dopo aver calcolato i valori di aggregazione cerco il melting layer
587  cout<<"applico ML criteria ad HCA"<<endl;
588  melting_layer_classification(ML, Ht, Hb);
589  class_designation(1,1);
590 /* unsigned elev=1;
591  unsigned azim=140;
592  cout<<"GC\tBS\tDS\tWS\tCR\tGR\tBD\tRA\tHR\tRH"<<endl;
593  for(unsigned rg=0;rg<vol_Ai[elev][azim].size();rg++)
594  {
595  cout.precision(5);
596  cout<<fixed<<vol_Ai[elev][azim][rg].Ai[GC_AP]<<"\t"<<vol_Ai[elev][azim][rg].Ai[BS]<<"\t"<<
597  vol_Ai[elev][azim][rg].Ai[DS]<<"\t"<<vol_Ai[elev][azim][rg].Ai[WS]<<"\t"<<
598  vol_Ai[elev][azim][rg].Ai[CR]<<"\t"<<vol_Ai[elev][azim][rg].Ai[GR]<<"\t"<<
599  vol_Ai[elev][azim][rg].Ai[BD]<<"\t"<<vol_Ai[elev][azim][rg].Ai[RA]<<"\t"<<
600  vol_Ai[elev][azim][rg].Ai[HR]<<"\t"<<vol_Ai[elev][azim][rg].Ai[RH]<<"\t"<<
601  vol_hca[elev](azim,rg)<<endl;
602  }
603 */
604 }
605 
607 {
608  gdal_init_once();
609  if(elev>=0)
610  {
611  ostringstream filename;
612  filename<<"class_"<<elev<<".png";
614  img=vol_hca[elev].cast<unsigned short>();
615  img*=6553;
616  write_image(img, filename.str(), "PNG");
617  }
618  else
619  {
620  for(unsigned el=0;el<vol_hca.size();el++)
621  {
622  ostringstream filename;
623  filename<<"class_"<<el<<".png";
625  img=vol_hca[el].cast<unsigned short>();
626  img*=6553;
627  write_image(img, filename.str(), "PNG");
628  }
629  }
630 }
void compute_lkdp()
Initialize vol_lkdp 10log10 of the moving average slope of phidp along beam path. ...
Definition: classifier.cpp:221
T offset
Data Offset.
Definition: volume.h:274
Volume< double > vol_lkdp_6km
Definition: classifier.h:416
PolarScan< T > & scan(unsigned idx)
Access a polar scan.
Definition: volume.h:311
void correct_for_snr()
correct rhohv and zdr for noise (Schuur et al. 2003)
Definition: classifier.cpp:315
Codice per il caricamento di volumi ODIM in radarelab.
Struttura che eredita da Loader e definisce i metodi per accedere ai dati ODIM.
Definition: odim.h:22
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
PolarScan - structure to describe a polarScan containing a matrix of data and conversion factors...
Definition: volume.h:112
void class_designation(unsigned win_rg=1, unsigned win_az=1)
Designate class echo Find the maximum of aggregation values.
Definition: classifier.cpp:446
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
Classes to compute hydrometeor classification.
const unsigned beam_count
Number of beam_count used ast each elevations.
Definition: volume.h:432
PROB(double z, double zdr, double rhohv, double lkdp, double sdz, double sdphidp, double vrad)
Definition: classifier.cpp:32
std::shared_ptr< LoadInfo > load_info
Polar volume information.
Definition: volume.h:270
Volume< double > vol_grad_phi_phi
Definition: classifier.h:444
void melting_layer_classification(MeltingLayer &ML, float Ht, float Hb)
Check consistency respect to Melting Layer height.
Definition: classifier.cpp:363
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 HCA_Park_2009(float Ht, float Hb)
Compute Echo Classes Park et al. (2009) HCA algorithm.
Definition: classifier.cpp:545
std::string units
Data units according to ODIM documentation.
Definition: volume.h:268
Volume< double > vol_phidp_6km
Definition: classifier.h:408
Volume< double > vol_grad_zdr_theta
Definition: classifier.h:440
Volume< double > vol_grad_z_phi
Definition: classifier.h:428
Volume< double > vol_sdz
Definition: classifier.h:420
std::string quantity
Odim quantity name.
Definition: volume.h:266
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:606
compute confidence vector of radar variables
Definition: classifier.h:181
Volume< double > vol_grad_zdr_phi
Definition: classifier.h:436
void gdal_init_once()
Initialize the GDAL library when called for the first time; does nothing all other times...
Definition: image.cpp:12
Volume< double > vol_grad_phi_theta
Definition: classifier.h:448
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