Elaboradar 0.1
Caricamento in corso...
Ricerca in corso...
Nessun risultato
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
27using namespace radarelab;
28using namespace volume;
29using namespace std;
30namespace odim = OdimH5v21;
31
32PROB::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
47Matrix2D<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
137double 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
145HCA_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
171classifier::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
446void 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
545void 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;
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{
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}
Classes to compute hydrometeor classification.
PolarScan - structure to describe a polarScan containing a matrix of data and conversion factors.
Definition: volume.h:113
const unsigned beam_count
Number of beam_count used ast each elevations.
Definition: volume.h:434
Homogeneous volume with a common beam count for all PolarScans.
Definition: volume.h:431
compute confidence vector of radar variables
Definition: classifier.h:182
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)
Definition: classifier.cpp:32
Given radar variables compute matrix of probability.
Definition: classifier.h:101
std::string units
Data units according to ODIM documentation.
Definition: volume.h:270
T offset
Data Offset.
Definition: volume.h:276
std::string quantity
Odim quantity name.
Definition: volume.h:268
std::shared_ptr< LoadInfo > load_info
Polar volume information.
Definition: volume.h:272
PolarScan< T > & scan(unsigned idx)
Access a polar scan.
Definition: volume.h:313
Sequence of PolarScans which can have a different beam count for each elevation.
Definition: volume.h:264
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
classifier(const std::string &file)
Constructor from odim file.
Definition: classifier.cpp:171
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
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
EchoClass
List classes of radar echoes Classes defined in Park et al. (2009) Overload << operator to print clas...
Definition: classifier.h:34
void gdal_init_once()
Initialize the GDAL library when called for the first time; does nothing all other times.
Definition: image.cpp:12
String functions.
Definition: cart.cpp:4
Codice per il caricamento di volumi ODIM in radarelab.
Base for all matrices we use, since we rely on row-major data.
Definition: matrix.h:37
std::map< std::string, Scans< double > * > to_load
Map used to specify quantity to be loaded.
Definition: loader.h:26
std::vector< double > get_nominal_elevations()
Get list of available elevation.
Definition: odim.h:41
void request_quantity(const std::string &name, Scans< double > *volume)
Define a request - Fill to_load attribute
Definition: odim.cpp:29
void load(const std::string &pathname)
Load method.
Definition: odim.cpp:34
Struttura che eredita da Loader e definisce i metodi per accedere ai dati ODIM.
Definition: odim.h:23