Elaboradar  0.1
 Tutto Classi Namespace File Funzioni Variabili Tipi enumerati (enum) Gruppi
steiner.cpp
1 #include "steiner.h"
2 #include "radarelab/par_class.h"
3 #include "radarelab/algo/dbz.h"
4 
5 //#ifdef __cplusplus
6 //extern "C" {
7 //#endif
8 //#include <func_Z_R.h>
9 //#ifdef __cplusplus
10 //}
11 //#endif
12 
13 // valore mancante
14 #define MISSING 0
15 
16 // fattore conversione gradi-radianti
17 #define DTOR M_PI/180.
18 
19 //Definizioni geometriche
20 static const double AMPLITUDE = 0.9; /* esternalizzo?*/ // ampiezza fascio radar
21 
22 using namespace std;
23 
24 namespace radarelab {
25 namespace algo {
26 
27 namespace steiner {
28 
29 void Point::add_sample(double sample)
30 {
31  if (sample <= MINVAL_DB) return;
32  Z_bckgr += DBZ::BYTEtoZ(DBZ::DBtoBYTE(sample));
33  bckgr += sample;
34  ++npoints;
35 }
36 
37 void Point::finalize()
38 {
39  if (npoints > 0){
40  Z_bckgr = Z_bckgr / npoints;
41  //bckgr[i]=bckgr[i]/npoints; //no
42  if (Z_bckgr > 0) bckgr = 10 * (log10(Z_bckgr));
43  }
44  //il valore del raggio convettivo varia a seconda del background, da 1 a 5 km
45  if (bckgr < 25.)
46  convective_radius = 1.;
47  else if (bckgr >= 25. && bckgr < 30.)
48  convective_radius = 2.;
49  else if (bckgr >= 30. && bckgr < 35.)
50  convective_radius = 3.;
51  else if (bckgr >= 35. && bckgr < 40.)
52  convective_radius = 4.;
53  else if (bckgr > 40.)
54  convective_radius = 5.;
55 }
56 
57 }
58 
59 CalcoloSteiner::CalcoloSteiner(
60  const Volume<double>& volume,
61  const volume::ElevFin<double>& elev_fin,
62  unsigned max_bin)
63  // TODO: evaluate cell_size scan by scan?
64  : volume(volume), elev_fin(elev_fin), max_bin(max_bin), size_cell(volume.scan(0).cell_size),
65  conv_STEINER(Matrix2D<unsigned char>::Constant(volume.beam_count, max_bin, MISSING))
66 {
67  using namespace steiner;
68 
69  logging_category = log4c_category_get("radar.vpr");
70 
71  // Traccio una lista dei punti che hanno valore non nullo e sotto base
72  // bright band (lista_bckg) contenente iaz e irange e conto i punti
73  for (unsigned i=0; i < volume.beam_count; ++i)
74  for (unsigned j=0; j < max_bin; ++j) // propongo max_bin visto che risoluzione è la stessa
75  //if ( volume.scan(0)[i][j] > 1 && (float)(quota[i][j])/1000. < hbbb ) //verifico che il dato usato per la ZLR cioè la Z al lowest level sia > soglia e la sua quota sia sotto bright band o sopra bright band
76  if (j < volume.scan(0).beam_size && volume.scan(0).get(i, j) > MINVAL_DB)
77  lista_bckg.push_back(Point(i, j)); // IAZIMUT, IRANGE
78 }
79 
80 void CalcoloSteiner::calcolo_background() // sui punti precipitanti calcolo bckgr . nb LA CLASSIFICAZIONE DI STEINER NON HA BISOGNO DI RICAMPIONAMENTO CILINDRICO PERCIÒ uso direttamente la matrice polare
81  // definisco una lista di punti utili per l'analisi, cioè con valore non nullo e sotto la bright band o sopra (NB. per l'analisi considero i punti usati per la ZLR che si suppongono non affetti da clutter e beam blocking<50%. questo ha tutta una serie di implicazioni..... tra cui la proiezione, il fatto che nel confronto entrano dei 'buchi'...cioè il confronto è fatto su una matrice pseudo-orizzontale bucata e quote variabili tuttavia, considerato che i gradienti verticali in caso convettivo e fuori dalla bright band non sono altissimi spero funzioni ( andro' a verificare che l'ipotesi sia soddisfatta)
82  // per trovare il background uso uno pseudo quadrato anzichè cerchio, mantenendo come semi-lato il raggio di steiner (11KM); devo perciò definire una semi-ampiezza delle finestre in azimut e range corrispondente al raggio di 11 km
83  // quella in azimut dipende dal range
84 
85 {
86  using namespace steiner;
87 
88  if (lista_bckg.size() < 2)
89  return;
90 
91  // Per il calcolo della finestra range su cui calcolare il background
92  // divido il raggio di Steiner (11km) per la dimensione della cella.
93  // Definisco ampiezza semi-finestra range corrispondente al raggio di
94  // steiner (11km), in unità matrice polare.
95  unsigned delta_nr = round(STEINER_RADIUS * 1000. / size_cell);
96  LOG_DEBUG("delta_n range per analisi Steiner = %u --- dimensione lista_bckg %d", delta_nr,lista_bckg.size());
97 
98  for (vector<Point>::iterator i = lista_bckg.begin(); i != lista_bckg.end(); ++i)
99  {
100  // estremi della finestra in range
101  int kmin = i->range - delta_nr;
102  unsigned kmax = min(i->range + delta_nr, max_bin);
103 
104  if (kmin>0)
105  {
106  //definisco ampiezza semi finestra nazimut corrispondente al raggio di steiner (11km) (11/distanzacentrocella)(ampiezzaangoloscansione)
107  unsigned delta_naz=ceil(STEINER_RADIUS/((i->range * size_cell/1000. + size_cell/2000.)/(AMPLITUDE*DTOR)));
108  if (delta_naz > volume.beam_count / 2)
109  delta_naz = volume.beam_count / 2;
110 
111  int jmin = i->azimut - delta_naz;
112  int jmax = i->azimut + delta_naz;
113 
114  for (int j = jmin; j < jmax; ++j)
115  for (unsigned k = kmin; k < kmax; ++k)
116  i->add_sample(elev_fin.db_at_elev_preci((j + volume.beam_count) % volume.beam_count, k));
117  } else {
118  // FIXME: questo fa mezzo scan tra 0 e kmax, e mezzo scan tra 0 e
119  // -kmin. Sempre gli stessi mezzi scan a prescindere dalla
120  // posizione di i. Ha senso?
121  for (unsigned j=0 ; j<volume.beam_count/2 ; j++)
122  for (unsigned k=0 ; k<kmax ; k++)
123  i->add_sample(elev_fin.db_at_elev_preci(j, k));
124 
125  for (unsigned j= volume.beam_count/2 ; j<volume.beam_count ; j++)
126  for (int k=0 ; k<-kmin ; k++)
127  i->add_sample(elev_fin.db_at_elev_preci(j, k));
128  }
129  i->finalize();
130  }
131 }
132 
133 void CalcoloSteiner::ingrasso_nuclei(float cr,int ja,int kr)
134 {
135  int dr=(int)(cr*1000./size_cell);//definisco ampiezza semi-finestra range corrispondente al raggio di steiner (11km), unità matrice polare
136 
137  int kmin=kr-dr;
138  unsigned kmax=kr+dr;
139 
140  int daz=ceil(cr/((kr*size_cell/1000.+size_cell/2000.)/(AMPLITUDE*DTOR)));
141  int jmin=ja-daz;
142  unsigned jmax=ja+daz;
143 
144  LOG_DEBUG("dr cr kmin kmax %d %f %d %d %d %d", dr,cr, kmin,kmax,jmin,jmax);
145 
146  if (kmin>0) {
147  if (kmax > max_bin) kmax = max_bin;
148 
149  if (jmin<0) {
150  jmin=volume.beam_count-jmin%volume.beam_count;
151  for (unsigned j=jmin; j< volume.beam_count ; j++) {
152  for (unsigned k=kmin ; k<kmax ; k++) {
153  conv_STEINER(j, k)=CONV_VAL;
154  }
155  }
156  LOG_DEBUG("jmin %d", jmin);
157  jmin=0;
158 
159  }
160 
161  if (jmax>=volume.beam_count) {
162  jmax=jmax%volume.beam_count;
163  for (unsigned j=0; j<jmax ; j++) {
164  for (unsigned k=kmin; k<kmax ; k++) {
165  conv_STEINER(j, k)=CONV_VAL;
166  }
167  }
168  LOG_DEBUG("jmax %d", jmax);
169  jmax=volume.beam_count;
170  }
171  for (unsigned j=jmin; j<jmax ; j++) {
172  for (unsigned k=kmin; k<kmax ; k++) {
173  conv_STEINER(j%volume.beam_count, k)=CONV_VAL;
174  }
175  }
176  }
177  else
178  {
179  for (unsigned j=0 ; j<volume.beam_count/2 ; j++)
180  for (unsigned k=0 ; k<kmax ; k++){
181  conv_STEINER(j, k)=CONV_VAL;
182  }
183  for (unsigned j= volume.beam_count/2 ; j<volume.beam_count ; j++)
184  for (unsigned k=0 ; k < (unsigned)-kmin ; k++){
185  conv_STEINER(j, k)=CONV_VAL;
186  }
187  }
188  LOG_DEBUG ("Finita ingrasso nuclei");
189  return;
190 }
191 
192 void CalcoloSteiner::classifico_STEINER()
193 {
194  using namespace steiner;
195 
196  for (vector<Point>::const_iterator i = lista_bckg.begin(); i != lista_bckg.end(); ++i)
197  {
198  int j = i->azimut;
199  int k = i->range;
200  if (j < 0 || k < 0) continue;
201 
202  double db = elev_fin.db_at_elev_preci(j, k);
203  // calcolo diff col background
204  float diff_bckgr = db - i->bckgr;
205  // test su differenza con bckground , se soddisfatto e simultaneamente il VIZ non ha dato class convettiva (?)
206  if ((db > 40.) ||
207  (i->bckgr < 0 && diff_bckgr > 10) ||
208  (i->bckgr < 42.43 && i->bckgr > 0 && diff_bckgr > 10. - i->bckgr * i->bckgr / 180.) ||
209  (i->bckgr > 42.43 && diff_bckgr > 0))
210  {
211  // assegno il punto nucleo di Steiner
212  conv_STEINER(j, k) = CONV_VAL;
213 
214  // ingrasso il nucleo
215  float cr = i->convective_radius;
216  LOG_DEBUG(" %f cr", cr);
217  ingrasso_nuclei(cr, j, k);
218  }
219  }
220  return;
221 }
222 
223 }
224 }