Elaboradar  0.1
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_nodata = volume[el_inf].nodata;
129  float bin_low = volume[el_inf].get(i, k);
130  float bin_high;
131  if (el_up >= volume.size() || k >= volume[el_up].beam_size){
132  el_up=el_inf;
133  bin_high = volume[el_up].get(i, k);
134  //bin_high=fondo_scala;
135  } else{
136  bin_high = volume[el_up].get(i, k);
137  }
138  // LOG_WARN(" utilizzo %d %d %d %d ",i,k,el_inf, el_up );
139  //float bin_high_nodata = volume[el_up].nodata;
140 
141  //----------questo serviva per evitare di tagliare la precipitazione shallow ma si dovrebbe trovare un metodo migliore p.es. v. prove su soglia
142  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
143  {
144  do_test_AP=false;
145  MAX_DIF_NEXT=DBZ::BYTEtoDB(255);
146  MAX_DIF=DBZ::BYTEtoDB(255);
147  MIN_VALUE=fondo_scala;
148  MIN_VALUE_NEXT= fondo_scala;
149  }
150  else
151  {
152  do_test_AP=true;
153  MAX_DIF=MAX_DIF_OR;
154  MAX_DIF_NEXT=MAX_DIF_NEXT_OR;
155  MIN_VALUE=MIN_VALUE_OR;
156  MIN_VALUE_NEXT=MIN_VALUE_NEXT_OR;
157  }
158  // ------------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
159  bool test_an;
160  if (cont_anap> THR_CONT_ANAP ||count_first_elev < 80 )
161  test_an=(bin_low > fondo_scala && bin_high >= fondo_scala );
162  else
163  test_an=(bin_low > fondo_scala && bin_high > fondo_scala );
164  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));
165 
166  double loc_conf_text_thre ;
167  if (k <= 20 ) loc_conf_text_thre = 20. ;
168  else loc_conf_text_thre= conf_texture_threshold ;
169 
170  //------------------ se ho qualcosa sia al livello base che sopra allora effettuo il confronto-----------------
171  if(test_an )
172  {
173  //------------------ se ho trovato anap prima nel raggio cambio le soglie le abbasso)-----------------
174  if(flag_anap)
175  {
176  //-----------caso di propagazione anomala presente nella cella precedente ---------
177  MAX_DIF_USED = MAX_DIF_NEXT;
178  MIN_VALUE_USED = MIN_VALUE_NEXT;
179  }
180  else
181  {
182  //-----------caso di propagazione anomala non presente nella cella precedente ---------
183  MAX_DIF_USED = MAX_DIF;
184  MIN_VALUE_USED = MIN_VALUE;
185  }
186 
187  if( do_test_AP &&
188  (
189  (
190  bin_low-bin_high >= MAX_DIF_USED
191  ||
192  (
193  bin_high <= MIN_VALUE_USED
194  &&
195  bin_low > MIN_VALUE + 5
196  )
197  )
198  ||
199  (
200  SD[el_inf].get(i,k) > loc_conf_text_thre
201  // &&
202 // (
203 // SD[el_inf].get((i+1)%NUM_AZ_X_PPI,k) < 3
204 // ||
205 // SD[el_inf].get((i-1+NUM_AZ_X_PPI)%NUM_AZ_X_PPI,k) < 3
206 // )
207  )
208  )
209  )
210  {
211  //--------ricopio valore a el_up su tutte elev inferiori--------------
212  // double loc_conf_text_thre ;
213  // if (k <= 100 ) loc_conf_text_thre =10 ;
214  // else loc_conf_text_thre= loc_conf_text_thre ;
215  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){
216  for(unsigned l=0; l<el_up; l++){
217  volume[l].set(i, k, bin_high); //ALTERN
218  }
219  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)) ;
220  } else {
221  for(unsigned l=0; l<el_up; l++){
222  volume[l].set(i, k, fondo_scala); //ALTERN
223  }
224  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)) ;
225  }
226  //---------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
227  flag_anap = true;
228  cont_anap=cont_anap+1;
229  grid_stats.incr_anap(i, k);
230  if (do_quality)
231  {
232  dato_corrotto(i, k)=ANAP_YES;/*matrice risultato test: propagazione anomala*/
233  elev_fin(i, k) = el_up;
234  }
235  if (el_up > first_level_static(i, k)) grid_stats.incr_elev(i, k);//incremento la statistica cambio elevazione
236  if (do_beamblocking)
237  beam_blocking(i, k)=0;
238  }
239 
240  //-----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)-----------
241  else
242  {
243  for (unsigned ii=0; ii<7; ii++){
244  int iaz_prox = (i + ii - 3 + volume.beam_count) % volume.beam_count;
245  if( SD[el_inf].get(iaz_prox,k) < loc_conf_text_thre && SD[el_inf].get(iaz_prox,k) > 0.01) count_low++;
246  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++;
247  }
248  if ( !(SD[el_inf].get(i,k) < 1.3*loc_conf_text_thre && SD[el_inf].get(i,k) > 0.01 && count_low >=5)) {
249  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)){
250  bin_low = fondo_scala;
251  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)) ;
252  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)) ;
253  }
254  else {
255  bin_low=bin_high;
256  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)) ;
257  if (do_quality)
258  {
259  elev_fin(i, k) = el_up;
260  }
261  if (el_up > first_level_static(i, k)) grid_stats.incr_elev(i, k);//incremento la statistica cambio elevazione
262  if (do_beamblocking)
263  beam_blocking(i, k)=0;
264  }
265  }
266  else {
267  if (do_beamblocking && do_bloccorr)
268  {
269  bin_low = DBZ::beam_blocking_correction(bin_low, beam_blocking(i, k));
270  grid_stats.incr_bloc(i, k, beam_blocking(i, k));
271  }
272  }
273  for(unsigned l=0; l<=el_inf; l++)
274  volume[l].set(i, k, bin_low);
275 sprintf(buf3, " - fin %6.2f - TA-NO_AN",volume[0].get(i,k)) ;
276 
277  if (do_quality)
278  {
279  dato_corrotto(i, k)=ANAP_OK;
280  elev_fin(i, k) = el_inf;
281  }
282  if (el_inf > first_level_static(i, k)) grid_stats.incr_elev(i, k);//incremento la statistica cambio elevazione
283  flag_anap = false;
284  }
285  }/* test_anap */
286  //----------------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
287  else if (bin_low < fondo_scala)
288  {
289  for(unsigned l=0; l<el_up; l++)
290  {
291  if (volume[l].beam_size > k && volume[el_up].beam_size > k)
292  volume[l].set(i, k, bin_high);
293  else if (volume[l].beam_size > k )
294  volume[l].set(i, k, fondo_scala);
295 
296  }
297 // 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));
298  //----------------controlli su bin_high nel caso in cui bin_low sia un no data per assegnare matrice anap (dato_corrotto(i, k))
299  if (do_quality)
300  {
301  if (bin_high<fondo_scala) dato_corrotto(i, k)=ANAP_NODAT;/*manca dato sotto e sopra*/
302  bool test_an1;
303  if (cont_anap< THR_CONT_ANAP )
304  test_an1=(bin_high>=fondo_scala); //modificato per contemplare > o >=
305  else
306  test_an1=(bin_high>fondo_scala);
307 
308  if (test_an1) dato_corrotto(i, k)=ANAP_NOCONTROL;
309  if (bin_high==fondo_scala) dato_corrotto(i, k)=ANAP_OK;/*non piove (oppure sono sopra livello preci...)*/
310  }
311 
312  if (do_beamblocking)
313  beam_blocking(i, k)=0;
314  }
315 
316  //----------------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)
317  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*/
318 
319  {
320  unsigned count =0;
321  for (unsigned ii=0; ii<7; ii++){
322  int iaz = (i + ii - 3 + volume.beam_count) % volume.beam_count;
323  if( SD[el_inf].get(iaz,k) < loc_conf_text_thre && SD[el_inf].get(iaz,k) > 0.01) count++;
324  }
325  if ( !(SD[el_inf].get(i,k) < loc_conf_text_thre && SD[el_inf].get(i,k) >0.01 && count >=5 ))
326  bin_low = fondo_scala;
327 
328  for(unsigned l=0; l<=el_inf; l++)//riempio con i valori di el_inf tutte le elevazioni sotto (ricostruisco il volume)
329  {
330  if (volume[l].beam_size > k)
331  volume[l].set(i, k, bin_low);
332  }
333 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)) ;
334 
335  if (do_quality)
336  {
337  dato_corrotto(i, k)=ANAP_OK; // dubbio
338  elev_fin(i, k) = el_inf;
339  }
340 
341  if (el_inf > first_level_static(i, k)) grid_stats.incr_elev(i, k);
342  }
343  /*-----------------------------------------------------------fine di tutti gli if-----------*/
344  //-----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.
345 
346  if (do_quality)
347  quota(i, k)=(unsigned short)PolarScanBase::sample_height(
348  elev_fin.elevation_at_elev_preci(i, k),
349  (k + 0.5) * volume[0].cell_size);
350 
351 // if (k < 20) LOG_WARN("%s %s %s",buf1,buf2, buf3);
352 // if (k == 20) LOG_WARN("b@(%3d,%3d) ------------------------",i,k);
353  }
354  } // end for over beam_count
355 }
356 
357 template<class T>
358 void Anaprop<T>::remove_without_SD(
359  Volume<T>& volume,
360  PolarScan<unsigned char>& beam_blocking,
361  const PolarScan<unsigned char>& first_level,
362  const PolarScan<unsigned char>& first_level_static,
363  const Volume<double>& SD)
364 {
365  const double fondo_scala = volume[0].undetect;
366 LOG_WARN("Anaprop remove without SD");
367 
368  //for(unsigned i=200; i<201; i++)
369  for(unsigned i=0; i<volume[0].beam_count; i++)
370  {
371  //------------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)
372  /* 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*/
373 
374  bool do_test_AP = true;
375  double MAX_DIF=MAX_DIF_OR;
376  double MAX_DIF_NEXT=MAX_DIF_NEXT_OR;
377  double MIN_VALUE=MIN_VALUE_OR;
378  double MIN_VALUE_NEXT=MIN_VALUE_NEXT_OR;
379 
380  double MIN_VALUE_USED, MAX_DIF_USED;
381 
382  bool flag_anap = false;
383  unsigned cont_anap=0;// aggiunto per risolvere problema di uso con preci shallow
384  unsigned count_first_elev=0;
385 
386  //for(unsigned k=150 ; k<volume[0].beam_size; k++)
387  for(unsigned k=0; k<volume[0].beam_size; k++)
388  {
389  //------------- incremento statistica tot ------------------
390  grid_stats.incr_tot(i, k);
391  // ------------assegno l'elevazione el_inf a first_level e elev_fin a el_inf---------
392  int loc_el_inf = first_level(i, k) < volume.size() ? first_level(i,k) : volume.size()-1;
393 
394  while ( k >= volume[loc_el_inf].beam_size)
395  {
396  // 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);
397  loc_el_inf--;
398  if (loc_el_inf < 0) throw std::runtime_error("loc_el_inf < 0");
399  }
400 
401 /* ---------------------------------
402  * PER IL MOMENTO NON BUTTO ANCORA QUESTO CODICE CI DEVO PENSARE
403  */
404  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)){
405  // 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));
406  loc_el_inf--;
407  }
408 /* */
409  const unsigned el_inf = loc_el_inf;
410  if (el_inf == 0) count_first_elev++;
411  if (do_quality)
412  elev_fin(i, k) = el_inf;
413 
414  // ------------assegno a el_up il successivo di el_inf e se >=NEL metto bin_high=fondo_scala
415  //const unsigned el_up = el_inf +1;
416  unsigned el_up = el_inf +1;
417 
418  // ------------assegno bin_low e bin_high
419 
420  float bin_low = volume[el_inf].get(i, k);
421  float bin_high;
422  if (el_up >= volume.size() || k >= volume[el_up].beam_size){
423  el_up=el_inf;
424  //cout<<" i="<<i<<" el_up= "<<el_up<<" volume.beam_size= "<<volume[el_up].beam_size<<endl;
425  bin_high = volume[el_up].get(i, k);
426  //bin_high=fondo_scala;
427  } else{
428  bin_high = volume[el_up].get(i, k);
429  }
430 
431  //----------questo serviva per evitare di tagliare la precipitazione shallow ma si dovrebbe trovare un metodo migliore p.es. v. prove su soglia
432  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
433  {
434  do_test_AP=false;
435  MAX_DIF_NEXT=DBZ::BYTEtoDB(255);
436  MAX_DIF=DBZ::BYTEtoDB(255);
437  MIN_VALUE=fondo_scala;
438  MIN_VALUE_NEXT= fondo_scala;
439  }
440  else
441  {
442  do_test_AP=true;
443  MAX_DIF=MAX_DIF_OR;
444  MAX_DIF_NEXT=MAX_DIF_NEXT_OR;
445  MIN_VALUE=MIN_VALUE_OR;
446  MIN_VALUE_NEXT=MIN_VALUE_NEXT_OR;
447  }
448  // ------------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
449  bool test_an;
450  if (cont_anap> THR_CONT_ANAP ||count_first_elev < 80 )
451  test_an=(bin_low > fondo_scala && bin_high >= fondo_scala );
452  else
453  test_an=(bin_low > fondo_scala && bin_high > fondo_scala );
454  // 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);
455  //------------------ se ho qualcosa sia al livello base che sopra allora effettuo il confronto-----------------
456  if(test_an )
457  {
458  //------------------ se ho trovato anap prima nel raggio cambio le soglie le abbasso)-----------------
459  if(flag_anap)
460  {
461  //-----------caso di propagazione anomala presente nella cella precedente ---------
462  MAX_DIF_USED = MAX_DIF_NEXT;
463  MIN_VALUE_USED = MIN_VALUE_NEXT;
464  }
465  else
466  {
467  //-----------caso di propagazione anomala non presente nella cella precedente ---------
468  MAX_DIF_USED = MAX_DIF;
469  MIN_VALUE_USED = MIN_VALUE;
470  }
471 
472  if( do_test_AP &&
473  ( bin_low-bin_high >= MAX_DIF_USED
474  ||
475  (
476  bin_high <= MIN_VALUE_USED
477  &&
478  bin_low > MIN_VALUE + 5
479  )
480  )
481  )
482  {
483  //--------ricopio valore a el_up su tutte elev inferiori--------------
484  if(k < volume[el_up].beam_size ){
485  for(unsigned l=0; l<el_up; l++){
486  volume[l].set(i, k, bin_high); //ALTERN
487  }
488  // 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 ) ;
489  } else {
490  for(unsigned l=0; l<el_up; l++){
491  volume[l].set(i, k, fondo_scala); //ALTERN
492  }
493  // 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);
494  }
495  //---------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
496  flag_anap = true;
497  cont_anap=cont_anap+1;
498  grid_stats.incr_anap(i, k);
499  if (do_quality)
500  {
501  dato_corrotto(i, k)=ANAP_YES;/*matrice risultato test: propagazione anomala*/
502  elev_fin(i, k) = el_up;
503  }
504  if (el_up > first_level_static(i, k)) grid_stats.incr_elev(i, k);//incremento la statistica cambio elevazione
505  if (do_beamblocking)
506  beam_blocking(i, k)=0;
507  }
508 
509  //-----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)-----------
510  else
511  {
512  if (do_beamblocking && do_bloccorr)
513  {
514  bin_low = DBZ::beam_blocking_correction(bin_low, beam_blocking(i, k));
515  grid_stats.incr_bloc(i, k, beam_blocking(i, k));
516  }
517 
518  for(unsigned l=0; l<=el_inf; l++)
519  volume[l].set(i, k, bin_low);
520 // 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 );
521 
522  if (do_quality)
523  {
524  dato_corrotto(i, k)=ANAP_OK;
525  elev_fin(i, k) = el_inf;
526  }
527  if (el_inf > first_level_static(i, k)) grid_stats.incr_elev(i, k);//incremento la statistica cambio elevazione
528  flag_anap = false;
529  }
530  }/* test_anap */
531  //----------------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
532  else if (bin_low < fondo_scala)
533  {
534  for(unsigned l=0; l<el_up; l++)
535  {
536  if (volume[l].beam_size > k && volume[el_up].beam_size > k)
537  volume[l].set(i, k, bin_high);
538  else if (volume[l].beam_size > k )
539  volume[l].set(i, k, fondo_scala);
540 
541  }
542 // 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);
543  //----------------controlli su bin_high nel caso in cui bin_low sia un no data per assegnare matrice anap (dato_corrotto(i, k))
544  if (do_quality)
545  {
546  if (bin_high<fondo_scala) dato_corrotto(i, k)=ANAP_NODAT;/*manca dato sotto e sopra*/
547  bool test_an1;
548  if (cont_anap< THR_CONT_ANAP )
549  test_an1=(bin_high>=fondo_scala); //modificato per contemplare > o >=
550  else
551  test_an1=(bin_high>fondo_scala);
552 
553  if (test_an1) dato_corrotto(i, k)=ANAP_NOCONTROL;
554  if (bin_high==fondo_scala) dato_corrotto(i, k)=ANAP_OK;/*non piove (oppure sono sopra livello preci...)*/
555  }
556 
557  if (do_beamblocking)
558  beam_blocking(i, k)=0;
559  }
560 
561  //----------------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)
562  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*/
563 
564  {
565  unsigned count =0;
566  for (unsigned ii=0; ii<7; ii++){
567  int iaz = (i + ii - 3 + volume.beam_count) % volume.beam_count;
568  if( SD[el_inf].get(iaz,k) < conf_texture_threshold && SD[el_inf].get(iaz,k) > 0.01) count++;
569  }
570  if ( !(SD[el_inf].get(i,k) < conf_texture_threshold && SD[el_inf].get(i,k) >0.01 && count >=5 ))
571  bin_low = fondo_scala;
572 
573  for(unsigned l=0; l<=el_inf; l++)//riempio con i valori di el_inf tutte le elevazioni sotto (ricostruisco il volume)
574  {
575  if (volume[l].beam_size > k)
576  volume[l].set(i, k, bin_low);
577  }
578 // 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);
579 
580  if (do_quality)
581  {
582  dato_corrotto(i, k)=ANAP_OK; // dubbio
583  elev_fin(i, k) = el_inf;
584  }
585 
586  if (el_inf > first_level_static(i, k)) grid_stats.incr_elev(i, k);
587  }
588  /*-----------------------------------------------------------fine di tutti gli if-----------*/
589  //-----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.
590 
591  if (do_quality)
592  quota(i, k)=(unsigned short)PolarScanBase::sample_height(
593  elev_fin.elevation_at_elev_preci(i, k),
594  (k + 0.5) * volume[0].cell_size);
595  }
596  } // end for over beam_count
597 }
598 
599 // Explicit template instantiation for T=double
600 template class Anaprop<double>;
601 
602 }
603 }
String functions.
Definition: cart.cpp:4