Elaboradar 0.1
Caricamento in corso...
Ricerca in corso...
Nessun risultato
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
8using namespace std;
9
10namespace radarelab {
11namespace algo {
12
13Livmin::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
26InstantaneousVPR::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/*
32idx funzione calcolo VPR istantaneo
33calcola il VPR istantaneo secondo il metodo di Germann e Joss (2003)
34Per il calcolo si considerano i punti con Z>THR_VPR, qualità>QMIN_VPR, BeamBlocking<20% e clutter free all'interno del volume scelto.
35Il 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.
36La 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)
37Il 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
43COSTANTI
44THR_VPR:soglia valore Z per calcolo vpr
45QMIN_VPR: qualità minima necessaria per entrare nel calcolo vpr
46MIN_AREA: area minima a un livello per avere un punto del vpr
47IAZ_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
48TCK_VPR: spessore dei livelli
49NMAXLAYER: n max. livelli
50VEXTMIN_VPR: estensione verticale minima del profilo
51THR_TDEVR: soglia di massima dev. standard di R per considerare buono il profilo
52CT_MIN: frazione di volume minimo precipitante per considerare buono il profilo
53
54VARIABILI
55int flag_vpr[][][]: flag che indica l'usabilità della cella polare in termini di posizione del bin, beam blocking e clutter
56long int *cv,*ct: volume del settore libero,volume precipitante
57float vpr1[]: vpr istantaneo
58long int vert_ext,vol_rain: estensione verticale profilo, volume pioggia del singolo bin
59long int area_vpr[NMAXLAYER]; area totale usata per calcolo vpr
60
61*/
62void 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
247namespace {
248float 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
261VPR 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}
Homogeneous volume with a common beam count for all PolarScans.
Definition volume.h:431
float comp_levels(float v0, float v1, float nodata, float peso)
combina livelli
Namespace per volume dati.
Definition elev_fin.h:12
String functions.
Definition cart.cpp:4