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  Ht=ML.top[az];
374  Hb=ML.bot[az];
375  bool flag=true;
376  unsigned rg=0;
377  while(flag&&rg<vol_z.scan(el).beam_size)
378  {
379  if(vol_z.scan(el).height(rg,+0.45)<Hb)
380  {
381  vol_Ai[el][az][rg].Ai[DS]=0.;
382  vol_Ai[el][az][rg].Ai[WS]=0.;
383  vol_Ai[el][az][rg].Ai[CR]=0.;
384  vol_Ai[el][az][rg].Ai[GR]=0.;
385  rg++;
386  }
387  else flag=false;
388  }
389  flag=true;
390  while(flag&&rg<vol_z.scan(el).beam_size)
391  {
392  if(vol_z.scan(el).height(rg)<Hb)
393  {
394  vol_Ai[el][az][rg].Ai[DS]=0.;
395  vol_Ai[el][az][rg].Ai[CR]=0.;
396  rg++;
397  //vol_Ai[el][az][rg].Ai[GR]=0.; // TODO: aggiunta a posteriori per vedere l'effetto che fa
398  }
399  else flag=false;
400  }
401  flag=true;
402  while(flag&&rg<vol_z.scan(el).beam_size)
403  {
404  if(vol_z.scan(el).height(rg)<Ht)
405  {
406  vol_Ai[el][az][rg].Ai[CR]=0.;
407  vol_Ai[el][az][rg].Ai[RA]=0.;
408  vol_Ai[el][az][rg].Ai[HR]=0.;
409  rg++;
410  }
411  else flag=false;
412  }
413  flag=true;
414  while(flag&&rg<vol_z.scan(el).beam_size)
415  {
416  if(vol_z.scan(el).height(rg,-0.45)<Ht)
417  {
418  vol_Ai[el][az][rg].Ai[RA]=0.;
419  vol_Ai[el][az][rg].Ai[HR]=0.;
420  rg++;
421  }
422  else flag=false;
423  }
424  flag=true;
425  while(rg<vol_z.scan(el).beam_size)
426  {
427  vol_Ai[el][az][rg].Ai[GC_AP]=0.;
428  vol_Ai[el][az][rg].Ai[BS]=0.;
429  vol_Ai[el][az][rg].Ai[WS]=0.;
430  vol_Ai[el][az][rg].Ai[BD]=0.;
431  vol_Ai[el][az][rg].Ai[RA]=0.;
432  vol_Ai[el][az][rg].Ai[HR]=0.;
433  rg++;
434  }
435 
436  }
437  }
438 }
439 
440 void classifier::class_designation(unsigned win_rg, unsigned win_az)
441 {
442  if(win_rg||win_az)
443  {
444  std::vector< std::vector< std::vector<HCA_Park> > > vol_Ai_filtered;
445  vol_Ai_filtered=vol_Ai;
446  unsigned half_rg=0.5*(win_rg-1);
447  unsigned half_az=0.5*(win_az-1);
448  int pre_rg,post_rg,pre_az,post_az;
449  unsigned count=0;
450  //cout<<half_az<<"\t"<<half_rg<<endl;
451  for(unsigned el=0;el<vol_Ai.size();el++)
452  {
453  vol_hca.push_back(PolarScan<EchoClass>(vol_z.scan(el).beam_count,vol_z.scan(el).beam_size, NC));
454  vol_hca[el].elevation=vol_z[el].elevation;
455  vol_hca[el].gain=1;
456  vol_hca[el].offset=0;
457  vol_hca[el].undetect=NC;
458  vol_hca[el].nodata=255;
459  for(unsigned az=0;az<vol_Ai[el].size();az++)
460  for(unsigned rg=0;rg<vol_Ai[el][az].size();rg++)
461  {
462  //cout<<el<<"\t"<<az<<"\t"<<rg<<endl;
463  vol_Ai_filtered[el][az][rg].clearAi();
464  vol_Ai_filtered[el][az][rg].Ai+=vol_Ai[el][az][rg].Ai;
465  count++;
466  //cout<<"faccio gli az"<<endl;
467  for(unsigned j=1;j<half_az+1;j++)
468  {
469  pre_az=az-j;
470  post_az=az+j;
471  //cout<<pre_az<<"\t"<<post_az<<endl;
472  if(pre_az<0)pre_az+=vol_Ai[el].size();
473  if(post_az>=vol_Ai[el].size())post_az-=vol_Ai[el].size();
474  //cout<<pre_az<<"\t"<<post_az<<endl;
475  vol_Ai_filtered[el][az][rg].Ai+=vol_Ai[el][pre_az][rg].Ai;
476  vol_Ai_filtered[el][az][rg].Ai+=vol_Ai[el][post_az][rg].Ai;
477  count+=2;
478  }
479  //cout<<"faccio gli rg"<<flush;
480  for(unsigned i=1;i<half_rg+1;i++)
481  {
482  pre_rg=rg-i;
483  post_rg=rg+i;
484  if(pre_rg>=0)
485  {
486  vol_Ai_filtered[el][az][rg].Ai+=vol_Ai[el][az][pre_rg].Ai;
487  count++;
488  }
489  if(post_rg<vol_Ai[el][az].size())
490  {
491  vol_Ai_filtered[el][az][rg].Ai+=vol_Ai[el][az][post_rg].Ai;
492  count++;
493  }
494  //cout<<"faccio gli az "<<i<<flush;
495  for(unsigned j=1;j<half_az+1;j++)
496  {
497  pre_az=az-j;
498  post_az=az+j;
499  pre_az=az-j;
500  post_az=az+j;
501  if(pre_az<0)pre_az+=vol_Ai[el].size();
502  if(post_az>=vol_Ai[el].size())post_az-=vol_Ai[el].size();
503  if(pre_rg>=0)
504  {
505  vol_Ai_filtered[el][az][rg].Ai+=vol_Ai[el][pre_az][pre_rg].Ai;
506  vol_Ai_filtered[el][az][rg].Ai+=vol_Ai[el][pre_az][pre_rg].Ai;
507  count+=2;
508  }
509  if(post_rg<vol_Ai[el][az].size())
510  {
511  vol_Ai_filtered[el][az][rg].Ai+=vol_Ai[el][post_az][post_rg].Ai;
512  vol_Ai_filtered[el][az][rg].Ai+=vol_Ai[el][post_az][post_rg].Ai;
513  count+=2;
514  }
515  }
516  }
517  //cout<<" fatto"<<endl;
518  vol_Ai_filtered[el][az][rg].Ai.array()/=(double)count;
519  //cout<<"normalizzato"<<endl;
520  vol_hca[el](az,rg)=vol_Ai_filtered[el][az][rg].echo(0.00001);
521  //cout<<"azzero"<<endl;
522  count=0;
523  }
524  }
525  }
526  else
527  {
528  for(unsigned el=0;el<vol_z.size();el++)
529  {
530  vol_hca.push_back(PolarScan<EchoClass>(vol_z.scan(el).beam_count,vol_z.scan(el).beam_size, NC));
531  vol_hca[el].elevation=vol_z[el].elevation;
532  for(unsigned az=0;az<vol_z.scan(el).beam_count;az++)
533  for(unsigned rg=0;rg<vol_z.scan(el).beam_size;rg++)
534  vol_hca[el](az,rg)=vol_Ai[el][az][rg].echo(0.00001);
535  }
536  }
537 }
538 
540 {
541  vector< vector< HCA_Park > > SCAN;
542  vector< HCA_Park > BEAM;
543  printf("inizio HCA\n");
544  vol_Ai.resize(vol_z.size());
545  double Z,Zdr,rhohv,lkdp,sdz,sdphidp,vrad;
546  double phidp,snr,gradphitheta,gradphiphi,gradZtheta,gradZphi,gradZdrtheta,gradZdrphi;
547  for(unsigned el=0;el<vol_z.size();el++)
548  {
549  cout<<"\tHCA el "<<el<<endl;
550  SCAN.resize(vol_z.scan(el).beam_count);
551  for(unsigned az=0;az<vol_z.scan(el).beam_count;az++)
552  {
553  BEAM.resize(vol_z.scan(el).beam_size);
554  for(unsigned rg=0;rg<vol_z.scan(el).beam_size;rg++)
555  {
556  Z=vol_z_1km.scan(el).get(az,rg);
557  Zdr=vol_zdr_2km.scan(el).get(az,rg);
558  rhohv=vol_rhohv_2km.scan(el).get(az,rg);
559  lkdp=Z>40?vol_lkdp_2km[el].get(az,rg):vol_lkdp_6km[el].get(az,rg);
560  sdz=vol_sdz.scan(el).get(az,rg);
561  sdphidp=vol_sdphidp.scan(el).get(az,rg);
562  vrad=vol_vrad.scan(el).get(az,rg);
563  phidp=vol_phidp[el](az,rg);
564  snr=vol_snr[el](az,rg);
565  gradphitheta=vol_grad_phi_theta[el](az,rg);
566  gradphiphi=vol_grad_phi_phi[el](az,rg);
567  gradZtheta=vol_grad_z_theta[el](az,rg);
568  gradZphi=vol_grad_z_phi[el](az,rg);
569  gradZdrtheta=vol_grad_zdr_theta[el](az,rg);
570  gradZdrphi=vol_grad_zdr_phi[el](az,rg);
571 
572  HCA_Park hca(Z,Zdr,rhohv,lkdp,sdz,sdphidp,vrad,phidp,snr,gradphitheta,gradphiphi,gradZtheta,gradZphi,gradZdrtheta,gradZdrphi);
573  BEAM[rg]=hca;
574  }
575  SCAN[az]=BEAM;
576  }
577  vol_Ai[el]=SCAN;
578  }
579  // Dopo aver calcolato i valori di aggregazione cerco il melting layer
581  cout<<"applico ML criteria ad HCA"<<endl;
583  class_designation(1,1);
584 /* unsigned elev=1;
585  unsigned azim=140;
586  cout<<"GC\tBS\tDS\tWS\tCR\tGR\tBD\tRA\tHR\tRH"<<endl;
587  for(unsigned rg=0;rg<vol_Ai[elev][azim].size();rg++)
588  {
589  cout.precision(5);
590  cout<<fixed<<vol_Ai[elev][azim][rg].Ai[GC_AP]<<"\t"<<vol_Ai[elev][azim][rg].Ai[BS]<<"\t"<<
591  vol_Ai[elev][azim][rg].Ai[DS]<<"\t"<<vol_Ai[elev][azim][rg].Ai[WS]<<"\t"<<
592  vol_Ai[elev][azim][rg].Ai[CR]<<"\t"<<vol_Ai[elev][azim][rg].Ai[GR]<<"\t"<<
593  vol_Ai[elev][azim][rg].Ai[BD]<<"\t"<<vol_Ai[elev][azim][rg].Ai[RA]<<"\t"<<
594  vol_Ai[elev][azim][rg].Ai[HR]<<"\t"<<vol_Ai[elev][azim][rg].Ai[RH]<<"\t"<<
595  vol_hca[elev](azim,rg)<<endl;
596  }
597 */
598 }
599 
601 {
602  gdal_init_once();
603  if(elev>=0)
604  {
605  ostringstream filename;
606  filename<<"class_"<<elev<<".png";
608  img=vol_hca[elev].cast<unsigned short>();
609  img*=6553;
610  write_image(img, filename.str(), "PNG");
611  }
612  else
613  {
614  for(unsigned el=0;el<vol_hca.size();el++)
615  {
616  ostringstream filename;
617  filename<<"class_"<<el<<".png";
619  img=vol_hca[el].cast<unsigned short>();
620  img*=6553;
621  write_image(img, filename.str(), "PNG");
622  }
623  }
624 }
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
PolarScan< T > & scan(unsigned idx)
Access a polar scan.
Definition: volume.h:298
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:440
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:419
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
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
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:600
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
void HCA_Park_2009()
Compute Echo Classes Park et al. (2009) HCA algorithm.
Definition: classifier.cpp:539