Elaboradar  0.1
vpr.cpp
1 #include "vpr.h"
2 #include "radarelab/vpr_par.h"
3 
4 // TODO: should aplitude be computed rather than hardcoded?
5 #define AMPLITUDE 0.9
6 #define DTOR (M_PI/180.)
7 
8 using namespace std;
9 
10 namespace radarelab {
11 namespace algo {
12 
13 Livmin::Livmin(const VPR& vpr)
14 {
15  for (unsigned ilay = 0; ilay < vpr.size(); ++ilay)
16  {
17  if (vpr.val[ilay] <= NODATAVPR) continue;
18  livmin = ilay * TCK_VPR + TCK_VPR / 2;
19  idx = ilay;
20  found = true;
21  break;
22  }
23 }
24 
25 
26 InstantaneousVPR::InstantaneousVPR(const Volume<double>& volume, const Volume<unsigned char>& qual, Volume<unsigned char>& flag_vpr, int az_min, int az_max)
27  : volume(volume), qual(qual), flag_vpr(flag_vpr), az_min(az_min), az_max(az_max), dbz(volume)
28 {
29 }
30 
31 /*
32 idx funzione calcolo VPR istantaneo
33 calcola il VPR istantaneo secondo il metodo di Germann e Joss (2003)
34 Per il calcolo si considerano i punti con Z>THR_VPR, qualità>QMIN_VPR, BeamBlocking<20% e clutter free all'interno del volume scelto.
35 Il punto del vpr corrispondente ad un livello è calcolato come la somma dei volumi di pioggia di ciascuna cella polare presente nell'intervallo di quota, diviso l'area totale. Per ogni livello devo avere un'area minima coperta di pioggia,se non si verifica, ricopio il valore del primo livello che raggiunge questa copertura.
36 La parte del profilo alta, che sta cioè al di sopra dell'ultimo livello calcolato, viene riempite tirando una retta con coeff. amgolre negativo e costante, pari a -0.03 (valore passibile di discussione)
37 Il profilo è poi soggetto a quality check e viene rigettato (return(1)) se:
38 - estensione verticale troppo bassa (devo avere almeno 2 km di profilo da un punto che sta non più alto di 700 m dal suolo)
39 - cv<CT_MIN*ct, cioè frazione di volume precipitante piccola
40 - la deviazione standard del volume di R a 1100 m normalizzata rispetto al volume totale precipitante supera soglia prefissata
41 
42 
43 COSTANTI
44 THR_VPR:soglia valore Z per calcolo vpr
45 QMIN_VPR: qualità minima necessaria per entrare nel calcolo vpr
46 MIN_AREA: area minima a un livello per avere un punto del vpr
47 IAZ_MIN_SPC,IAZ_MAX_SPC,IAZ_MIN_GAT,I_AZ_MAX_GAT,R_MIN_VPR,R_MAX_VPR:limiti dell'area considerata per il calcolo VPR
48 TCK_VPR: spessore dei livelli
49 NMAXLAYER: n max. livelli
50 VEXTMIN_VPR: estensione verticale minima del profilo
51 THR_TDEVR: soglia di massima dev. standard di R per considerare buono il profilo
52 CT_MIN: frazione di volume minimo precipitante per considerare buono il profilo
53 
54 VARIABILI
55 int flag_vpr[][][]: flag che indica l'usabilità della cella polare in termini di posizione del bin, beam blocking e clutter
56 long int *cv,*ct: volume del settore libero,volume precipitante
57 float vpr1[]: vpr istantaneo
58 long int vert_ext,vol_rain: estensione verticale profilo, volume pioggia del singolo bin
59 long int area_vpr[NMAXLAYER]; area totale usata per calcolo vpr
60 
61 */
62 void InstantaneousVPR::compute()
63 {
64  LOG_CATEGORY("radar.vpr");
65  int i,iA,ilay,il,ilast;
66  long int dist,dist_plain,vol_rain;
67  float area,quota_true_st;
68  float grad;
69 
70  //----------------inizializzazioni----------------
71  long int vert_ext=0;
72 
73  /*SECTION A: cv and ct retrieval and calculation of current area- weighted vpr*/
74 
75 
76  //------------riconoscimento sito per definizione limiti azimut---------
77  int iaz_min = az_min;
78  int iaz_max = az_max;
79  LOG_DEBUG(" Iaz_min %d iaz_max %d",iaz_min,iaz_max);
80 
81  for (unsigned l=0; l < volume.size(); ++l)//ciclo elevazioni
82  {
83  const PolarScan<double>& scan = volume[l];
84 
85  for (unsigned k=0; k < scan.beam_size; k++)/*ciclo range*/
86  {
87  //-------------calcolo distanza-----------
88  // TODO: why the casts to int?
89  dist=k*(long int)(scan.cell_size)+(int)(scan.cell_size)/2.;
90 
91  //-----ciclo settore azimut???
92  for (iA=iaz_min; iA<iaz_max; iA++)//ciclo sulle unità di azimut (0,9°) ---------
93  {
94  //----------ottengo il valore dell'indice in unità di azimut (0,9°) ---------
95  i = (iA + volume.beam_count) % volume.beam_count;
96 
97  //--------calcolo elevazione e quota---------
98  quota_true_st = scan.sample_height_real(i, k);
99 
100  //--------trovo ilay---------
101  ilay=floor(quota_true_st/TCK_VPR);// in teoria quota indipendente da azimuth , in realtà no (se c'è vento)
102 
103  if (ilay <0 || ilay >= NMAXLAYER) {
104  //fprintf(log_vpr,"ilay %d errore\n",ilay);
105  break;
106  }
107 
108 
109  vol_rain=0;
110  // dist=(long int)(dist*cos((float)(cum_bac.volume_at_elev_preci(i, k).teta_true)*CONV_RAD));
111 
112  /* //---------calcolo la distanza proiettata sul piano------------- */
113 
114  const double elevaz = scan.elevations_real(i) * M_PI / 180.;
115  dist_plain=(long int)(dist*cos(elevaz));
116  if (dist_plain <RMIN_VPR || dist_plain > RMAX_VPR )
117  flag_vpr.scan(l).set(i, k, 0);
118 // if (iA == iaz_min) LOG_DEBUG(" k %3d dist %6d dist_plain %6d quota_true_st %8.2f ilay %3d elevaz %5.2f %f", k, dist, dist_plain,quota_true_st, ilay,scan.elevation,elevaz);
119 
120  if (qual.scan(l).get(i, k) < QMIN_VPR) flag_vpr.scan(l).set(i, k, 0);
121 
122  //AGGIUNTA PER CLASS
123 #if 0
124  if(cum_bac.do_class){
125  if(conv(i,k)>= CONV_VAL){
126  flag_vpr.scan(l).set(i, k, 0);
127  }
128  }
129 #endif
130 
131  // ------per calcolare l'area del pixel lo considero un rettangolo dim bin x ampiezzamediafascio x flag vpr/1000 per evitare problemi di memoria?
132  area = scan.cell_size * dist_plain * AMPLITUDE * DTOR * flag_vpr.scan(l).get(i, k)/1000.; // divido per mille per evitare nr troppo esagerato
133 
134  // ------incremento il volume totale di area
135 
136  cv = cv + (long int)(area);
137 
138 
139  //---------------------condizione per incrementare VPR contributo: valore sopra 13dbz, qualità sopra 20 flag>0 (no clutter e dentro settore)------------------
140  double sample = scan.get(i, k);
141  if (sample > THR_VPR && flag_vpr.scan(l).get(i, k) > 0 )
142  {
143  //-------incremento il volume di pioggia = pioggia x area
144  vol_rain=(long int)(dbz.DBZ_to_mp_func(sample)*area);//peso ogni cella con la sua area
145 
146  //-------incremento l'area precipitante totale ct,aggiungendo però,cosa che avevo messo male una THR solo per ct, cioè per il peso
147  if (sample > THR_PDF)
148  ct = ct + (long int)(area);
149 
150  //------se l'area in quello strato è già maggiore di 0 allora incremento il volume dello strato altrimenti lo scrivo ex novo. poi vpr1 andrà diviso per l'area
151  if (vpr.area[ilay]> 0) vpr.val[ilay]=vpr.val[ilay]+(float)(vol_rain);
152  else vpr.val[ilay]=(float)(vol_rain);
153 
154  //------incremento l'area dello strato----------
155  vpr.area[ilay]=vpr.area[ilay]+area;
156  }
157  }
158  }
159  }
160 
161  LOG_INFO("calcolati ct e cv ct= %li cv= %li", ct, cv);
162 
163  /*SECTION B: vpr quality checks and re-normalisation of vpr*/
164 
165  //--------------CONTROLLO DI QUALITA' E NORMALIZZAZIONE DEL PROFILO ISTANTANEO CALCOLATO
166  //-------- se il volume supera quello minimo------
167  if (ct > CT_MIN * cv) {
168 
169  ilast=0;
170  vert_ext=0;
171 
172 
173  //----- calcolo 'estensione verticale del profilo , negli strati dove l'area è troppo piccola assegno NODATAVPR, NORMALIZZO il profilo, e se l'estensione è minore di VEXTMIN_VPR esco-
174 
175 
176  for (ilay=0; ilay<NMAXLAYER; ilay++){
177 
178 
179  LOG_INFO(" ilay %d area_vpr= %ld ct= %ld cv= %ld", ilay, vpr.area[ilay], ct, cv);
180 
181  if (vpr.area[ilay]>=MIN_AREA) {
182  vert_ext=vert_ext+TCK_VPR;
183  vpr.val[ilay]=vpr.val[ilay]/(float)(vpr.area[ilay]);
184 
185  }
186  else
187  {
188  vpr.val[ilay]=NODATAVPR;
189 
190 
191  //---- se incontro un punto vuoto oltre 700 m ( o se sono arrivata alla fine) assegno ilast ed esco dal ciclo
192 
193  /* if (ilast > 0 && vert_ext>VEXTMIN_VPR){ */
194 
195  if (ilay*TCK_VPR+TCK_VPR/2 > MINPOINTBASE || ilay== NMAXLAYER -1){ //cambio il criterio per calcolare estensione minima. devo avere almeno 2 km consecutivi a partire da un punto che stia almeno a 700 m (MINPOINTBASE),
196  LOG_INFO("raggiunta cima profilo");
197  ilast=ilay-1;// c'era errore!!!
198 
199  //---------- raggiunta la cima profilo faccio check immediato sull'estensione verticale
200  if (vert_ext<VEXTMIN_VPR ){
201  LOG_INFO("estensione profilo verticale troppo bassa");
202  ct = 0;
203  ilast=0;
204  for (il=0; il<ilast; il++) vpr.val[il]=NODATAVPR;
205  return;
206  }
207 
208  break; // esco dal ciclo..modifica
209  }
210  }
211  }
212  }// fine se volumeprecipitante sufficiente
213 
214  // ---------se il volume non supera quello minimo esco---------
215  else {
216  LOG_INFO("volume precipitante troppo piccolo");
217  ct = 0;
218  ilast=0;
219  for (il=0; il<NMAXLAYER; il++) vpr.val[il]=NODATAVPR; //--devo riassegnare o mi rimane 'sporco' forse si potrebbe usare una ver diversa
220  return;
221  }
222 
223 
224  //------calcolo il gradiente del profilo come media del gradiente negli ultimi 3 strati per assegnare la parte 'alta' (novità)
225 
226  grad=((vpr.val[ilast]-vpr.val[ilast-1]) + (vpr.val[ilast-1]-vpr.val[ilast-2])/(2.)+ (vpr.val[ilast-2]-vpr.val[ilast-3])/(3.) ) /3.;
227  if (grad > 0.002)
228  grad=-0.03 ;
229  LOG_INFO(" %f", grad);
230 
231  //--riempio la parte alta del profilo decrementando di grad il profilo in ogni strato fino a raggiunere 0, SI PUÒ TOGLIERE E METTERE NODATA
232 
233  for (ilay=ilast+1; ilay<NMAXLAYER; ilay++) {
234  if (vpr.val[ilay-1] + grad > 0.002)
235  vpr.val[ilay]= vpr.val[ilay-1]+grad;
236  else
237  vpr.val[ilay]=0;
238  }
239 
240  // HO CAMBIATO DA GRADIENTE FISSO PARI A (V. VECCHIO) A GRADIENTE RICAVATO DAL PROFILO PER LA PARTE ALTA
241  //---HO TOLTO TUTTA LA PARTE CHE FA IL CHECK SULLA STDEV A 1100 M SI PUO' RIVEDERE SE METTERLA SE SERVE.
242 
243 
244  success = true;
245 }
246 
247 namespace {
248 float comp_levels(float v0, float v1, float nodata, float peso)
249 {
250  float result;
251  /* if ((v0<nodata+1)&&(v1<nodata+1)) result=nodata; */
252  /* if (v0<nodata+1) result=v1; */
253  /* if (v1<nodata+1) result=v0; */
254  if ((v0>nodata) && (v1>nodata) ) result=((1.-peso)*v0+peso*v1); /* in questa configurazione il vpr è di altezza costante nel tempo ma un po' 'sconnesso' in alto*/
255 
256  else result=nodata;
257  return(result);
258 }
259 }
260 
261 VPR combine_profiles(const VPR& _vpr0, const VPR& _vpr1, long int cv, long int ct)
262 {
263  VPR vpr;
264  VPR vpr0 = _vpr0;
265  VPR vpr1 = _vpr1;
266  vpr.area = vpr1.area;
267  /*----calcolo il peso c0 per la combinazione dei profili*/
268  long int c0 = 2 * cv;
269  int n=0,diff=0;
270  //----------------se l'istantaneo c'è o ho trovato un file con cui combinare
271  //-----se ho i due profili riempio parte bassa con differenza media allineandoli e combino poi
272  // calcolo la diff media
273  for (unsigned ilay=0; ilay<NMAXLAYER; ilay++){
274  if ( vpr0.val[ilay]> NODATAVPR && vpr1.val[ilay]>NODATAVPR ){
275  diff=diff + vpr0.val[ilay]-vpr1.val[ilay];
276  n=n+1;
277  }
278  }
279  if (n>0){
280  //------------- trovo livello minimo -------
281  Livmin livmin1(vpr0);
282  Livmin livmin2(vpr1);
283  diff=diff/n;
284  for (unsigned ilay=0; ilay<max(livmin1.livmin/TCK_VPR,livmin2.livmin/TCK_VPR); ilay++){
285  if (vpr0.val[ilay]<= NODATAVPR && vpr1.val[ilay] > NODATAVPR)
286  vpr0.val[ilay]=vpr1.val[ilay]-diff;
287  vpr.area[ilay]=vpr1.area[ilay];
288  if (vpr1.val[ilay]<= NODATAVPR && vpr0.val[ilay] > NODATAVPR)
289  vpr1.val[ilay]=vpr0.val[ilay]+diff;
290  vpr.area[ilay]=vpr0.area[ilay];
291 
292  }
293  }
294  // peso vpr corrente per combinazione
295  float alfat = (float)ct / (c0 + ct);
296  for (unsigned ilay=0; ilay<NMAXLAYER; ilay++){
297  if (vpr0.val[ilay] > NODATAVPR && vpr1.val[ilay] > NODATAVPR)
298  vpr.val[ilay]=comp_levels(vpr0.val[ilay],vpr1.val[ilay],NODATAVPR,alfat);// combino livelli
299  vpr.area[ilay]=vpr1.area[ilay];
300  }
301 
302  return vpr;
303 }
304 
305 }
306 }
float comp_levels(float v0, float v1, float nodata, float peso)
combina livelli
String functions.
Definition: cart.cpp:4