38 std::vector<bool>
Cleaner::clean_beam(
const Eigen::VectorXd& beam_z,
const Eigen::VectorXd& beam_w,
const Eigen::VectorXd& beam_v,
int iray)
const
40 const unsigned beam_size = beam_z.rows();
41 vector<bool> res(beam_size,
false);
42 bool in_a_segment =
false;
44 unsigned segment_length;
48 for (
unsigned ibin = 0; ibin < beam_size; ++ibin)
54 if (beam_w(ibin) == W_threshold && beam_v(ibin) == bin_wind_magic_number)
65 if (beam_w(ibin) != W_threshold || beam_v(ibin) != bin_wind_magic_number || ibin == (beam_size - 1))
69 if (ibin == (beam_size - 1)) end = ibin;
71 segment_length = end - start;
72 counter = counter + (unsigned)(segment_length);
77 for (
int ib = ibin - 12; ib < (signed)ibin; ++ib)
78 if (ib >= 0 && beam_z(ib) > Z_missing)
80 if (c_b > 0.25*12) before =
true;
83 for (
unsigned ia = ibin + 1; ia <= ibin + 12; ++ia)
84 if (ia < beam_size && beam_z(ia) >= Z_missing)
86 if (c_a > 0.25*12) after =
true;
89 if ((segment_length >= min_segment_length && !before && !after) ||
90 segment_length >= max_segment_length || counter > 100)
94 for (
unsigned ib = start; ib <= end; ++ib)
95 if( beam_z(ib) > Z_missing)res[ib] =
true;
105 std::vector<unsigned char>
Cleaner::eval_clean_beam(
const Eigen::VectorXd& beam_z,
const Eigen::VectorXd& beam_w,
const Eigen::VectorXd& beam_v,
int iray)
const
107 const unsigned beam_size = beam_z.rows();
108 vector<unsigned char> res(beam_size, 0);
109 bool in_a_segment =
false;
110 unsigned start = 0, end = 0;
111 unsigned segment_length;
112 bool before =
false, after =
false;
113 unsigned counter = 0;
115 for (
unsigned ibin = 0; ibin < beam_size; ++ibin)
121 if (beam_w(ibin) == W_threshold && beam_v(ibin) == bin_wind_magic_number)
131 if (beam_w(ibin) != W_threshold || beam_v(ibin) != bin_wind_magic_number || ibin == (beam_size - 1))
133 in_a_segment =
false;
135 if (ibin == (beam_size - 1)) end = ibin;
137 segment_length = end - start;
138 counter = counter + (unsigned)(segment_length);
143 for (
int ib = ibin - 12; ib < (signed)ibin; ++ib)
144 if (ib >= 0 && beam_z(ib) > Z_missing)
146 if (c_b > 0.25*12) before =
true;
149 for (
unsigned ia = ibin + 1; ia <= ibin + 12; ++ia)
150 if (ia < beam_size && beam_z(ia) >= Z_missing)
152 if (c_a > 0.25*12) after =
true;
155 if ((segment_length >= min_segment_length && !before && !after) ||
156 segment_length >= max_segment_length || counter > 100)
160 for (
unsigned ib = start; ib <= end; ++ib)
161 if( beam_z(ib) > Z_missing)res[ib] = 1;
177 const unsigned beam_size = beam_z.rows();
178 vector<bool> res(beam_size,
false);
179 bool in_a_segment =
false;
180 unsigned start = 0, end;
181 unsigned segment_length;
183 unsigned counter = 0;
184 unsigned counter_trash = 0;
185 unsigned counter_clutter =0;
186 for (
unsigned ibin = 0; ibin < beam_size; ++ibin)
188 bool is_clutter =
false;
189 bool is_trash =
false;
194 if ( beam_w(ibin) == W_threshold && beam_v(ibin) == bin_wind_magic_number && beam_z (ibin) != Z_missing ) {
196 if( beam_sdzdr(ibin) <= 0.01 ){
200 if (beam_z (ibin) >= 45. ){
202 if ((ibin >100 &&
double(counter_trash)/
double(ibin) >=0.5 && ( beam_sdzdr(ibin) >1 || beam_sd (ibin) > 5. )) ||
203 (beam_sdzdr(ibin) >4.0 && beam_sd (ibin) > 20.) ) {
210 }
else if ( (ibin >100 &&
double(counter_trash)/
double(ibin) >=0.5 && ( beam_sdzdr(ibin) >1 || beam_sd (ibin) > 5. )) ||
211 (beam_sd (ibin) >2. && (beam_sdzdr(ibin) >2.0 || beam_sd (ibin) > 10. )) ) {
219 if ((beam_w(ibin) * fabs(beam_v(ibin))) <= 0.25 && beam_z (ibin) != Z_missing ) {
224 if( is_clutter) counter_clutter ++;
225 if( is_trash ) counter_trash ++;
226 if(ibin <40 &&
false){
227 printf(
" %4d %4d %6.2f %6.2f %10.6f %6.2f %6.2f ",iray,ibin , beam_z(ibin),beam_v(ibin),beam_w(ibin), beam_sd(ibin),beam_sdzdr(ibin));
228 printf(
" ----- %2x %2x %2x %2x ",(
unsigned char)((beam_z(ibin)-scan_z.
offset)/scan_z.
gain/256),
229 (
unsigned char)((beam_v(ibin)-scan_v.
offset)/scan_v.
gain/256),
230 (
unsigned char)((beam_w(ibin)-scan_w.
offset)/scan_w.
gain/256),
231 (
unsigned char)((beam_sd(ibin)-SD.
offset)/SD.
gain/256));
236 if ( is_clutter || is_trash )
247 if ( ! (is_clutter || is_trash ) || ibin == (beam_size - 1))
249 in_a_segment =
false;
251 if (ibin == (beam_size - 1)) end = ibin;
253 segment_length = end - start+1;
254 counter = counter + (unsigned)(segment_length);
257 if (segment_length <= 2*min_segment_length ){
260 for (
int ib = ibin - 2*min_segment_length; ib < (signed)ibin; ++ib)
261 if (ib >= 0 && (beam_z(ib) > Z_missing && beam_w(ib) != W_threshold && ( beam_w(ib) > 0.5 || fabs(beam_v(ib)) > 0.5) ) )
263 if (
double(count)/double(min(
int(ibin),
int(2*min_segment_length))) >=0.25) before =
true;
267 for (
unsigned ia = ibin + 1; ia <= ibin + 2*min_segment_length; ++ia)
268 if (ia < beam_size && (beam_z(ia) > Z_missing && (beam_w(ia) != W_threshold && ( beam_w(ia) > 0.5 || fabs(beam_v(ia)) > 0.5)) ))
270 if (
double(count)/double(min(
int(beam_size - ibin),
int(2*min_segment_length))) >=0.25) after =
true;
273 if ((segment_length >= min_segment_length && (!before || !after) ) || segment_length >= max_segment_length)
278 for (
unsigned ib = start; ib <= end; ++ib)
296 const unsigned beam_size = beam_z.rows();
297 vector<bool> res(beam_size,
false);
298 bool in_a_segment =
false;
300 unsigned segment_length;
302 unsigned counter = 0;
304 for (
unsigned ibin = 0; ibin < beam_size; ++ibin)
314 if ( ((beam_w(ibin) == W_threshold && beam_v(ibin) == bin_wind_magic_number) ||(beam_w(ibin) * fabs(beam_v(ibin)) <= 0.25) ) && beam_z (ibin) != Z_missing && beam_sd(ibin) > sd_threshold )
325 if ( ( ( beam_w(ibin) != W_threshold || beam_v(ibin) != bin_wind_magic_number) && (beam_w(ibin) * fabs(beam_v(ibin)) > 0.25) ) || ibin == (beam_size - 1) || beam_z(ibin) == Z_missing || beam_sd(ibin) <= sd_threshold )
327 in_a_segment =
false;
329 if (ibin == (beam_size - 1)) end = ibin;
331 segment_length = end - start+1;
332 counter = counter + (unsigned)(segment_length);
335 if (segment_length <= 2*min_segment_length ){
338 for (
int ib = ibin - 2*min_segment_length; ib < (signed)ibin; ++ib)
339 if (ib >= 0 && (beam_z(ib) > Z_missing && beam_w(ib) != W_threshold && ( beam_w(ib) > 0.5 || fabs(beam_v(ib)) > 0.5) ) )
341 if (
double(count)/double(min(
int(ibin),
int(2*min_segment_length))) >=0.25) before =
true;
345 for (
unsigned ia = ibin + 1; ia <= ibin + 2*min_segment_length; ++ia)
346 if (ia < beam_size && (beam_z(ia) > Z_missing && (beam_w(ia) != W_threshold && ( beam_w(ia) > 0.5 || fabs(beam_v(ia)) > 0.5)) ))
348 if (
double(count)/double(min(
int(beam_size - ibin),
int(2*min_segment_length))) >=0.25) after =
true;
351 if ((segment_length >= min_segment_length && (!before || !after) ) || segment_length >= max_segment_length)
356 for (
unsigned ib = start; ib <= end; ++ib)
369 tuple<std::vector<unsigned char>,std::vector<double>> Cleaner::eval_classID_beam(
const Eigen::VectorXd& beam_z,
const Eigen::VectorXd& beam_w,
const Eigen::VectorXd& beam_v,
const Eigen::VectorXd& beam_sd,
const Eigen::VectorXd& beam_zdr,
const Eigen::VectorXd& beam_rohv,
const Eigen::VectorXd& beam_sqi,
const Eigen::VectorXd& beam_snr,
const Eigen::VectorXd& beam_zvd,
const Eigen::VectorXd& beam_sdray,
const Eigen::VectorXd& beam_sdaz,
const Eigen::VectorXd& beam_zdr_sd,
int iray,
const string radar,
double v_ny,
const char* fuzzy_path,
bool stamp,
bool force_meteo)
const
372 const unsigned beam_size = beam_z.rows();
373 vector<unsigned char> res(beam_size, 0);
374 vector<double> diffs(beam_size, 0);
384 string f_dir = fuzzy_path;
388 string fin_w = f_dir+
"/matrix-"+radar+
".txt";
390 vector<string> w_vector;
391 w_vector = read_matrix_from_txt(fin_w);
392 Num_entries = w_vector.size()/Num_echoes;
395 for(
int i=0;i<Num_echoes;i++){
396 for(
int j=0;j<Num_entries;j++){
397 Wij(i,j) = stod( w_vector[i*Num_entries+j]);
405 string fin_t = f_dir+
"/Trap-"+radar+
".txt";
406 vector<string> t_vector;
407 t_vector = read_matrix_from_txt(fin_t);
408 double Traps[Num_entries][Num_echoes][Ntraps];
409 for(
int i=0;i<Num_entries;i++){
410 for(
int j=0;j<Num_echoes;j++){
411 for(
int k=0;k<Ntraps;k++){
412 Traps[i][j][k] = stod( t_vector[i*(Num_echoes*Ntraps)+j*Num_echoes+k]);
422 vector<unsigned> counter (Num_entries,0) ;
423 for (
unsigned ibin = 0; ibin < beam_size; ++ibin)
430 ArrayXd Class_WP(Num_echoes);
432 if (beam_z(ibin) == Z_missing) {
443 if (beam_v(ibin) != bin_wind_magic_number ) {
446 double pv0 = trap(v_ny*Traps[0][0][0], v_ny*Traps[0][0][1], v_ny*Traps[0][0][2], v_ny*Traps[0][0][3], beam_v(ibin));
447 double pv3 = trap(v_ny*Traps[0][0][4],v_ny*Traps[0][0][4],v_ny*Traps[0][0][4],v_ny*Traps[0][0][4],beam_v(ibin));
448 Pij(0,1)=std::max(pv0,pv3);
449 }
else Pij(0,1)=trap(v_ny*Traps[0][0][0], v_ny*Traps[0][0][1], v_ny*Traps[0][0][2], v_ny*Traps[0][0][3], beam_v(ibin),v_ny*Traps[0][0][4]);
451 double prob_v = trap (v_ny, v_ny, v_ny, v_ny, beam_v(ibin));
455 double pv0 = trap(v_ny*Traps[0][1][0], v_ny*Traps[0][1][1], v_ny*Traps[0][1][2], v_ny*Traps[0][1][3], beam_v(ibin));
456 double pv1 = trap(-8.,-8.,-7.,-7.,beam_v(ibin));
457 double pv2 = trap(5.5,5.5,6.5,6.5,beam_v(ibin));
458 double pv3 = trap(v_ny*Traps[0][1][4],v_ny*Traps[0][1][4],v_ny*Traps[0][1][4],v_ny*Traps[0][1][4],beam_v(ibin));
459 Pij(1,1)=std::max({pv0,pv1,pv2,pv3});
460 }
else Pij(1,1)=trap(v_ny*Traps[0][1][0], v_ny*Traps[0][1][1], v_ny*Traps[0][1][2], v_ny*Traps[0][1][3], beam_v(ibin),v_ny*Traps[0][1][4]);
474 Pij(0,2)= trap(Traps[1][0][0],Traps[1][0][1],Traps[1][0][2],Traps[1][0][3],beam_w(ibin));
475 Pij(1,2)= trap(Traps[1][1][0],Traps[1][1][1],Traps[1][1][2],Traps[1][1][3],beam_w(ibin));
476 double prob_w = trap(Traps[1][2][0],Traps[1][2][1],Traps[1][2][2],Traps[1][2][3],beam_w(ibin));
482 Pij(0,0) = trap(Traps[2][0][0],Traps[2][0][1],Traps[2][0][2],Traps[2][0][3],beam_z(ibin), Traps[2][0][4]);
483 Pij(0,3) = trap (Traps[3][0][0],Traps[3][0][1],Traps[3][0][2],Traps[3][0][3], beam_sd(ibin));
484 Pij(0,4) = trap (Traps[4][0][0],Traps[4][0][1],Traps[4][0][2],Traps[4][0][3], beam_sdray(ibin));
485 Pij(0,5) = trap (Traps[5][0][0],Traps[5][0][1],Traps[5][0][2],Traps[5][0][3], beam_sdaz(ibin));
487 Pij(0,6) = trap(Traps[6][0][0],Traps[6][0][1],Traps[6][0][2],Traps[6][0][3],beam_zdr(ibin));
488 Pij(0,7) = trap(Traps[7][0][0],Traps[7][0][1],Traps[7][0][2],Traps[7][0][3],beam_rohv(ibin),Traps[7][0][4]);
489 Pij(0,8) = trap(Traps[8][0][0],Traps[8][0][1],Traps[8][0][2],Traps[8][0][3],beam_zdr_sd(ibin));
490 Pij(0,9) = trap(Traps[9][0][0],Traps[9][0][1],Traps[9][0][2],Traps[9][0][3], beam_sqi(ibin));
491 Pij(0,10) = trap(Traps[10][0][0],Traps[10][0][1],Traps[10][0][2],Traps[10][0][3], beam_snr(ibin));
492 if((radar==
"GAT")&&(ibin<=120))
493 Pij(0,11) = trap(Traps[11][0][0],0.,Traps[11][0][2],Traps[11][0][3],beam_zvd(ibin));
495 Pij(0,11) = trap(Traps[11][0][0],Traps[11][0][1],Traps[11][0][2],Traps[11][0][3],beam_zvd(ibin));
498 Pij(1,0) = trap (Traps[2][1][0],Traps[2][1][1],Traps[2][1][2],Traps[2][1][3], beam_z(ibin), Traps[2][1][4]);
499 Pij(1,3) = trap (Traps[3][1][0],Traps[3][1][1],Traps[3][1][2],Traps[3][1][3], beam_sd(ibin));
500 Pij(1,4) = trap (Traps[4][1][0],Traps[4][1][1],Traps[4][1][2],Traps[4][1][3], beam_sdray(ibin));
501 Pij(1,5) = trap (Traps[5][1][0],Traps[5][1][1],Traps[5][1][2],Traps[5][1][3], beam_sdaz(ibin));
503 Pij(1,6) = trap(Traps[6][1][0],Traps[6][1][1],Traps[6][1][2],Traps[6][1][3],beam_zdr(ibin),Traps[6][1][4] );
504 Pij(1,7) = trap(Traps[7][1][0],Traps[7][1][1],Traps[7][1][2],Traps[7][1][3],beam_rohv(ibin),Traps[7][1][4]);
505 Pij(1,8) = trap(Traps[8][1][0],Traps[8][1][1],Traps[8][1][2],Traps[8][1][3],beam_zdr_sd(ibin));
506 Pij(1,9) = trap(Traps[9][1][0],Traps[9][1][1],Traps[9][1][2],Traps[9][1][3],beam_sqi(ibin));
507 Pij(1,10) = trap(Traps[10][1][0],Traps[10][1][1],Traps[10][1][2],Traps[10][1][3], beam_snr(ibin));
509 Pij(1,11) = std::max(trap(Traps[11][1][0],Traps[11][1][1],Traps[11][1][2],Traps[11][1][3],beam_zvd(ibin)),trap(35.,40.,60.,60.,beam_zvd(ibin)));
511 Pij(1,11) = trap(Traps[11][1][0],Traps[11][1][1],Traps[11][1][2],Traps[11][1][3],beam_zvd(ibin));
514 Pij(2,0) = trap(Traps[2][2][0],Traps[2][2][1],Traps[2][2][2],Traps[2][2][3],beam_z(ibin), Traps[2][2][4]);
515 Pij(2,3) = trap (Traps[3][2][0],Traps[3][2][1],Traps[3][2][2],Traps[3][2][3], beam_sd(ibin));
516 Pij(2,4) = trap (Traps[4][2][0],Traps[4][2][1],Traps[4][2][2],Traps[4][2][3], beam_sdray(ibin));
517 Pij(2,5) = trap (Traps[5][2][0],Traps[5][2][1],Traps[5][2][2],Traps[5][2][3], beam_sdaz(ibin));
519 Pij(2,6) = trap(Traps[6][2][0],Traps[6][2][1],Traps[6][2][2],Traps[6][2][3],beam_zdr(ibin),Traps[6][2][4]);
520 Pij(2,7) = trap(Traps[7][2][0],Traps[7][2][1],Traps[7][2][2],Traps[7][2][3],beam_rohv(ibin));
521 Pij(2,8) = trap(Traps[8][2][0],Traps[8][2][1],Traps[8][2][2],Traps[8][2][3],beam_zdr_sd(ibin));
522 Pij(2,9) = trap(Traps[9][2][0],Traps[9][2][1],Traps[9][2][2],Traps[9][2][3],beam_sqi(ibin),Traps[9][2][4]);
523 Pij(2,10) = trap(Traps[10][2][0],Traps[10][2][1],Traps[10][2][2],Traps[10][2][3], beam_snr(ibin));
525 Pij(2,11) = std::max(trap(Traps[11][2][0],Traps[11][2][1],Traps[11][2][2],Traps[11][2][3],beam_zvd(ibin)),trap(10.,20.,40.,55.,beam_zvd(ibin)));
527 Pij(2,11) = trap(Traps[11][1][0],Traps[11][1][1],Traps[11][1][2],Traps[11][1][3],beam_zvd(ibin));
529 Pij(3,0) = trap(Traps[2][3][0],Traps[2][3][1],Traps[2][3][2],Traps[2][3][3],beam_z(ibin), Traps[2][3][4]);
530 Pij(3,3) = trap (Traps[3][3][0],Traps[3][3][1],Traps[3][3][2],Traps[3][3][3], beam_sd(ibin));
531 Pij(3,4) = trap (Traps[4][3][0],Traps[4][3][1],Traps[4][3][2],Traps[4][3][3], beam_sdray(ibin));
532 Pij(3,5) = trap (Traps[5][3][0],Traps[5][3][1],Traps[5][3][2],Traps[5][3][3], beam_sdaz(ibin));
534 Pij(3,6) = trap(Traps[6][3][0],Traps[6][3][1],Traps[6][3][2],Traps[6][3][3],beam_zdr(ibin),Traps[6][3][4]);
535 Pij(3,7) = trap(Traps[7][3][0],Traps[7][3][1],Traps[7][3][2],Traps[7][3][3],beam_rohv(ibin));
536 Pij(3,8) = trap(Traps[8][3][0],Traps[8][3][1],Traps[8][3][2],Traps[8][3][3],beam_zdr_sd(ibin));
537 Pij(3,9) = trap(Traps[9][3][0],Traps[9][3][1],Traps[9][3][2],Traps[9][3][3],beam_sqi(ibin),Traps[9][3][4]);
538 Pij(3,10) = trap(Traps[10][3][0],Traps[10][3][1],Traps[10][3][2],Traps[10][3][3], beam_snr(ibin));
539 Pij(3,11) = trap(Traps[11][3][0],Traps[11][3][1],Traps[11][3][2],Traps[11][3][3],beam_zvd(ibin));
541 Pij(4,0) = trap(Traps[2][4][0],Traps[2][4][1],Traps[2][4][2],Traps[2][4][3],beam_z(ibin), Traps[2][4][4]);
542 Pij(4,3) = trap (Traps[3][4][0],Traps[3][4][1],Traps[3][4][2],Traps[3][4][3], beam_sd(ibin));
543 Pij(4,4) = trap (Traps[4][4][0],Traps[4][4][1],Traps[4][4][2],Traps[4][4][3], beam_sdray(ibin));
544 Pij(4,5) = trap (Traps[5][4][0],Traps[5][4][1],Traps[5][4][2],Traps[5][4][3], beam_sdaz(ibin));
546 Pij(4,6) = trap(Traps[6][4][0],Traps[6][4][1],Traps[6][4][2],Traps[6][4][3],beam_zdr(ibin),Traps[6][4][4]);
547 Pij(4,7) = trap(Traps[7][4][0],Traps[7][4][1],Traps[7][4][2],Traps[7][4][3],beam_rohv(ibin),Traps[7][4][4]);
548 Pij(4,8) = trap(Traps[8][4][0],Traps[8][4][1],Traps[8][4][2],Traps[8][4][3],beam_zdr_sd(ibin));
549 Pij(4,9) = trap(Traps[9][4][0],Traps[9][4][1],Traps[9][4][2],Traps[9][4][3],beam_sqi(ibin),Traps[9][4][4]);
550 Pij(4,10) = trap(Traps[10][4][0],Traps[10][4][1],Traps[10][4][2],Traps[10][4][3], beam_snr(ibin));
551 Pij(4,11) = trap(Traps[11][4][0],Traps[11][4][1],Traps[11][4][2],Traps[11][4][3],beam_zvd(ibin));
555 Class_WP = ((Wij.array()*Pij.array()).matrix()*VectorXd::Ones(Num_entries)).array()/(Wij*VectorXd::Ones(Num_entries)).array();
557 Class_WP.maxCoeff(&i);
559 if (Class_WP(i) < 0.1 ) ID=5;
561 double diff = Class_WP[i]-Class_WP[0];
566 printf(
"bin %4d ",ibin);
567 printf(
"ID %d -- diff=%8.3f\n",ID, diff);
570 for(
unsigned c=0;c<5;c++) printf(
"%5.3f ",Class_WP(c));
571 printf(
" ID %d \n",ID);
572 printf(
" %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f \n",beam_z(ibin), beam_v(ibin), beam_w(ibin), beam_sd(ibin), beam_sdray(ibin), beam_sdaz(ibin), beam_zdr(ibin),
573 beam_rohv(ibin), beam_zdr_sd(ibin), beam_sqi(ibin), beam_snr(ibin), beam_zvd(ibin));
574 for (
unsigned c=0;c<5; c++){
575 for (
unsigned d=0; d<12; d++) printf(
" %8.3f", Pij(c,d));
584 if((force_meteo)&&(diff<=8.0)&&(prevID==0)&&(ID!=1)) res[ibin] = 0;
591 else if(radar==
"SPC"){
592 if((force_meteo)&&(diff<=10.)&&(prevID==0)) res[ibin] = 0;
594 if((force_meteo)&&(diff<=20.0)&&(prevID==0)&&(ID==4)) res[ibin] = 0;
595 }
else cout<<
"radar non specificato"<<endl;
607 std::vector<unsigned char> Cleaner::eval_classID_beam(
const Eigen::VectorXd& beam_z,
const Eigen::VectorXd& beam_w,
const Eigen::VectorXd& beam_v,
const Eigen::VectorXd& beam_sd,
const Eigen::VectorXd& beam_sdray,
const Eigen::VectorXd& beam_sdaz,
int iray,
const string radar,
double v_ny,
const char* fuzzy_path)
const
610 const unsigned beam_size = beam_z.rows();
611 vector<unsigned char> res(beam_size, 0);
617 string f_dir = fuzzy_path;
620 string fin_w = f_dir+
"/matrix-"+radar+
"-nozdr.txt";
621 vector<string> w_vector;
622 w_vector = read_matrix_from_txt(fin_w);
623 Num_entries = w_vector.size()/Num_echoes;
626 for(
int i=0;i<Num_echoes;i++){
627 for(
int j=0;j<Num_entries;j++){
628 Wij(i,j) = stod( w_vector[i*Num_entries+j]);
633 string fin_t = f_dir+
"/Trap-"+radar+
"-nozdr.txt";
634 vector<string> t_vector;
635 t_vector = read_matrix_from_txt(fin_t);
636 double Traps[Num_entries][Num_echoes][Ntraps];
637 for(
int i=0;i<Num_entries;i++){
638 for(
int j=0;j<Num_echoes;j++){
639 for(
int k=0;k<Ntraps;k++){
640 Traps[i][j][k] = stod( t_vector[i*(Num_echoes*Ntraps)+j*Num_echoes+k]);
649 vector<unsigned> counter (Num_entries,0) ;
650 for (
unsigned ibin = 0; ibin < beam_size; ++ibin)
657 ArrayXd Class_WP(Num_echoes);
659 if (beam_z(ibin) == Z_missing) {
668 if (beam_v(ibin) != bin_wind_magic_number ) {
669 Pij(0,1)=trap(v_ny*Traps[0][0][0], v_ny*Traps[0][0][1], v_ny*Traps[0][0][2], v_ny*Traps[0][0][3], beam_v(ibin));
670 double prob_v = trap (v_ny, v_ny, v_ny, v_ny, beam_v(ibin),v_ny);
673 Pij(1,1)=trap(Traps[0][1][0],Traps[0][1][1],Traps[0][1][2],Traps[0][1][3], beam_v(ibin));
686 Pij(0,2)= trap(Traps[1][0][0],Traps[1][0][1],Traps[1][0][2],Traps[1][0][3],beam_w(ibin));
687 Pij(1,2)= trap(Traps[1][1][0],Traps[1][1][1],Traps[1][1][2],Traps[1][1][3],beam_w(ibin));
688 double prob_w = trap(Traps[1][2][0],Traps[1][2][1],Traps[1][2][2],Traps[1][2][3],beam_w(ibin));
694 Pij(0,0) = trap(Traps[2][0][0],Traps[2][0][1],Traps[2][0][2],Traps[2][0][3],beam_z(ibin), Traps[2][0][4]);
695 Pij(0,3) = trap (Traps[3][0][0],Traps[3][0][1],Traps[3][0][2],Traps[3][0][3], beam_sd(ibin));
696 Pij(0,4) = trap (Traps[4][0][0],Traps[4][0][1],Traps[4][0][2],Traps[4][0][3], beam_sdray(ibin));
697 Pij(0,5) = trap (Traps[5][0][0],Traps[5][0][1],Traps[5][0][2],Traps[5][0][3], beam_sdaz(ibin));
700 Pij(1,0) = trap (Traps[2][1][0],Traps[2][1][1],Traps[2][1][2],Traps[2][1][3], beam_z(ibin), Traps[2][1][4]);
701 Pij(1,3) = trap (Traps[3][1][0],Traps[3][1][1],Traps[3][1][2],Traps[3][1][3], beam_sd(ibin));
702 Pij(1,4) = trap (Traps[4][1][0],Traps[4][1][1],Traps[4][1][2],Traps[4][1][3], beam_sdray(ibin));
703 Pij(1,5) = trap (Traps[5][1][0],Traps[5][1][1],Traps[5][1][2],Traps[5][1][3], beam_sdaz(ibin));
705 Pij(2,0) = trap(Traps[2][2][0],Traps[2][2][1],Traps[2][2][2],Traps[2][2][3],beam_z(ibin), Traps[2][2][4]);
706 Pij(2,3) = trap (Traps[3][2][0],Traps[3][2][1],Traps[3][2][2],Traps[3][2][3], beam_sd(ibin));
707 Pij(2,4) = trap (Traps[4][2][0],Traps[4][2][1],Traps[4][2][2],Traps[4][2][3], beam_sdray(ibin));
708 Pij(2,5) = trap (Traps[5][2][0],Traps[5][2][1],Traps[5][2][2],Traps[5][2][3], beam_sdaz(ibin));
710 Pij(3,0) = trap(Traps[2][3][0],Traps[2][3][1],Traps[2][3][2],Traps[2][3][3],beam_z(ibin), Traps[2][3][4]);
711 Pij(3,3) = trap (Traps[3][3][0],Traps[3][3][1],Traps[3][3][2],Traps[3][3][3], beam_sd(ibin));
712 Pij(3,4) = trap (Traps[4][3][0],Traps[4][3][1],Traps[4][3][2],Traps[4][3][3], beam_sdray(ibin));
713 Pij(3,5) = trap (Traps[5][3][0],Traps[5][3][1],Traps[5][3][2],Traps[5][3][3], beam_sdaz(ibin));
715 Pij(4,0) = trap(Traps[2][4][0],Traps[2][4][1],Traps[2][4][2],Traps[2][4][3],beam_z(ibin), Traps[2][4][4]);
716 Pij(4,3) = trap (Traps[3][4][0],Traps[3][4][1],Traps[3][4][2],Traps[3][4][3], beam_sd(ibin));
717 Pij(4,4) = trap (Traps[4][4][0],Traps[4][4][1],Traps[4][4][2],Traps[4][4][3], beam_sdray(ibin));
718 Pij(4,5) = trap (Traps[5][4][0],Traps[5][4][1],Traps[5][4][2],Traps[5][4][3], beam_sdaz(ibin));
723 Class_WP = ((Wij.array()*Pij.array()).matrix()*VectorXd::Ones(Num_entries)).array()/(Wij*VectorXd::Ones(Num_entries)).array();
725 Class_WP.maxCoeff(&i);
727 if (Class_WP(i) < 0.1 ) ID=5;
739 std::vector<unsigned char>
Cleaner::eval_clean_beam(
const Eigen::VectorXd& beam_z,
const Eigen::VectorXd& beam_w,
const Eigen::VectorXd& beam_v,
const Eigen::VectorXd& beam_sd,
const Eigen::VectorXd& beam_sdray,
const Eigen::VectorXd& beam_sdaz,
int iray)
const
742 const unsigned beam_size = beam_z.rows();
743 vector<unsigned char> res(beam_size, 0);
744 unsigned counter = 0;
745 unsigned countIntWeak = 0;
746 unsigned countIntStrong = 0;
747 unsigned countIntMed = 0;
748 unsigned countClutter = 0;
749 unsigned countNoise = 0;
751 for (
unsigned ibin = 0; ibin < beam_size; ++ibin)
753 printf(
" %4d %4d %6.2f %6.2f %10.6f %10.6f %10.6f %10.6f",iray,ibin , beam_z(ibin),beam_v(ibin),beam_w(ibin),beam_sd(ibin),beam_sdray(ibin),beam_sdaz(ibin));
754 if ((beam_w(ibin) > W_threshold && beam_v(ibin) != bin_wind_magic_number && beam_sd(ibin) >= 1. &&
755 beam_sd(ibin) <= 10. ) || beam_z(ibin) == Z_missing) {
761 if (beam_w(ibin) == W_threshold && beam_v(ibin) == bin_wind_magic_number)
763 if (beam_sd(ibin) >= 5. && beam_sdray(ibin) >= 2 && beam_sdaz(ibin) > 4.){
768 if (beam_sd(ibin) >= 1. && beam_sd(ibin) <= 5. && beam_sdray(ibin) < 2. && beam_sdaz(ibin) > 4.){
773 if (beam_sd(ibin) >= 3. && beam_sd(ibin) <= 7. && beam_sdray(ibin) < 2. && beam_sdaz(ibin) > 4.&& beam_sdaz(ibin) <7 ){
779 if (beam_sd(ibin) >= 5. && beam_sdray(ibin) < 3. && beam_sdaz(ibin) < 3.){
784 if (beam_sd(ibin) >= 10. && beam_sdray(ibin) < 2. && beam_sdaz(ibin) < 5. && beam_z(ibin) <10.){
791 printf(
"%2d %4d %4d %4d %4d %4d\n",res[ibin],countClutter,countIntStrong, countIntMed, countIntWeak, countNoise);
803 return evaluateCleanID(scan_z, scan_w, scan_v, scan_cleanID, scan_v.
undetect,iel);
811 throw std::runtime_error(
"scan_z beam_count no equal to scan_w beam_count");
813 throw std::runtime_error(
"scan_z beam_size no equal to scan_w beam_size");
816 throw std::runtime_error(
"scan_z beam_count no equal to scan_v beam_count");
818 throw std::runtime_error(
"scan_z beam_size no equal to scan_v beam_size");
822 const unsigned beam_count = scan_z.
beam_count;
823 const unsigned beam_size = scan_z.
beam_size;
833 Z_S.push_back(scan_z);
834 radarelab::volume::textureSD( Z_S,SD2D, 1000. , 3,
false);
836 radarelab::volume::textureSD( Z_S,SD_Ray, scan_z.
cell_size*9 , 360./scan_z.
beam_count,
true);
837 radarelab::volume::textureSD( Z_S,SD_Az, scan_z.
cell_size , 5*360./scan_z.
beam_count,
true);
839 for (
unsigned i = 0; i < beam_count; ++i)
842 vector<unsigned char> corrected = cleaner.eval_clean_beam(scan_z.row(i), scan_w.row(i), scan_v.row(i),SD2D[0].row(i), SD_Ray[0].row(i), SD_Az[0].row(i), i);
844 for (
unsigned ib = 0; ib < beam_size; ++ib)
845 scan_cleanID(i,ib)=corrected[ib];
849 void Cleaner::evaluateClassID(
PolarScan<double>& scan_z,
PolarScan<double>& scan_w,
PolarScan<double>& scan_v,
PolarScan<double>& scan_zdr,
PolarScan<double>& scan_rohv,
PolarScan<double>& scan_sqi,
PolarScan<double>& scan_snr,
PolarScan<double>& scan_zvd,
PolarScan<unsigned char>& scan_cleanID,
PolarScan<double>& scan_DiffProb,
double bin_wind_magic_number,
const string radar,
const char* fuzzy_path,
unsigned iel,
bool force_meteo)
853 throw std::runtime_error(
"scan_z beam_count no equal to scan_w beam_count");
855 throw std::runtime_error(
"scan_z beam_size no equal to scan_w beam_size");
858 throw std::runtime_error(
"scan_z beam_count no equal to scan_v beam_count");
860 throw std::runtime_error(
"scan_z beam_size no equal to scan_v beam_size");
864 const unsigned beam_count = scan_z.
beam_count;
865 const unsigned beam_size = scan_z.
beam_size;
867 cout<<
"VRAD OFFSET="<<scan_v.
offset<<
" , GAIN="<<scan_v.
gain<<endl;
871 Z_S.push_back(scan_z);
872 ZDR_S.push_back(scan_zdr);
873 radarelab::volume::textureSD( Z_S,SD2D, 1000. , 3,
false);
874 radarelab::volume::textureSD( Z_S,SD_Ray, scan_z.
cell_size*21 , 360./scan_z.
beam_count,
true);
875 radarelab::volume::textureSD( Z_S,SD_Az, scan_z.
cell_size , 5*360./scan_z.
beam_count,
true);
877 radarelab::volume::textureSD( ZDR_S,ZDR_SD2D, 1000. , 3,
false);
880 for (
unsigned i = 0; i <beam_count ; ++i)
885 auto [corrected, diff_prob] = cleaner.eval_classID_beam(scan_z.row(i), scan_w.row(i), scan_v.row(i), SD2D[0].row(i), scan_zdr.row(i), scan_rohv.row(i), scan_sqi.row(i), scan_snr.row(i), scan_zvd.row(i), SD_Ray[0].row(i), SD_Az[0].row(i), ZDR_SD2D[0].row(i), i, radar, scan_v.
offset, fuzzy_path, stamp, force_meteo);
887 for (
unsigned ib = 0; ib < beam_size; ++ib){
888 scan_cleanID(i,ib)=corrected[ib];
889 scan_DiffProb(i,ib)=diff_prob[ib];
901 throw std::runtime_error(
"scan_z beam_count no equal to scan_w beam_count");
903 throw std::runtime_error(
"scan_z beam_size no equal to scan_w beam_size");
906 throw std::runtime_error(
"scan_z beam_count no equal to scan_v beam_count");
908 throw std::runtime_error(
"scan_z beam_size no equal to scan_v beam_size");
912 const unsigned beam_count = scan_z.
beam_count;
913 const unsigned beam_size = scan_z.
beam_size;
917 Z_S.push_back(scan_z);
918 radarelab::volume::textureSD( Z_S,SD2D, 1000. , 3,
false);
919 radarelab::volume::textureSD( Z_S,SD_Ray, scan_z.
cell_size*21 , 360./scan_z.
beam_count,
true);
920 radarelab::volume::textureSD( Z_S,SD_Az, scan_z.
cell_size , 5*360./scan_z.
beam_count,
true);
923 for (
unsigned i = 0; i <beam_count ; ++i)
925 vector<unsigned char> corrected = cleaner.eval_classID_beam(scan_z.row(i), scan_w.row(i), scan_v.row(i), SD2D[0].row(i), SD_Ray[0].row(i), SD_Az[0].row(i), i, radar, scan_v.
offset, fuzzy_path);
926 for (
unsigned ib = 0; ib < beam_size; ++ib)
927 scan_cleanID(i,ib)=corrected[ib];
933 return clean(scan_z, scan_w, scan_v, scan_v.
undetect,iel);
939 throw std::runtime_error(
"scan_z beam_count no equal to scan_w beam_count");
941 throw std::runtime_error(
"scan_z beam_size no equal to scan_w beam_size");
944 throw std::runtime_error(
"scan_z beam_count no equal to scan_v beam_count");
946 throw std::runtime_error(
"scan_z beam_size no equal to scan_v beam_size");
950 const unsigned beam_count = scan_z.
beam_count;
951 const unsigned beam_size = scan_z.
beam_size;
957 Z_S.push_back(scan_z);
958 radarelab::volume::textureSD( Z_S,SD2D, 1000. , 3,
false);
959 radarelab::volume::textureSD( Z_S,SD_Ray, scan_z.
cell_size*9 , 360./scan_z.
beam_count,
false);
960 radarelab::volume::textureSD( Z_S,SD_Az, scan_z.
cell_size , 5*360./scan_z.
beam_count,
false);
991 if (set_undetect) new_val=scan_z.
nodata;
993 for (
unsigned i = 0; i < beam_count; ++i)
997 vector<bool> corrected = cleaner.
clean_beam(scan_z.row(i), scan_w.row(i), scan_v.row(i),i);
999 for (
unsigned ib = 0; ib < beam_size; ++ib)
1003 scan_z(i, ib) = new_val;
1020 return clean(scan_z, scan_w, scan_v, scan_zdr, scan_v.
undetect,iel,set_undetect);
1024 cout<<
"Chiamato cleaner "<<set_undetect<<endl;
1026 throw std::runtime_error(
"scan_z beam_count no equal to scan_w beam_count");
1028 throw std::runtime_error(
"scan_z beam_size no equal to scan_w beam_size");
1031 throw std::runtime_error(
"scan_z beam_count no equal to scan_v beam_count");
1033 throw std::runtime_error(
"scan_z beam_size no equal to scan_v beam_size");
1034 double z_val=scan_z.
nodata;
1035 if(set_undetect) z_val=scan_z.
undetect;
1039 const unsigned beam_count = scan_z.
beam_count;
1040 const unsigned beam_size = scan_z.
beam_size;
1046 Z_S.push_back(scan_z);
1047 radarelab::volume::textureSD( Z_S,SD2D, 1000. , 3,
false);
1049 ZDR_S.push_back(scan_zdr);
1050 radarelab::volume::textureSD( ZDR_S,SDZDR2D, 1000. , 3,
false);
1066 printf(
"scrivo Z ");
1067 img = (scan_z.array() - scan_z.
offset )/ scan_z.
gain /256 ;
1068 sprintf(pippo,
"_%02d.png",iel);
1070 img_tmp=img.cast<
unsigned char>();
1072 radarelab::write_image(img_tmp,
"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_Z"+ext,
"PNG");
1083 img = (SD2D[0].array()-SD2D[0].
offset)/SD2D[0].gain/256 ;
1084 img_tmp=img.cast<
unsigned char>();
1085 radarelab::write_image(img_tmp,
"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_SD2d"+ext,
"PNG");
1087 img = (SDZDR2D[0].array()-SDZDR2D[0].
offset)/SDZDR2D[0].gain/256 ;
1088 img_tmp=img.cast<
unsigned char>();
1089 radarelab::write_image(img_tmp,
"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_SDZDR2d"+ext,
"PNG");
1092 double new_val=cleaner.Z_missing;
1093 if (set_undetect) new_val=scan_z.
undetect;
1095 for (
unsigned i = 0; i < beam_count; ++i)
1098 vector<bool> corrected = cleaner.clean_beam(scan_z.row(i), scan_w.row(i), scan_v.row(i),SD2D[0].row(i), SDZDR2D[0].row(i), scan_z, scan_w, scan_v, SD2D[0],i);
1100 for (
unsigned ib = 0; ib < beam_size; ++ib)
1104 scan_z(i, ib) = new_val;
1119 std::string Z_Quantity;
1123 double Cleaner::trap(
double x1,
double x2,
double x3,
double x4,
double val,
double x5)
const
1125 if((val<=x3&&val>=x2))
return 1.;
1126 else if(val<x2&&val>x1)
return val/(x2-x1)-x1/(x2-x1);
1127 else if (val<x4&&val>x3)
return val/(x3-x4)-x4/(x3-x4);
1128 else if(val<=x5)
return 1.;
1133 vector<string> Cleaner::read_matrix_from_txt(
string fin)
const
1139 vector<string> myVector;
1140 ifstream f(fin, ifstream::in);
1144 while(getline(f,line)){
1145 stringstream stream (line);
1146 while( getline(stream, line,
' ')){
1148 if(line[0]==
'\\'){
continue;}
1150 myVector.push_back(line);
double nodata
Value used as 'no data' value.
double undetect
Minimum amount that can be measured.
double offset
Conversion factor.
double gain
Conversion factor.
PolarScan - structure to describe a polarScan containing a matrix of data and conversion factors.
void gdal_init_once()
Initialize the GDAL library when called for the first time; does nothing all other times.
unsigned beam_size
Number of samples in each beam.
unsigned beam_count
Count of beams in this scan.
double cell_size
Size of a beam cell in meters.
std::vector< unsigned char > eval_clean_beam(const Eigen::VectorXd &beam_z, const Eigen::VectorXd &beam_w, const Eigen::VectorXd &beam_v, int i) const
Funzione per ripulire raggio.Utilizza (sigmaV, V) Analoga a clean_beam(const Eigen::VectorXd& beam_z,...
double trap(double x1, double x2, double x3, double x4, double val, double x5=-9999.) const
std::vector< bool > clean_beam(const Eigen::VectorXd &beam_z, const Eigen::VectorXd &beam_w, const Eigen::VectorXd &beam_v, int i) const
Funzione per ripulire raggio.Utilizza (sigmaV, V)
const double Z_missing
Valore dato mancante DBZH.
Struttura per cleaner dati grezzi sulla base dei valori di V, W e la deviazione standard di Z.
Struttura che contiene mappa per caricamento dati.