Elaboradar  0.1
 Tutto Classi Namespace File Funzioni Variabili Tipi enumerati (enum) Gruppi
anaprop.cpp
1 #include "anaprop.h"
2 #include "radarelab/algo/dbz.h"
3 
4 // Soglie algoritmi
5 #define MAX_DIF_OR 30. /* differenzio limiti controllo anap */
6 #define MIN_VALUE_OR -10. /* a seconda che sia alla prima o success.*/
7 #define MAX_DIF_NEXT_OR 15. /*/ elevazione */
8 #define MIN_VALUE_NEXT_OR 0.
9 #define THR_CONT_ANAP 1 /* limite in numero occorrenze anaprop sul raggio dopo i 30 km per non togliere =*/
10 
11 namespace radarelab {
12 namespace algo {
13 
14 using namespace std;
15 
16 namespace anaprop {
17 
18 GridStats::GridStats()
19 {
20 }
21 
22 GridStats::~GridStats()
23 {
24  if (stat_anap) delete[] stat_anap;
25  if (stat_tot) delete[] stat_tot;
26  if (stat_bloc) delete[] stat_bloc;
27  if (stat_elev) delete[] stat_elev;
28 }
29 
30 void GridStats::init(const Volume<double>& volume)
31 {
32  size_az = volume[0].beam_count / step_stat_az + 1;
33  size_beam = volume[0].beam_size / step_stat_range + 1;
34  stat_anap = new unsigned[size_az * size_beam];
35  stat_tot = new unsigned[size_az * size_beam];
36  stat_bloc = new unsigned[size_az * size_beam];
37  stat_elev = new unsigned[size_az * size_beam];
38 
39  for (unsigned i = 0; i < size_az * size_beam; ++i)
40  stat_anap[i] = stat_tot[i] = stat_bloc[i] = stat_elev[i] = 0;
41 }
42 
43 }
44 
45 template<class T>
46 Anaprop<T>::Anaprop(const Volume<T>& volume)
47  : elev_fin(volume), dato_corrotto(volume.beam_count, volume.max_beam_size(), 0),
48  quota(volume.beam_count, volume.max_beam_size(), 0)
49 {
50  logging_category = log4c_category_get("radar.anaprop");
51  LOG_WARN("Anaprop init");
52  grid_stats.init(volume);
53 }
54 
55 template<class T>
56 void Anaprop<T>::init_elev_fin_static(const Volume<T>& volume, const PolarScan<unsigned char>& first_level_static)
57 {
58  for(unsigned i=0; i < volume[0].beam_count; ++i)
59  for(unsigned k=0; k < volume[0].beam_size; ++k)
60  //---assegno el_inf a mappa statica
61  elev_fin(i, k) = first_level_static(i, k);
62 }
63 
64 template<class T>
65 void Anaprop<T>::remove(
66  Volume<T>& volume,
67  PolarScan<unsigned char>& beam_blocking,
68  const PolarScan<unsigned char>& first_level,
69  const PolarScan<unsigned char>& first_level_static,
70  const Volume<double>& SD)
71 {
72  const double fondo_scala = volume[0].undetect;
73  unsigned NUM_AZ_X_PPI=volume[0].beam_count;
74 
75  //for(unsigned i=200; i<201; i++)
76  for (unsigned i = 0; i < volume.beam_count; ++i)
77  {
78  //------------assegno le soglie per anaprop : se sono oltre 60 km e se la differenza tra il bin sotto il base e quello sopra <10 non applico test (cambio i limiti per renderli inefficaci)
79  /* differenza massima tra le due elevazioni successive perchè non sia clutter e valore minimo a quella superiore pe il primo e per i successivi (NEXT) bins*/
80 
81  bool do_test_AP = true;
82  double MAX_DIF=MAX_DIF_OR;
83  double MAX_DIF_NEXT=MAX_DIF_NEXT_OR;
84  double MIN_VALUE=MIN_VALUE_OR;
85  double MIN_VALUE_NEXT=MIN_VALUE_NEXT_OR;
86 
87  double MIN_VALUE_USED, MAX_DIF_USED;
88 
89  bool flag_anap = false;
90  unsigned cont_anap=0;// aggiunto per risolvere problema di uso con preci shallow
91  unsigned count_first_elev=0;
92 
93  //for(unsigned k=150 ; k<volume[0].beam_size; k++)
94  for(unsigned k=0; k<volume[0].beam_size; k++)
95  {
96 
97  char buf1 [256], buf2[256], buf3[256];
98  unsigned count_low =0;
99  unsigned count_high=0;
100  buf3[0]='\0';
101  //------------- incremento statistica tot ------------------
102  grid_stats.incr_tot(i, k);
103  // ------------assegno l'elevazione el_inf a first_level e elev_fin a el_inf---------
104  //int loc_el_inf = first_level(i, k);
105  int loc_el_inf = first_level(i, k) < volume.size() ? first_level(i,k) : volume.size()-1;
106 // LOG_WARN(" Decremento %d %d %d %d ",i,k,loc_el_inf, volume.size() );
107  while ( k >= volume[loc_el_inf].beam_size)
108  {
109 // LOG_INFO("Decremento el_inf per k fuori range (i,k,beam_size,el_inf_dec) (%d,%d,%d,%d)",i,k,volume[loc_el_inf].beam_size,loc_el_inf-1);
110  loc_el_inf--;
111  if (loc_el_inf < 0) throw std::runtime_error("loc_el_inf < 0");
112  }
113  while (loc_el_inf > 0 && SD[loc_el_inf-1].get(i,k) < conf_texture_threshold && SD[loc_el_inf-1].get(i,k)>=0.01 && volume[loc_el_inf-1].get(i,k) > volume[loc_el_inf].get(i,k)){
114 // LOG_WARN("Decremento el_inf Sotto esiste qualcosa %2d %3d %3d %6.2f %6.2f %6.2f",loc_el_inf, i, k , SD[loc_el_inf-1].get(i,k),volume[loc_el_inf-1].get(i,k),volume[loc_el_inf].get(i,k));
115  loc_el_inf--;
116  }
117  const unsigned el_inf = loc_el_inf;
118  if (el_inf == 0) count_first_elev++;
119  if (do_quality)
120  elev_fin(i, k) = el_inf;
121 
122  // ------------assegno a el_up il successivo di el_inf e se >=NEL metto bin_high=fondo_scala
123  //const unsigned el_up = el_inf +1;
124  unsigned el_up = el_inf +1;
125 
126  // ------------assegno bin_low e bin_high
127 
128  float bin_low = volume[el_inf].get(i, k);
129  float bin_high;
130  if (el_up >= volume.size() || k >= volume[el_up].beam_size){
131  el_up=el_inf;
132  bin_high = volume[el_up].get(i, k);
133  //bin_high=fondo_scala;
134  } else{
135  bin_high = volume[el_up].get(i, k);
136  }
137  // LOG_WARN(" utilizzo %d %d %d %d ",i,k,el_inf, el_up );
138 
139  //----------questo serviva per evitare di tagliare la precipitazione shallow ma si dovrebbe trovare un metodo migliore p.es. v. prove su soglia
140  if(bin_high == fondo_scala && SD[el_inf].get(i,k)<= conf_texture_threshold && SD[el_inf].get(i,k) > 0.01) //-----------ANNULLO EFFETTO TEST ANAP
141  {
142  do_test_AP=false;
143  MAX_DIF_NEXT=DBZ::BYTEtoDB(255);
144  MAX_DIF=DBZ::BYTEtoDB(255);
145  MIN_VALUE=fondo_scala;
146  MIN_VALUE_NEXT= fondo_scala;
147  }
148  else
149  {
150  do_test_AP=true;
151  MAX_DIF=MAX_DIF_OR;
152  MAX_DIF_NEXT=MAX_DIF_NEXT_OR;
153  MIN_VALUE=MIN_VALUE_OR;
154  MIN_VALUE_NEXT=MIN_VALUE_NEXT_OR;
155  }
156  // ------------separo i diversi casi x analisi anaprop: ho dati sia al livello base che sopra o no e ho trovato anaprop in precedenza sul raggio o no
157  bool test_an;
158  if (cont_anap> THR_CONT_ANAP ||count_first_elev < 80 )
159  test_an=(bin_low > fondo_scala && bin_high >= fondo_scala );
160  else
161  test_an=(bin_low > fondo_scala && bin_high > fondo_scala );
162  sprintf(buf1, "b@(%3d,%3d) - el_inf %2d - el_up %2d -low %6.2f - up %6.2f - cont %3d %1d %1d %6.2f %6.2f %6.2f %6.2f -- %6.2f %6.2f ",i,k,el_inf,el_up,bin_low,bin_high, cont_anap,test_an, flag_anap, MAX_DIF, MIN_VALUE, MAX_DIF_NEXT, MIN_VALUE_NEXT, SD[el_inf].get(i,k),SD[el_up].get((i),k));
163 
164  double loc_conf_text_thre ;
165  if (k <= 20 ) loc_conf_text_thre = 20. ;
166  else loc_conf_text_thre= conf_texture_threshold ;
167 
168  //------------------ se ho qualcosa sia al livello base che sopra allora effettuo il confronto-----------------
169  if(test_an )
170  {
171  //------------------ se ho trovato anap prima nel raggio cambio le soglie le abbasso)-----------------
172  if(flag_anap)
173  {
174  //-----------caso di propagazione anomala presente nella cella precedente ---------
175  MAX_DIF_USED = MAX_DIF_NEXT;
176  MIN_VALUE_USED = MIN_VALUE_NEXT;
177  }
178  else
179  {
180  //-----------caso di propagazione anomala non presente nella cella precedente ---------
181  MAX_DIF_USED = MAX_DIF;
182  MIN_VALUE_USED = MIN_VALUE;
183  }
184 
185  if( do_test_AP &&
186  (
187  (
188  bin_low-bin_high >= MAX_DIF_USED
189  ||
190  (
191  bin_high <= MIN_VALUE_USED
192  &&
193  bin_low > MIN_VALUE + 5
194  )
195  )
196  ||
197  (
198  SD[el_inf].get(i,k) > loc_conf_text_thre
199  // &&
200 // (
201 // SD[el_inf].get((i+1)%NUM_AZ_X_PPI,k) < 3
202 // ||
203 // SD[el_inf].get((i-1+NUM_AZ_X_PPI)%NUM_AZ_X_PPI,k) < 3
204 // )
205  )
206  )
207  )
208  {
209  //--------ricopio valore a el_up su tutte elev inferiori--------------
210  // double loc_conf_text_thre ;
211  // if (k <= 100 ) loc_conf_text_thre =10 ;
212  // else loc_conf_text_thre= loc_conf_text_thre ;
213  if(k < volume[el_up].beam_size && SD[el_up].get(i,k) <= loc_conf_text_thre && SD[el_up].get(i,k) >= 0.01){
214  for(unsigned l=0; l<el_up; l++){
215  volume[l].set(i, k, bin_high); //ALTERN
216  }
217  sprintf(buf2, " tex_th %4.1f cl %4d cu %4d - fin %6.2f - TA-AN",loc_conf_text_thre, count_low, count_high, volume[0].get(i,k)) ;
218  } else {
219  for(unsigned l=0; l<el_up; l++){
220  volume[l].set(i, k, fondo_scala); //ALTERN
221  }
222  if (k < volume[el_up].beam_size) sprintf(buf2, " tex_th %4.1f cl %4d cu %4d - fin %6.2f - TA-AN set to fondo_scala",loc_conf_text_thre,count_low, count_high, volume[0].get(i,k)) ;
223  }
224  //---------assegno l'indicatore di presenza anap nel raggio e incremento statistica anaprop, assegno matrici che memorizzano anaprop e elevazione_finale e azzero beam blocking perchè ho cambiato elevazione
225  flag_anap = true;
226  cont_anap=cont_anap+1;
227  grid_stats.incr_anap(i, k);
228  if (do_quality)
229  {
230  dato_corrotto(i, k)=ANAP_YES;/*matrice risultato test: propagazione anomala*/
231  elev_fin(i, k) = el_up;
232  }
233  if (el_up > first_level_static(i, k)) grid_stats.incr_elev(i, k);//incremento la statistica cambio elevazione
234  if (do_beamblocking)
235  beam_blocking(i, k)=0;
236  }
237 
238  //-----non c'è propagazione anomala:ricopio su tutte e elevazioni il valore di el_inf e correggo il beam blocking, incremento la statistica beam_blocking, assegno matrice anaprop a 0 nel punto , assegno a 0 indicatore anap nel raggio, assegno elevazione finale e incremento statisica cambio elevazione se el_inf > first_level_static(i, k)-----------
239  else
240  {
241  for (unsigned ii=0; ii<7; ii++){
242  int iaz_prox = (i + ii - 3 + volume.beam_count) % volume.beam_count;
243  if( SD[el_inf].get(iaz_prox,k) < loc_conf_text_thre && SD[el_inf].get(iaz_prox,k) > 0.01) count_low++;
244  if( k < SD[el_up].beam_size && SD[el_up].get(iaz_prox,k) < 1.7*loc_conf_text_thre&& SD[el_up].get(iaz_prox,k) > 0.01) count_high++;
245  }
246  if ( !(SD[el_inf].get(i,k) < 1.3*loc_conf_text_thre && SD[el_inf].get(i,k) > 0.01 && count_low >=5)) {
247  if ( k >= SD[el_up].beam_size || !(SD[el_up].get(i,k) < 1.7* loc_conf_text_thre && SD[el_up].get(i,k) > 0.01 && count_high >=3)){
248  bin_low = fondo_scala;
249  if ( k < SD[el_up].beam_size ) sprintf(buf2, " tex_th %4.1f cl %4d cu %4d - fin %6.2f - TA-NO_AN low = fondo",loc_conf_text_thre,count_low, count_high, volume[0].get(i,k)) ;
250  else sprintf(buf2, " tex_th %4.1f cl %4d cu %4d - fin %6.2f - TA-NO_AN low = fondo",loc_conf_text_thre,count_low, count_high, volume[0].get(i,k)) ;
251  }
252  else {
253  bin_low=bin_high;
254  sprintf(buf2, " tex_th %4.1f cl %4d cu %4d - fin %6.2f - TA-NO_AN low = high",loc_conf_text_thre,count_low, count_high, volume[0].get(i,k)) ;
255  if (do_quality)
256  {
257  elev_fin(i, k) = el_up;
258  }
259  if (el_up > first_level_static(i, k)) grid_stats.incr_elev(i, k);//incremento la statistica cambio elevazione
260  if (do_beamblocking)
261  beam_blocking(i, k)=0;
262  }
263  }
264  else {
265  if (do_beamblocking && do_bloccorr)
266  {
267  bin_low = DBZ::beam_blocking_correction(bin_low, beam_blocking(i, k));
268  grid_stats.incr_bloc(i, k, beam_blocking(i, k));
269  }
270  }
271  for(unsigned l=0; l<=el_inf; l++)
272  volume[l].set(i, k, bin_low);
273 sprintf(buf3, " - fin %6.2f - TA-NO_AN",volume[0].get(i,k)) ;
274 
275  if (do_quality)
276  {
277  dato_corrotto(i, k)=ANAP_OK;
278  elev_fin(i, k) = el_inf;
279  }
280  if (el_inf > first_level_static(i, k)) grid_stats.incr_elev(i, k);//incremento la statistica cambio elevazione
281  flag_anap = false;
282  }
283  }/* test_anap */
284  //----------------se al livello base non ho dato riempio con i valori di el_up tutte le elevazioni sotto (ricostruisco il volume) e assegno beam_blocking 0
285  else if (bin_low < fondo_scala)
286  {
287  for(unsigned l=0; l<el_up; l++)
288  {
289  if (volume[l].beam_size > k && volume[el_up].beam_size > k)
290  volume[l].set(i, k, bin_high);
291  else if (volume[l].beam_size > k )
292  volume[l].set(i, k, fondo_scala);
293 
294  }
295 // LOG_WARN("b@(%3d,%3d) - el_inf %2d - el_up %2d -low %6.2f - up %6.2f fin %6.2f- cont %3d %1d %1d %6.2f %6.2f %6.2f %6.2f -- %6.2f %6.2f %6.2f NO_TA-low <fondo",i,k,el_inf,el_up,bin_low,bin_high, volume[0].get(i,k),cont_anap,test_an, flag_anap, MAX_DIF, MIN_VALUE, MAX_DIF_NEXT, MIN_VALUE_NEXT, SD[el_inf].get(i,k),SD[el_inf].get((i+1)%NUM_AZ_X_PPI,k) ,SD[el_inf].get((i-1+NUM_AZ_X_PPI)%NUM_AZ_X_PPI,k));
296  //----------------controlli su bin_high nel caso in cui bin_low sia un no data per assegnare matrice anap (dato_corrotto(i, k))
297  if (do_quality)
298  {
299  if (bin_high<fondo_scala) dato_corrotto(i, k)=ANAP_NODAT;/*manca dato sotto e sopra*/
300  bool test_an1;
301  if (cont_anap< THR_CONT_ANAP )
302  test_an1=(bin_high>=fondo_scala); //modificato per contemplare > o >=
303  else
304  test_an1=(bin_high>fondo_scala);
305 
306  if (test_an1) dato_corrotto(i, k)=ANAP_NOCONTROL;
307  if (bin_high==fondo_scala) dato_corrotto(i, k)=ANAP_OK;/*non piove (oppure sono sopra livello preci...)*/
308  }
309 
310  if (do_beamblocking)
311  beam_blocking(i, k)=0;
312  }
313 
314  //----------------se bin_low == fondo_scala riempio matrice volume.vol_pol con dato a el_inf (mi resta il dubbio di quest'if se seve o basti un else ) azzero matrice anap (dato ok)
315  else if (bin_low == fondo_scala || bin_high <= fondo_scala)/* quel che resta da (bin_low > fondo_scala && bin_high >= fondo_scala) e (bin_low < fondo_scala) ; messo =per bin_high*/
316 
317  {
318  unsigned count =0;
319  for (unsigned ii=0; ii<7; ii++){
320  int iaz = (i + ii - 3 + volume.beam_count) % volume.beam_count;
321  if( SD[el_inf].get(iaz,k) < loc_conf_text_thre && SD[el_inf].get(iaz,k) > 0.01) count++;
322  }
323  if ( !(SD[el_inf].get(i,k) < loc_conf_text_thre && SD[el_inf].get(i,k) >0.01 && count >=5 ))
324  bin_low = fondo_scala;
325 
326  for(unsigned l=0; l<=el_inf; l++)//riempio con i valori di el_inf tutte le elevazioni sotto (ricostruisco il volume)
327  {
328  if (volume[l].beam_size > k)
329  volume[l].set(i, k, bin_low);
330  }
331 sprintf(buf2, " tex_th %4.1f cl %4d cu %4d - fin %6.2f - NO_TA low = fondo",loc_conf_text_thre,count_low, count_high, volume[0].get(i,k)) ;
332 
333  if (do_quality)
334  {
335  dato_corrotto(i, k)=ANAP_OK; // dubbio
336  elev_fin(i, k) = el_inf;
337  }
338 
339  if (el_inf > first_level_static(i, k)) grid_stats.incr_elev(i, k);
340  }
341  /*-----------------------------------------------------------fine di tutti gli if-----------*/
342  //-----finiti tutti i controlli assegno le varibili di qualita definitive: elevazione, quota calcolata sull'elevazione reale con propagazione standard , e quota relativa al suolo calcolata con elevazione nominale e propagazione da radiosondaggio.
343 
344  if (do_quality)
345  quota(i, k)=(unsigned short)PolarScanBase::sample_height(
346  elev_fin.elevation_at_elev_preci(i, k),
347  (k + 0.5) * volume[0].cell_size);
348 
349 // if (k < 20) LOG_WARN("%s %s %s",buf1,buf2, buf3);
350 // if (k == 20) LOG_WARN("b@(%3d,%3d) ------------------------",i,k);
351  }
352  } // end for over beam_count
353 }
354 
355 template<class T>
356 void Anaprop<T>::remove_without_SD(
357  Volume<T>& volume,
358  PolarScan<unsigned char>& beam_blocking,
359  const PolarScan<unsigned char>& first_level,
360  const PolarScan<unsigned char>& first_level_static,
361  const Volume<double>& SD)
362 {
363  const double fondo_scala = volume[0].undetect;
364 LOG_WARN("Anaprop remove without SD");
365 
366  //for(unsigned i=200; i<201; i++)
367  for(unsigned i=0; i<volume[0].beam_count; i++)
368  {
369  //------------assegno le soglie per anaprop : se sono oltre 60 km e se la differenza tra il bin sotto il base e quello sopra <10 non applico test (cambio i limiti per renderli inefficaci)
370  /* differenza massima tra le due elevazioni successive perchè non sia clutter e valore minimo a quella superiore pe il primo e per i successivi (NEXT) bins*/
371 
372  bool do_test_AP = true;
373  double MAX_DIF=MAX_DIF_OR;
374  double MAX_DIF_NEXT=MAX_DIF_NEXT_OR;
375  double MIN_VALUE=MIN_VALUE_OR;
376  double MIN_VALUE_NEXT=MIN_VALUE_NEXT_OR;
377 
378  double MIN_VALUE_USED, MAX_DIF_USED;
379 
380  bool flag_anap = false;
381  unsigned cont_anap=0;// aggiunto per risolvere problema di uso con preci shallow
382  unsigned count_first_elev=0;
383 
384  //for(unsigned k=150 ; k<volume[0].beam_size; k++)
385  for(unsigned k=0; k<volume[0].beam_size; k++)
386  {
387  //------------- incremento statistica tot ------------------
388  grid_stats.incr_tot(i, k);
389  // ------------assegno l'elevazione el_inf a first_level e elev_fin a el_inf---------
390  int loc_el_inf = first_level(i, k) < volume.size() ? first_level(i,k) : volume.size()-1;
391  while ( k >= volume[loc_el_inf].beam_size)
392  {
393  // LOG_INFO("Decremento el_inf per k fuori range (i,k,beam_size,el_inf_dec) (%d,%d,%d,%d)",i,k,volume[loc_el_inf].beam_size,loc_el_inf-1);
394  loc_el_inf--;
395  if (loc_el_inf < 0) throw std::runtime_error("loc_el_inf < 0");
396  }
397 /* ---------------------------------
398  * PER IL MOMENTO NON BUTTO ANCORA QUESTO CODICE CI DEVO PENSARE
399  */
400  while (loc_el_inf > 0 && SD[loc_el_inf-1].get(i,k) < conf_texture_threshold && SD[loc_el_inf-1].get(i,k)>=0.01 && volume[loc_el_inf-1].get(i,k) > volume[loc_el_inf].get(i,k)){
401  // LOG_WARN("Decremento el_inf Sotto esiste qualcosa %2d %3d %3d %6.2f %6.2f %6.2f",loc_el_inf, i, k , SD[loc_el_inf-1].get(i,k),volume[loc_el_inf-1].get(i,k),volume[loc_el_inf].get(i,k));
402  loc_el_inf--;
403  }
404 /* */
405  const unsigned el_inf = loc_el_inf;
406  if (el_inf == 0) count_first_elev++;
407  if (do_quality)
408  elev_fin(i, k) = el_inf;
409 
410  // ------------assegno a el_up il successivo di el_inf e se >=NEL metto bin_high=fondo_scala
411  //const unsigned el_up = el_inf +1;
412  unsigned el_up = el_inf +1;
413 
414  // ------------assegno bin_low e bin_high
415 
416  float bin_low = volume[el_inf].get(i, k);
417  float bin_high;
418  if (el_up >= volume.size() || k >= volume[el_up].beam_size){
419  el_up=el_inf;
420  bin_high = volume[el_up].get(i, k);
421  //bin_high=fondo_scala;
422  } else{
423  bin_high = volume[el_up].get(i, k);
424  }
425 
426  //----------questo serviva per evitare di tagliare la precipitazione shallow ma si dovrebbe trovare un metodo migliore p.es. v. prove su soglia
427  if(bin_high == fondo_scala && SD[el_inf].get(i,k)<= conf_texture_threshold && SD[el_inf].get(i,k) > 0.01) //-----------ANNULLO EFFETTO TEST ANAP
428  {
429  do_test_AP=false;
430  MAX_DIF_NEXT=DBZ::BYTEtoDB(255);
431  MAX_DIF=DBZ::BYTEtoDB(255);
432  MIN_VALUE=fondo_scala;
433  MIN_VALUE_NEXT= fondo_scala;
434  }
435  else
436  {
437  do_test_AP=true;
438  MAX_DIF=MAX_DIF_OR;
439  MAX_DIF_NEXT=MAX_DIF_NEXT_OR;
440  MIN_VALUE=MIN_VALUE_OR;
441  MIN_VALUE_NEXT=MIN_VALUE_NEXT_OR;
442  }
443  // ------------separo i diversi casi x analisi anaprop: ho dati sia al livello base che sopra o no e ho trovato anaprop in precedenza sul raggio o no
444  bool test_an;
445  if (cont_anap> THR_CONT_ANAP ||count_first_elev < 80 )
446  test_an=(bin_low > fondo_scala && bin_high >= fondo_scala );
447  else
448  test_an=(bin_low > fondo_scala && bin_high > fondo_scala );
449  // LOG_WARN("b@(%3d,%3d) - el_inf %2d - el_up %2d -low %6.2f - up %6.2f - cont %3d %1d %1d %6.2f %6.2f %6.2f %6.2f ",i,k,el_inf,el_up,bin_low,bin_high, cont_anap,test_an, flag_anap, MAX_DIF, MIN_VALUE, MAX_DIF_NEXT, MIN_VALUE_NEXT);
450  //------------------ se ho qualcosa sia al livello base che sopra allora effettuo il confronto-----------------
451  if(test_an )
452  {
453  //------------------ se ho trovato anap prima nel raggio cambio le soglie le abbasso)-----------------
454  if(flag_anap)
455  {
456  //-----------caso di propagazione anomala presente nella cella precedente ---------
457  MAX_DIF_USED = MAX_DIF_NEXT;
458  MIN_VALUE_USED = MIN_VALUE_NEXT;
459  }
460  else
461  {
462  //-----------caso di propagazione anomala non presente nella cella precedente ---------
463  MAX_DIF_USED = MAX_DIF;
464  MIN_VALUE_USED = MIN_VALUE;
465  }
466 
467  if( do_test_AP &&
468  ( bin_low-bin_high >= MAX_DIF_USED
469  ||
470  (
471  bin_high <= MIN_VALUE_USED
472  &&
473  bin_low > MIN_VALUE + 5
474  )
475  )
476  )
477  {
478  //--------ricopio valore a el_up su tutte elev inferiori--------------
479  if(k < volume[el_up].beam_size ){
480  for(unsigned l=0; l<el_up; l++){
481  volume[l].set(i, k, bin_high); //ALTERN
482  }
483  // LOG_WARN("b@(%3d,%3d) - el_inf %2d - el_up %2d -low %6.2f - up %6.2f fin %6.2f- cont %3d %1d %1d %6.2f %6.2f %6.2f %6.2f TA-AN",i,k,el_inf,el_up,bin_low,bin_high, volume[0].get(i,k),cont_anap,test_an, flag_anap, MAX_DIF, MIN_VALUE, MAX_DIF_NEXT, MIN_VALUE_NEXT ) ;
484  } else {
485  for(unsigned l=0; l<el_up; l++){
486  volume[l].set(i, k, fondo_scala); //ALTERN
487  }
488  // if (k < volume[el_up].beam_size) LOG_WARN("b@(%3d,%3d) - el_inf %2d - el_up %2d -low %6.2f - up %6.2f fin %6.2f- cont %3d %1d %1d %6.2f %6.2f %6.2f %6.2f TA-AN set to fondo_scala",i,k,el_inf,el_up,bin_low,bin_high, volume[0].get(i,k),cont_anap,test_an, flag_anap, MAX_DIF, MIN_VALUE, MAX_DIF_NEXT, MIN_VALUE_NEXT);
489  }
490  //---------assegno l'indicatore di presenza anap nel raggio e incremento statistica anaprop, assegno matrici che memorizzano anaprop e elevazione_finale e azzero beam blocking perchè ho cambiato elevazione
491  flag_anap = true;
492  cont_anap=cont_anap+1;
493  grid_stats.incr_anap(i, k);
494  if (do_quality)
495  {
496  dato_corrotto(i, k)=ANAP_YES;/*matrice risultato test: propagazione anomala*/
497  elev_fin(i, k) = el_up;
498  }
499  if (el_up > first_level_static(i, k)) grid_stats.incr_elev(i, k);//incremento la statistica cambio elevazione
500  if (do_beamblocking)
501  beam_blocking(i, k)=0;
502  }
503 
504  //-----non c'è propagazione anomala:ricopio su tutte e elevazioni il valore di el_inf e correggo il beam blocking, incremento la statistica beam_blocking, assegno matrice anaprop a 0 nel punto , assegno a 0 indicatore anap nel raggio, assegno elevazione finale e incremento statisica cambio elevazione se el_inf > first_level_static(i, k)-----------
505  else
506  {
507  if (do_beamblocking && do_bloccorr)
508  {
509  bin_low = DBZ::beam_blocking_correction(bin_low, beam_blocking(i, k));
510  grid_stats.incr_bloc(i, k, beam_blocking(i, k));
511  }
512 
513  for(unsigned l=0; l<=el_inf; l++)
514  volume[l].set(i, k, bin_low);
515 // LOG_WARN("b@(%3d,%3d) - el_inf %2d - el_up %2d -low %6.2f - up %6.2f fin %6.2f- cont %3d %1d %1d %6.2f %6.2f %6.2f %6.2f TA-NO_AN",i,k,el_inf,el_up,bin_low,bin_high, volume[0].get(i,k),cont_anap,test_an, flag_anap, MAX_DIF, MIN_VALUE, MAX_DIF_NEXT, MIN_VALUE_NEXT );
516 
517  if (do_quality)
518  {
519  dato_corrotto(i, k)=ANAP_OK;
520  elev_fin(i, k) = el_inf;
521  }
522  if (el_inf > first_level_static(i, k)) grid_stats.incr_elev(i, k);//incremento la statistica cambio elevazione
523  flag_anap = false;
524  }
525  }/* test_anap */
526  //----------------se al livello base non ho dato riempio con i valori di el_up tutte le elevazioni sotto (ricostruisco il volume) e assegno beam_blocking 0
527  else if (bin_low < fondo_scala)
528  {
529  for(unsigned l=0; l<el_up; l++)
530  {
531  if (volume[l].beam_size > k && volume[el_up].beam_size > k)
532  volume[l].set(i, k, bin_high);
533  else if (volume[l].beam_size > k )
534  volume[l].set(i, k, fondo_scala);
535 
536  }
537 // LOG_WARN("b@(%3d,%3d) - el_inf %2d - el_up %2d -low %6.2f - up %6.2f fin %6.2f- cont %3d %1d %1d %6.2f %6.2f %6.2f %6.2f NO_TA-low <fondo",i,k,el_inf,el_up,bin_low,bin_high, volume[0].get(i,k),cont_anap,test_an, flag_anap, MAX_DIF, MIN_VALUE, MAX_DIF_NEXT, MIN_VALUE_NEXT);
538  //----------------controlli su bin_high nel caso in cui bin_low sia un no data per assegnare matrice anap (dato_corrotto(i, k))
539  if (do_quality)
540  {
541  if (bin_high<fondo_scala) dato_corrotto(i, k)=ANAP_NODAT;/*manca dato sotto e sopra*/
542  bool test_an1;
543  if (cont_anap< THR_CONT_ANAP )
544  test_an1=(bin_high>=fondo_scala); //modificato per contemplare > o >=
545  else
546  test_an1=(bin_high>fondo_scala);
547 
548  if (test_an1) dato_corrotto(i, k)=ANAP_NOCONTROL;
549  if (bin_high==fondo_scala) dato_corrotto(i, k)=ANAP_OK;/*non piove (oppure sono sopra livello preci...)*/
550  }
551 
552  if (do_beamblocking)
553  beam_blocking(i, k)=0;
554  }
555 
556  //----------------se bin_low == fondo_scala riempio matrice volume.vol_pol con dato a el_inf (mi resta il dubbio di quest'if se seve o basti un else ) azzero matrice anap (dato ok)
557  else if (bin_low == fondo_scala || bin_high <= fondo_scala)/* quel che resta da (bin_low > fondo_scala && bin_high >= fondo_scala) e (bin_low < fondo_scala) ; messo =per bin_high*/
558 
559  {
560  unsigned count =0;
561  for (unsigned ii=0; ii<7; ii++){
562  int iaz = (i + ii - 3 + volume.beam_count) % volume.beam_count;
563  if( SD[el_inf].get(iaz,k) < conf_texture_threshold && SD[el_inf].get(iaz,k) > 0.01) count++;
564  }
565  if ( !(SD[el_inf].get(i,k) < conf_texture_threshold && SD[el_inf].get(i,k) >0.01 && count >=5 ))
566  bin_low = fondo_scala;
567 
568  for(unsigned l=0; l<=el_inf; l++)//riempio con i valori di el_inf tutte le elevazioni sotto (ricostruisco il volume)
569  {
570  if (volume[l].beam_size > k)
571  volume[l].set(i, k, bin_low);
572  }
573 // LOG_WARN("b@(%3d,%3d) - el_inf %2d - el_up %2d -low %6.2f - up %6.2f fin %6.2f- cont %3d %1d %1d %6.2f %6.2f %6.2f %6.2f NO_TA-low ==fondo",i,k,el_inf,el_up,bin_low,bin_high, volume[0].get(i,k),cont_anap,test_an, flag_anap, MAX_DIF, MIN_VALUE, MAX_DIF_NEXT, MIN_VALUE_NEXT);
574 
575  if (do_quality)
576  {
577  dato_corrotto(i, k)=ANAP_OK; // dubbio
578  elev_fin(i, k) = el_inf;
579  }
580 
581  if (el_inf > first_level_static(i, k)) grid_stats.incr_elev(i, k);
582  }
583  /*-----------------------------------------------------------fine di tutti gli if-----------*/
584  //-----finiti tutti i controlli assegno le varibili di qualita definitive: elevazione, quota calcolata sull'elevazione reale con propagazione standard , e quota relativa al suolo calcolata con elevazione nominale e propagazione da radiosondaggio.
585 
586  if (do_quality)
587  quota(i, k)=(unsigned short)PolarScanBase::sample_height(
588  elev_fin.elevation_at_elev_preci(i, k),
589  (k + 0.5) * volume[0].cell_size);
590  }
591  } // end for over beam_count
592 }
593 
594 // Explicit template instantiation for T=double
595 template class Anaprop<double>;
596 
597 }
598 }