Elaboradar 0.1
Caricamento in corso...
Ricerca in corso...
Nessun risultato
cum_bac.cpp
1#include "cum_bac.h"
2#include <radarelab/logging.h>
3#include <radarelab/utils.h>
4#include <radarelab/image.h>
5#include <radarelab/sp20.h>
6#include <radarelab/odim.h>
10#include <radarelab/algo/dbz.h>
11#include <radarelab/algo/vpr.h>
12#include "site.h"
14#include <radarelab/algo/viz.h>
16#include "cartproducts.h"
17#include <radarelab/algo/top.h>
20#include <radarelab/cart.h>
21#include <radarlib/radar.hpp>
22#include <cstring>
23#include <cstdlib>
24#include <stdexcept>
25#include <cmath>
26#include <iostream>
27#include <unistd.h>
28#include "setwork.h"
29#include <sstream>
30
31//#ifdef __cplusplus
32//extern "C" {
33//#endif
34//#include <func_Z_R.h>
35//#ifdef __cplusplus
36//}
37//#endif
38
39#include <func_Q3d.h>
40
41#include <qual_par.h>
42#include <radarelab/par_class.h>
43
44#ifdef NEL
45#undef NEL
46#endif
47
48// Soglie algoritmi
49#define OVERBLOCKING 51 /* minimo BB non accettato*/
50#define SOGLIA_TOP 45.0 // soglia per trovare top
51#define MISSING 0 /*valore mancante*/
52
53//Definizioni geometriche
54#define AMPLITUDE 0.9 /* esternalizzo?*/ // ampiezza fascio radar
55
56// anaprop
57#define LIMITE_ANAP 240/* LIMITE in numero bins per cambiare controllo anaprop*/
58
59#define DTOR M_PI/180. /* esternalizzo?*/ //fattore conversione gradi-radianti
60#define CONV_RAD 360./4096.*DTOR // fattore conversione unità angolare radar-radianti
61
62using namespace std;
63using namespace radarelab;
64using namespace radarelab::algo;
65
66namespace elaboradar {
67
68namespace {
74struct HRay : public Matrix2D<double>
75{
76 static const int NSCAN = 6;
77
78 // distanza temporale radiosondaggio
79 double dtrs;
80
81 template<typename T>
82 HRay(const Volume<T>& vol) : Matrix2D<double>(vol.size(), vol.max_beam_size())
83 {
84 const double radius = 6378137.0;
85 const double kea = 4. / 3. * radius;
86
87 for (unsigned iel = 0; iel < vol.size(); ++iel)
88 {
89 const double elev = vol.scan(iel).elevation;
90 const double cell_size = vol.scan(iel).cell_size;
91
92 for (unsigned ibin = 0; ibin < cols(); ++ibin)
93 {
94 double range = (ibin + 0.5) * cell_size;
95 (*this)(iel, ibin) = ::sqrt(::pow(range, 2.) + ::pow(kea, 2.) + 2 * kea * range * ::sin(DTOR * elev)) - kea;
96 }
97 }
98 }
99
100 void load_hray(Assets& assets)
101 {
102 // quota centro fascio in funzione della distanza e elevazione
103 dtrs = assets.read_file_hray([this](unsigned el, unsigned bin, double value) {
104 if (el >= rows() || bin >= cols()) return;
105 (*this)(el, bin) = value;
106 });
107 }
108 void load_hray_inf(Assets& assets)
109 {
110 // quota limite inferiore fascio in funzione della distanza e elevazione
111 dtrs = assets.read_file_hray_inf([this](unsigned el, unsigned bin, double value) {
112 if (el >= rows() || bin >= cols()) return;
113 (*this)(el, bin) = value;
114 });
115 }
116};
117}
118
119CUM_BAC::CUM_BAC(Volume<double>& volume, const Config& cfg, const Site& site, bool medium, unsigned max_bin)
120 : MyMAX_BIN(max_bin), cfg(cfg), site(site), assets(cfg),
121 do_medium(medium), volume(volume), SD_Z6(volume.beam_count), cil(volume, 0, RES_HOR_CIL, RES_VERT_CIL),
122 dbz(volume), flag_vpr(volume, 0),
123 first_level(NUM_AZ_X_PPI, MyMAX_BIN), first_level_static(NUM_AZ_X_PPI, MyMAX_BIN),
124 bb_first_level(NUM_AZ_X_PPI, 1024), beam_blocking(NUM_AZ_X_PPI, 1024),
125 anaprop(volume), dem(NUM_AZ_X_PPI, 1024),
126 qual(volume, 0),
127 top(volume.beam_count, volume.max_beam_size())
128{
129 logging_category = log4c_category_get("radar.cum_bac");
131
132 //definisco stringa data in modo predefinito
133 time_t Time = volume.load_info->acq_date;
134 struct tm* tempo = gmtime(&Time);
135 // scrivo la variabile char date con la data in formato aaaammgghhmm
136 strftime(date, 20, "%Y%m%d%H%M", tempo);
137
138 // ------definisco i coeff MP in base alla stagione( mese) che servono per calcolo VPR e attenuazione--------------
139 algo::compute_top(volume, SOGLIA_TOP, top);
140}
141
142CUM_BAC::~CUM_BAC()
143{
144 if (calcolo_vpr) delete calcolo_vpr;
145}
146
148{
149 calcolo_vpr = new CalcoloVPR(*this);
150}
151
152void CUM_BAC::read_sp20_volume(Volume<double>& volume, const Site& site, const char* nome_file, bool do_clean, bool do_medium)
153{
154 using namespace radarelab::volume;
155 LOG_CATEGORY("radar.io");
156 LOG_INFO("Reading %s for site %s", nome_file, site.name.c_str());
157
158 SP20Loader loader;
159
160 Scans<double> z_volume;
161 Scans<double> w_volume;
162 Scans<double> v_volume;
163 loader.vol_z = &z_volume;
164 loader.vol_w = &w_volume;
165 loader.vol_v = &v_volume;
166 loader.load(nome_file);
167
168 // Normalise the scan elevations to match the elevations requested in Site
169 auto elev_array = site.get_elev_array(do_medium);
170 z_volume.normalize_elevations(elev_array);
171 w_volume.normalize_elevations(elev_array);
172 v_volume.normalize_elevations(elev_array);
173
174 if (do_clean)
175 {
176 for (unsigned i = 0; i < z_volume.size(); ++i) {
177 double bin_wind_magic_number = site.get_bin_wind_magic_number(v_volume.load_info->acq_date)
178 * v_volume.at(i).gain + v_volume.at(i).offset;
179 algo::Cleaner::clean(z_volume.at(i), w_volume.at(i), v_volume.at(i), bin_wind_magic_number);
180 }
181 }
182
183 algo::azimuthresample::MaxOfClosest<double> resampler;
184 resampler.resample_volume(z_volume, volume, 1.0);
185// Copy all radar site information to volume data
187
188}
189
190 void CUM_BAC::read_odim_volume(Volume<double>& volume, const Site& site, const char* nome_file, char* fuzzypath, bool do_clean, bool do_medium, bool set_undetect)
191{
192 using namespace radarelab::volume;
193 LOG_CATEGORY("radar.io");
194 namespace odim = OdimH5v21;
195 LOG_INFO("Reading %s for site %s", nome_file, site.name.c_str());
196
197 volume::ODIMLoader loader;
198
199 Scans<double> dbzh_volume;
200 Scans<double> th_volume;
201 Scans<double> v_volume;
202 Scans<double> w_volume;
203 Scans<double> zdr_volume;
204 Scans<double> rhohv_volume;
205 Scans<double> sqi_volume;
206 Scans<double> snr_volume;
207 //Scans<unsigned char> full_volume_cleanID;
208 //Scans<double> full_volume_diffprob;
209 //full_volume_diffprob.quantity="Diffprob";
210
211 string radar_name = site.name.c_str(); // da elaboradar/src/site.cpp : 'SPC' o 'GAT'
212 bool init_sqi = false;
213
214 loader.request_quantity(odim::PRODUCT_QUANTITY_DBZH, &dbzh_volume);
215 loader.request_quantity(odim::PRODUCT_QUANTITY_TH, &th_volume);
216
217 if (do_clean)
218 {
219 loader.request_quantity(odim::PRODUCT_QUANTITY_VRAD, &v_volume);
220 loader.request_quantity(odim::PRODUCT_QUANTITY_WRAD, &w_volume);
221 loader.request_quantity(odim::PRODUCT_QUANTITY_ZDR, &zdr_volume);
222 loader.request_quantity(odim::PRODUCT_QUANTITY_RHOHV,&rhohv_volume);
223 loader.request_quantity(odim::PRODUCT_QUANTITY_SQI,&sqi_volume);
224 loader.request_quantity(odim::PRODUCT_QUANTITY_SNR,&snr_volume);
225 }
226 loader.load(nome_file);
227
228 // FIXME: are they really empty? isn't make_scan called on all of them?
229 if (dbzh_volume.empty() && th_volume.empty())
230 {
231 LOG_ERROR("neither DBZH nor TH were found in %s", nome_file);
232 throw runtime_error("neither DBZH nor TH were found");
233 }
234
235 // Normalise the scan elevations to match the elevations requested in Site
236 auto elev_array = site.get_elev_array(do_medium);
237 for (auto i: loader.to_load){
238 if(!i.second->empty() ) i.second->normalize_elevations(elev_array);
239 }
240 Scans<double>* z_volume;
241 if (!dbzh_volume.empty()) {
242 LOG_WARN(" DBZH found");
243 z_volume = &dbzh_volume;
244 }
245 else {
246 LOG_WARN("no DBZH found: using TH");
247 z_volume = &th_volume;
248 }
249
250 //cout<<"z max ="<<z_volume->at(0).maxCoeff()<<endl;
251 if (do_clean && !w_volume.empty() && !v_volume.empty())
252 {
253 if(sqi_volume.empty()) init_sqi = true;
254 //cout<<"init_sqi"<<init_sqi<<endl;
255 if (zdr_volume.empty())
256 {
257 //caso ZDR e grandezze polarimetriche assenti -> lascio versione master al 16/2/2022
258 volume::Scans<unsigned char> full_volume_cleanID;
259 //for (unsigned i = 0; i < 1; ++i){
260 for (unsigned i = 0; i < z_volume->size(); ++i){
261// radarelab::algo::Cleaner::clean(z_volume->at(i), w_volume.at(i), v_volume.at(i),i);
262 full_volume_cleanID.append_scan(z_volume->at(i).beam_count,z_volume->at(i).beam_size,z_volume->at(i).elevation, z_volume->at(i).cell_size, 0);
263 radarelab::algo::Cleaner::evaluateCleanID(z_volume->at(i), w_volume.at(i), v_volume.at(i),full_volume_cleanID.at(i),i);
264 for (unsigned ibeam=0;ibeam<z_volume->at(i).beam_count; ibeam++)
265 for (unsigned j=0; j<z_volume->at(i).beam_size; j++){
266 if (full_volume_cleanID.at(i)(ibeam,j) != 0) z_volume->at(i)(ibeam,j)=z_volume->at(i).undetect;
267 }
268 }
269 } else {
270 cout<<"applico logica fuzzy"<<endl;
271 //caso ZDR e grandezze polarimetriche presenti--> uso fuzzy logic del branch elaboradar_updating al 16/2/2022
272 //volume::Scans<unsigned char> full_volume_cleanID;
273 //volume::Scans<double>full_volume_diffprob;
274 volume::Scans<unsigned char> full_volume_cleanID;
275
276 unsigned last = z_volume->size() -1;
277 for (unsigned i = 0; i < z_volume->size(); ++i){
278 //cout<<"it="<<i<<endl;
279 //volume::Scans<unsigned char> full_volume_cleanID;
280 volume::Scans<double>full_volume_diffprob;
281 full_volume_cleanID.append_scan(z_volume->at(i).beam_count,z_volume->at(i).beam_size,z_volume->at(i).elevation, z_volume->at(i).cell_size);
282 full_volume_diffprob.append_scan(z_volume->at(i).beam_count,z_volume->at(i).beam_size,z_volume->at(i).elevation, z_volume->at(i).cell_size);
283 //full_volume_cleanID.at(0).setZero();
284
285 volume::Scans<double> Texture;
286 //cout<<"init_sqi = "<<init_sqi<<endl;
287 if(init_sqi){
288 sqi_volume.append_scan(z_volume->at(i).beam_count,z_volume->at(i).beam_size,z_volume->at(i).elevation, z_volume->at(i).cell_size);
289 sqi_volume.at(i).setZero();
290 }
291
292 //calcolo texture di dbzh sulla verticale tra prima elevazione e seconda elevazione:
293 if(i< last){
294 volume::Scans<double> Input,Input2;
295 Input.push_back(z_volume->at(i));
296 Input2.push_back(z_volume->at(i+1));
297 radarelab::volume::textureVD(Input, Input2, Texture, true);
298 Texture.at(0).nodata=65535.;
299 Texture.at(0).undetect=0.;
300 //cout<<"it="<<i<<", Texture size = "<<Texture.size()<<" "<<Texture.at(0).size()<<endl;
301 }
302 else{
303 Texture.clear();
304 //cout<<"it="<<i<<", Texture size = "<<Texture.size()<<endl;
305 Texture.append_scan(z_volume->at(i).beam_count,z_volume->at(i).beam_size,z_volume->at(i).elevation, z_volume->at(i).cell_size);
306 Texture.at(0).setZero();
307 Texture.at(0).nodata=65535.;
308 Texture.at(0).undetect=0.;
309 //cout<<"it="<<i<<", Texture size = "<<Texture.size()<<" "<<Texture.at(0).size()<<endl;
310 }
311
312 //algo::Cleaner::evaluateClassID(z_volume->at(i), w_volume.at(i), v_volume.at(i), zdr_volume.at(i), rhohv_volume.at(i), sqi_volume.at(i), snr_volume.at(i), Texture.at(0), full_volume_cleanID.at(i), v_volume.at(i).undetect , radar_name, i);
313 //modifico il force_bool a true (ultimo parametro)
314
315 algo::Cleaner::evaluateClassID(z_volume->at(i), w_volume.at(i), v_volume.at(i), zdr_volume.at(i), rhohv_volume.at(i), sqi_volume.at(i), snr_volume.at(i), Texture.at(0), full_volume_cleanID.at(i), full_volume_diffprob.at(0), v_volume.at(i).undetect , radar_name, fuzzypath, i, true);
316
317 double new_value=z_volume->at(0).nodata;
318 //provo a commentare per test di uso solo nodata dei punti clssificati come non meteo
319 if (set_undetect) new_value=z_volume->at(i).undetect;
320
321 for (unsigned ii = 0; ii < z_volume->at(i).beam_count; ++ii)
322 for (unsigned ib = 0; ib < z_volume->at(i).beam_size; ++ib) {
323
324 if(full_volume_cleanID.at(i)(ii,ib) ) z_volume->at(i)(ii,ib)= new_value;
325
326 }
327 }
328
329 // commento doppia ripulitura tramite clean() della versione master al 16/2/2022
330 //algo::Cleaner::clean(z_volume->at(i), w_volume.at(i), v_volume.at(i),zdr_volume.at(i),i,true);
331 //algo::Cleaner::clean(z_volume->at(i), w_volume.at(i), v_volume.at(i),zdr_volume.at(i),i+100,true);
332 }
333 }
334
335 //cout<<"arrivo al ponghino"<<endl;
336
337 algo::azimuthresample::MaxOfClosest<double> resampler;
338 resampler.resample_volume(*z_volume, volume, 1.0);
339 cout<<"resampler fatto!!!"<<endl;
340
341 /*
342 printf("fbeam ϑ%f α%f", this->volume.scan(0)[0].teta, this->volume.scan(0)[0].alfa);
343 for (unsigned i = 0; i < 20; ++i)
344 printf(" %d", (int)this->volume.scan(0).get_raw(0, i));
345 printf("\n");
346 */
347
348 /*
349 int numRaggi»···»···= scan->getNumRays();
350 NUM_AZ_X_PPI
351
352 NEL
353
354 se due scan per stessa elecvazione, prendo il primo
355
356 guardare se il passo di griglia è 0.9 o dare errore
357 sennò prendere il beam che ha l'angolo piú vicino
358
359 fill_bin in sp20lib
360
361 leggere DBZH o TH (fare poi DBtoBYTE)
362 */
363
364 /*
365 struct VOL_POL volume.scan(NEL)[NUM_AZ_X_PPI];
366 T_MDB_data_header old_data_header;
367
368 //--------lettura volume------
369 int tipo_dati_richiesti = INDEX_Z;
370 int ier = read_dbp_SP20((char*)nome_file,volume.vol_pol,&old_data_header,
371 tipo_dati_richiesti,volume.nbeam_elev);
372
373 if (ier != OK)
374 LOG_ERROR("Reading %s returned error code %d", nome_file, ier);
375
376 // ----- Test sul volume test_volume....... --------
377 if (!test_volume(file_type))
378 {
379 LOG_ERROR("test_volume failed");
380 return false;
381 }
382 */
383
384 // TODO: look for the equivalent of declutter_rsp and check its consistency
385 // like in test_volume
386}
387
388
390{
391 //-------------leggo mappa statica ovvero first_level (funzione leggo_first_level)------------
393
394 //-------------se definita qualita' leggo dem e altezza fascio (mi servono per calcolare qualità)
395 if (do_quality)
396 {
399 }
400
401 //------------se definito DECLUTTER , non rimuovo anap e riscrivo volume polare facedndo declutter solo con mappa statica.... ancora valido?
402
403 if (do_declutter)
404 {
405 for(unsigned i=0; i<NUM_AZ_X_PPI; i++)
406 for(unsigned k=0; k<volume[0].beam_size; k++)
407 {
408 //---assegno el_inf a mappa statica
409 unsigned el_inf = first_level_static(i, k);
410 //---ricopio valori a mappa statica sotto
411 for(unsigned l=0; l<=el_inf; l++)
412 {
413 // Enrico: cerca di non leggere/scrivere fuori dal volume effettivo
414 if (k >= volume[l].beam_size) continue;
415 if (k < volume[el_inf].beam_size)
416 volume[l].set(i, k, volume[el_inf].get(i, k));
417 else
418 volume[l].set(i, k, MISSING_DB);
419 //------------se definito BEAM BLOCKING e non definito BLOCNOCORR (OPZIONE PER non correggere il beam blocking a livello di mappa statica PUR SAPENDO QUANT'È)
421 volume[l].set(i, k, algo::DBZ::beam_blocking_correction(volume[l].get(i, k), beam_blocking(i, k)));
422 }
423 }
424
425 anaprop.init_elev_fin_static(volume, first_level_static);
426 LOG_INFO("declutter_anaprop completed with static declutter");
427 }
428
429 //------------se non definito DECLUTTER inizio rimozione propagazione anomala al livello mappa dinamica e elaborazioni accessorie
430 else if (do_anaprop)
431 {
432 /* 26-5-2004 : se sono alla 1 o successive elevazioni
433 e range > 60 km cambio le soglie, in modo
434 da evitare di riconoscere come anaprop una pioggia shallow
435 Il criterio diventa: - se la differenza tra Z all'elevazione più bassa della
436 corrente e la Z corrente è <10 dbZ allora
437 rendo inefficaci i limiti di riconoscimento anaprop. */
438
439 //--------ciclo sugli azimut e bins per trovare punti con propagazione anomala----------------
440
441 textureSD(volume,SD_Z6,6000., false);
442
443 // test to define the more appropriate value for textture_threshold for rainy and clutter data
444 std::vector <double> above_0 (4,0);
445 std::vector <double> above_15 (4,0);
446 std::vector <double> above_30 (4,0);
447 std::vector <double> above_40 (4,0);
448 for( unsigned el =0; el <4; el++)
449 for (unsigned iaz=0; iaz<volume[el].beam_count; iaz++)
450 for (unsigned k=80; k < volume[el].beam_size; k ++)
451 if (volume[el].get(iaz,k) > 40.){
452 above_0[el]++;
453 above_15[el]++;
454 above_30[el]++;
455 above_40[el]++;
456 } else if (volume[el].get(iaz,k) > 30.){
457 above_0[el]++;
458 above_15[el]++;
459 above_30[el]++;
460 } else if (volume[el].get(iaz,k) > 15.){
461 above_0[el]++;
462 above_15[el]++;
463 } else if (volume[el].get(iaz,k) > 0.){
464 above_0[el]++;
465 }
466
467 anaprop.do_quality = do_quality;
468 anaprop.do_beamblocking = do_beamblocking;
469 anaprop.do_bloccorr = do_bloccorr;
470 if ( above_15[2]/above_15[0] >= 0.025){
471 if (above_0[1]/above_0[0] >= 0.6 && above_30[2]/above_15[2] <0.15 && above_0[1] >=50000){
472 anaprop.conf_texture_threshold = 5.;
473 LOG_WARN("TEXTURE THRESHOLD USED %4.1f -- 0. %6d %6d %6d %6d -- 15. %6d %6d %6d %6d -- 30. %6d %6d %6d %6d -- 40. %6d %6d %6d %6d", anaprop.conf_texture_threshold, (int)above_0[0], (int)above_0[1], (int)above_0[2], (int)above_0[3], (int)above_15[0], (int)above_15[1], (int)above_15[2], (int)above_15[3], (int)above_30[0], (int)above_30[1], (int)above_30[2], (int)above_30[3], (int)above_40[0], (int)above_40[1], (int)above_40[2], (int)above_40[3] );
475 } else {
476 // anaprop.conf_texture_threshold = 5.;
478 LOG_WARN("THUNDERSTORM %4.1f -- 0. %6d %6d %6d %6d -- 15. %6d %6d %6d %6d -- 30. %6d %6d %6d %6d -- 40. %6d %6d %6d %6d", -9.9, (int)above_0[0], (int)above_0[1], (int)above_0[2], (int)above_0[3], (int)above_15[0], (int)above_15[1], (int)above_15[2], (int)above_15[3], (int)above_30[0], (int)above_30[1], (int)above_30[2], (int)above_30[3], (int)above_40[0], (int)above_40[1], (int)above_40[2], (int)above_40[3] );
479 }
480 } else {
482 LOG_WARN("TEXTURE THRESHOLD USED %4.1f -- 0. %6d %6d %6d %6d -- 15. %6d %6d %6d %6d -- 30. %6d %6d %6d %6d -- 40. %6d %6d %6d %6d", anaprop.conf_texture_threshold, (int)above_0[0], (int)above_0[1], (int)above_0[2], (int)above_0[3], (int)above_15[0], (int)above_15[1], (int)above_15[2], (int)above_15[3], (int)above_30[0], (int)above_30[1], (int)above_30[2], (int)above_30[3], (int)above_40[0], (int)above_40[1], (int)above_40[2], (int)above_40[3] );
483
484 }
485 LOG_INFO("declutter_anaprop completed with anaprop");
486 ScrivoStatistica(anaprop.grid_stats);
487 }
488 else
489 {
490 LOG_WARN("declutter_anaprop completed without doing anything");
491 }
492
493 //---------------------------- Code to plot data from polarMatrix
494 /* Image <unsigned char> toBePlotted (volume[0].beam_size, volume[0].beam_count);
495 for(unsigned i=0; i<volume[0].beam_count; i++)
496 for(unsigned k=0 ; k<volume[0].beam_size; k++){
497 toBePlotted(i,k)= DBtoBYTE(volume[0].get(i, k));
498 }
499 radarelab::write_image(toBePlotted, "/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/Polarplot.png", "PNG");*/
500 LOG_INFO("elabora_Dato completata");
501}
502
504{
506 {
507 // Leggo mappa statica
509 // Allargo per coprire la dimensione del volume
512
513 // copio mappa statica su matrice first_level
515 LOG_INFO("Letta mappa statica");
516 }
517
518 if (do_beamblocking)
519 {
520 // Leggo file elevazioni per BB
522 // Allargo per coprire la dimensione del volume
523 if (bb_first_level.cols() < volume.max_beam_size())
525
526 // Leggo file valore di BB
528 // Allargo per coprire la dimensione del volume
529 if (beam_blocking.cols() < volume.max_beam_size())
531
532 /* Se elevazione clutter statico < elevazione BB, prendi elevazione BB,
533 altrimeti prendi elevazione clutter statico e metti a 0 il valore di BB*/
534 for(unsigned i=0; i < first_level.rows(); ++i)
535 for (unsigned j=0; j < first_level.cols(); ++j)
536 {
537 if (do_bloccorr)
538 {
539 if (first_level_static(i, j)<=bb_first_level(i, j))
540 first_level(i, j)=bb_first_level(i, j);
541 else
542 {
543 beam_blocking(i, j)=0;
545 }
546 } else {
547 if (first_level_static(i, j)>bb_first_level(i, j))
548 beam_blocking(i, j)=0;
549 if (first_level_static(i, j)<bb_first_level(i, j))
550 beam_blocking(i, j)=OVERBLOCKING;
551 }
552 }
553 }
554
555 /*-------------------------------
556 patch per espandere il clutter
557 -------------------------------*/
558 if(do_medium){
559 PolarScan<unsigned char> first_level_tmp(first_level);
560 LOG_INFO(" Dentro patch %u ",MyMAX_BIN);
561 for (unsigned i=NUM_AZ_X_PPI; i<800; i++)
562 for (unsigned j=0; j<MyMAX_BIN; j++)
563 for (unsigned k=i-1; k<i+2; k++)
564 if(first_level(i%NUM_AZ_X_PPI, j) < first_level_tmp(k%NUM_AZ_X_PPI, j))
565 first_level(i%NUM_AZ_X_PPI, j)=first_level_tmp(k%NUM_AZ_X_PPI, j);
566 LOG_INFO(" fine patch %u ",MyMAX_BIN);
567 }
568}
569
570void CUM_BAC::ScrivoStatistica(const algo::anaprop::GridStats& grid_stats)
571{
572 //Definizioni per statistica anap
573 static const int DIM1_ST = 16;
574 static const int DIM2_ST = 13;
575 /*--- numero minimo di celle presenti in un
576 settore per la statistica ---*/
577 static const int N_MIN_BIN = 500;
578
579 int az,ran;
580 unsigned char statistica[DIM1_ST][DIM2_ST];
581 unsigned char statistica_bl[DIM1_ST][DIM2_ST];
582 unsigned char statistica_el[DIM1_ST][DIM2_ST];
583
584 memset(statistica,255,DIM1_ST*DIM2_ST);
585 memset(statistica_bl,255,DIM1_ST*DIM2_ST);
586 memset(statistica_el,255,DIM1_ST*DIM2_ST);
587
588 for(az=0; az<DIM1_ST; az++)
589 for(ran=0; ran<DIM2_ST; ran++){
590 if (grid_stats.count(az, ran) >= N_MIN_BIN)
591 {
592 statistica[az][ran] = grid_stats.perc_anap(az, ran);
593 statistica_bl[az][ran] = grid_stats.perc_bloc(az, ran);
594 statistica_el[az][ran] = grid_stats.perc_elev(az, ran);
595 }
596 }
597 File f_stat;
598
599 if (f_stat.open_from_env("ANAP_STAT_FILE", "a"))
600 {
601 fwrite(date,12,1,f_stat);
602 fwrite(statistica,DIM1_ST*DIM2_ST,1,f_stat);
603 }
604
605 if (f_stat.open_from_env("BLOC_STAT_FILE", "a"))
606 {
607 fwrite(date,12,1,f_stat);
608 fwrite(statistica_bl,DIM1_ST*DIM2_ST,1,f_stat);
609 }
610
611 if (f_stat.open_from_env("ELEV_STAT_FILE", "a"))
612 {
613 fwrite(date,12,1,f_stat);
614 fwrite(statistica_el,DIM1_ST*DIM2_ST,1,f_stat);
615 }
616 return ;
617}
618
619/*
620 comstart caratterizzo_volume
621 idx calcola qualita' volume polare
622 calcola qualita' volume polare
623 NB il calcolo è fatto considerando q=0 al di sotto della mappa dinamica.
624 per ora drrs=dist nche nel caso di Gattatico, mentre dtrs è letto da file
625 si puo' scegliere tra qualita' rispetto a Z e rispetto a R, in realtà per ora sono uguali.
626
627 float quo: quota bin
628 float el: elevazione
629 float rst: raggio equivalente in condizioni standard
630 float dZ: correzione vpr
631 float sdevZ: stand. dev. correzione vpr
632
633//--- nb. non ho il valore di bb sotto_bb_first_level
634comend
635*/
637{
638 LOG_DEBUG("start caratterizzo_volume");
639
640 HRay hray_inf(volume); /*quota limite inferiore fascio in funzione della distanza e elevazione*/
641 hray_inf.load_hray_inf(assets);
642
643 cout<<"sono in caratterizzo volume"<<endl;
644
645 // path integrated attenuation
646 double PIA;
647 // dimensione verticale bin calcolata tramite approcio geo-ottico
648 float dh=1.;
649 // distanza radiosondaggio,
650 float dhst=1.;
651 // tempo dal radiosondaggio
652 float drrs=1.;
653 // distanza dal radar
654 float dist=1.;
655 // beam blocking
656 unsigned char bb=0;
657 // indice clutter da anaprop
658 unsigned char cl=0;
659
660 //----------ciclo su NSCAN(=6), cioè sul numero di elevazioni (nominali) per le quali ho calcolato il beam blocking
661 /* a questo punto servono: bb, cl, PIA, dtrs e drrs radiosond, quota, hsup e hinf beam-----------------*/
662
663 for (unsigned l=0; l<volume.size(); l++)/*ciclo elevazioni*/// VERIFICARE CHE VADA TUTTO OK
664 {
665 const auto& scan = volume[l];
666 for (int i=0; i<NUM_AZ_X_PPI; i++)/*ciclo azimuth*/
667 {
668 const double elevaz = scan.elevations_real(i); //--- elev reale in gradi
669
670 //--assegno PIA=0 lungo il raggio NB: il ciclo nn va cambiato in ordine di indici!
671 PIA=0.;
672
673 for (unsigned k=0; k<scan.beam_size; k++)/*ciclo range*/
674 {
675 double sample = scan.get(i, k);
676
677 //---------distanza in m dal radar (250*k+125 x il corto..)
678 dist = k * scan.cell_size + scan.cell_size / 2.;/*distanza radar */
679
680 //-----distanza dal radiosondaggio (per GAT si finge che sia colocato ..), perchè? (verificare che serva )
681 drrs=dist;
682 /* if (!(strcmp(sito,"GAT")) ) { */
683 /* drrs=dist; */
684 /* } */
685 /* if (!(strcmp(sito,"SPC")) ) { */
686 /* drrs=dist; */
687 /* } */
688
689
690 //assegno la PIA (path integrated attenuation) nel punto e POI la incremento (è funzione dell'attenuazione precedente e del valore nel punto)
691 PIA=dbz.attenuation(DBZ::DBtoBYTE(sample),PIA);
692
693 //------calcolo il dhst ciè l'altezza dal bin in condizioni standard utilizzando la funzione quota_f e le elevazioni reali
694 dhst = PolarScanBase::sample_height(elevaz + 0.45, dist)
695 - PolarScanBase::sample_height(elevaz - 0.45, dist);
696
697 //----qui si fa un po' di mischione: finchè ho il dato dal programma di beam blocking uso il dh con propagazione da radiosondaggio, alle elevazioni superiori assegno dh=dhst e calcolo quota come se fosse prop. standard, però uso le elevazioni nominali
698
699 if (l<volume.size() -1 ) {
700 // differenza tra limite sup e inf lobo centrale secondo appoccio geo-ott
701 dh = hray_inf(l + 1, k) - hray_inf(l, k);
702 }
703 else {
704 // non ho le altezze oltre nscan-1 pero' suppongo che a tali elevazioni la prop. si possa considerare standard
705 dh = dhst;
706 }
707
708 if (l < anaprop.elev_fin(i, k)) {
709 cl=algo::ANAP_YES;
710 bb=BBMAX;
711 } else if (l == anaprop.elev_fin(i, k)) {
712 cl=anaprop.dato_corrotto(i, k); /*cl al livello della mappa dinamica*/
713 bb=beam_blocking(i, k); /*bb al livello della mappa dinamica *///sarebbe da ricontrollare perchè con la copia sopra non è più così
714 } else if (l > anaprop.elev_fin(i, k)) {
715 cl=0; /*per come viene scelta la mappa dinamica si suppone che al livello superiore cl=0 e bb=0*/
716 bb=0; // sarebbe if (l-bb_first_level(i, k) >0 bb=0; sopra all'elevazione per cui bb<soglia il bb sia =0 dato che sono contigue o più però condiz. inclusa
717 }
718
719 //------dato che non ho il valore di beam blocking sotto i livelli che ricevo in ingresso ada progrmma beam blocking e
720 //--------dato che sotto elev_fin rimuovo i dati come fosse anaprop ( in realtà c'è da considerare che qui ho pure bb>50%)
721 //--------------assegno qualità zero sotto il livello di elev_fin (si può discutere...), potrei usare first_level_static confrontare e in caso sia sotto porre cl=1
722 if (l < anaprop.elev_fin(i, k)) {
723 qual.scan(l).set(i, k, 0);
724 cl=2;
725 } else {
726 //--bisogna ragionare di nuovo su definizione di qualità con clutter se si copia il dato sopra.--
727
728 //--calcolo la qualità--
729 // FIXME: qui tronca: meglio un round?
730 qual.scan(l).set(i, k, (unsigned char)(func_q_Z(cl,bb,dist,drrs,hray_inf.dtrs,dh,dhst,PIA)*100));
731 }
732
733 if (qual.scan(l).get(i, k) ==0) qual.scan(l).set(i, k, 1);//????a che serve???
734 if (calcolo_vpr)
735 {
736 /* sezione PREPARAZIONE DATI VPR*/
737 if(cl==0 && bb<BBMAX_VPR ) /*pongo le condizioni per individuare l'area visibile per calcolo VPR, riduco il bb ammesso (BBMAX_VPR=20)*/ //riveder.....?????
738 flag_vpr.scan(l).set(i, k, 1);
739 }
740 }
741 }
742 }
743
744 LOG_DEBUG("End caratterizzo_volume");
745 return;
746}
747
749{
750 LOG_CATEGORY("radar.class");
751 int hmax=-9999;
752
753 /* ;---------------------------------- */
754 /* ; FASE 0 : */
755 /* ;---------------------------------- */
756 // DEFINISCO QUOTE DELLA BASE E DEL TOP DELLA BRIGHT BAND USANDO IL DATO quota del picco DEL PRECEDENTE RUN O, SE NON PRESENTE LA QUOTA DELLO ZERO DA MODELLO
757
758 // Lettura quota massimo da VPR calcolo base e top bright band
759 LOG_INFO("data= %s",cum_bac.date);
760 // calcolo il gap
762 //-- se gap < memory leggo hmax da VPR
763 if (gap<=MEMORY){
765 //---suppongo una semiampiezza massima della bright band di 600 m e definisco htopbb e hbasebb come hmassimo +600 m (che da clima ci sta) e hmassimo -600 m
766 }
767
768 if (hmax >= 0)
769 {
770 hbbb=(hmax-600.)/1000.;
771 htbb=(hmax+600.)/1000.;
772 } else {
773 //-- se gap > memory o se non ho trovato il file
774 // Lettura 0 termico da modello, e calcolo base e top bright band
775 LOG_INFO("leggo 0termico per class da file %s",getenv("FILE_ZERO_TERMICO"));
776 // leggo informazioni di temperatura da modello*/
777 float zeroterm;//zerotermico
778 if (cum_bac.assets.read_0term(zeroterm))
779 {
780 //-- considerato che lo shift medio tra il picco e lo zero è tra 200 e 300 m, che il modello può avere un errore, definisco cautelativamente htbb come quota zero + 400 m e hbbb come quota zero -700 m .
781 htbb=zeroterm/1000. + 0.4; // se non ho trovato il vpr allora uso un range più ristretto, potrebbe essere caso convettivo
782 hbbb=zeroterm/1000. - 1.0;
783 } else {
784 LOG_ERROR("non ho trovato il file dello zero termico");
785 LOG_INFO("attenzione, non ho trovat zero termico ne da vpr ne da radiosondaggio");
786 htbb=0.; // discutibile così faccio tutto con VIZ
787 hbbb=0.;
788 }
789 }
790
791 // se hbasebb è <0 metto 0
792 if (hbbb<0.) hbbb=0.;
793
794 LOG_INFO("calcolati livelli sopra e sotto bright band hbbb=%f htbb=%f",hbbb,htbb);
795
796 const CylindricalVolume& cil = cum_bac.cil;
797
798 // ricampionamento del volume in coordinate cilindriche
799 LOG_DEBUG ("Matrice cilindrica Naz %3d Nrange %4d Nheight %4d", cil.slices.size(), cil.x_size, cil.z_size);
800 //-------------------------------------------------------------------------------------------------------------------------
801 // faccio la classificazione col metodo Vertical Integrated Reflectivity
802 algo::CalcoloVIZ viz(cil, htbb, hbbb, t_ground);
803 viz.classifico_VIZ();
804
805 //classificazione con STEINER
806 // if (hmax > 2000.) {// per evitare contaminazioni della bright band, si puo' tunare
807 // if (hbbb > 500.) {// per evitare contaminazioni della bright band, si puo' tunare
808
809 algo::CalcoloSteiner steiner(cum_bac.volume, cum_bac.anaprop.elev_fin, cil.x_size);
810 steiner.calcolo_background();
811 steiner.classifico_STEINER();
812 // }
813 merge_metodi(steiner, viz);
814 return ;
815}
816
817void CalcoloVPR::merge_metodi(const algo::CalcoloSteiner& steiner, const algo::CalcoloVIZ& viz)
818{
819 for (unsigned j = 0; j < NUM_AZ_X_PPI; ++j)
820 for (unsigned k = 0; k < steiner.max_bin; ++k)
821 if (cum_bac.anaprop.quota(j, k) < hbbb*1000. && steiner.conv_STEINER(j, k) == viz.conv_VIZ(j, k) && steiner.conv_STEINER(j, k) > 0)
822 conv(j,k) = steiner.conv_STEINER(j, k);
823 else
824 if (steiner.conv_STEINER(j, k) == viz.conv_VIZ(j, k) && steiner.conv_STEINER(j, k) > 0 && viz.stratiform(j, k) < 1)
825 conv(j,k) = viz.conv_VIZ(j, k);
826}
827
828//----------ALGORITMO
829/* combina il profilo verticale corrente con quello precedente tramite il metodo di Germann (2003)
830 a) calcolo gap tra ultimo profilo e istante corrente
831 b) se gap < MEMORY leggo profilo
832 c) se gap> MEMORY peso 0 il profilo storico e cerco oltre la data (per casi vecchi)
833 d) faccio func_vpr
834 f) cerco il profilo con cui combinare (->proprio, se gap<MEMORY ->dell'altro radar se gap_res<MEMORY e profile_heating_res=WARM)
835 g) Combino livelli con peso sottostante
836 Dati cv e ct, volume totale e volume precipitante il peso del vpr istantaneo è calcolato come segue:
837 c0=2*cv;
838 peso=(float)ct/(c0+ct)
839 long int c0,cv,ct; costanti di combinazione (v. ref.)
840 h) trovo livello minimo, se livello minimo profilo combinato più alto del precedente calcolo la diff media e sommo al vecchio
841 e) ricalcolo livello minimo
842 float vpr0.val[NMAXLAYER],vpr1.val[NMAXLAYER],vpr.val[NMAXLAYER]; profilo precedente, ultimo e combinato
843 float alfat,noval; peso, nodata
844 FILE *file;
845 int mode,ilay; modalità calcolo profilo (0=combinazione, 1=istantaneo),indice di strato
846*/
847int CalcoloVPR::combina_profili(const InstantaneousVPR& inst_vpr)
848{
849 LOG_CATEGORY("radar.vpr");
850
851 LOG_DEBUG (" modalita %d", MOD_VPR);
852 VPR vpr0;
853 bool combinante; // combinante: variabile che contiene presenza vpr alternativo
854 if (MOD_VPR == 0)
855 {
856 /* MOD_VPR=0: VPR combinato */
857 combinante = cum_bac.assets.find_vpr0(cum_bac.dbz, vpr0, gap);
858
859 for (unsigned i=0; i<vpr0.size(); i++)
860 LOG_DEBUG (" Profilo vecchio - livello %2d valore %6.2f",i,vpr0.val[i]);
861 //----a fine calcolo sul sito in esame stampo il valore del gap
862 LOG_INFO("gap %li",gap);
863 } else {
864 /* MOD_VPR=1: VPR istantaneo */
865 combinante = false;
866 }
867
868 if (combinante)
869 {
870 if (inst_vpr.success)
871 {
872 vpr = combine_profiles(vpr0, inst_vpr.vpr, inst_vpr.cv, inst_vpr.ct);
873 } else {
874 // se il calcolo dell'istantaneo non è andato bene , ricopio l'altro vpr e la sua area
875 vpr = vpr0;
876 }
877 } else {
878 if (inst_vpr.success)
879 {
880 // se il calcolo dell'istantaneo è andato bene ricopio il profilo
881 vpr = inst_vpr.vpr;
882 } else {
883 //-----se è andata male la ricerca dell'altro e anche il calcolo dell'istantaneo esco
884 return 1;
885 }
886 }
887
888 //------------- trovo livello minimo -------
889 Livmin livmin(vpr);
890 LOG_INFO(" livmin %i", livmin.livmin);
891
892 if (livmin.idx >= vpr.size() - 1 || !livmin.found)
893 return (1);
894
895 this->livmin = livmin.livmin;
896
897
898 //-----scrivo il profilo e la sua area-----
900 for (unsigned i=0; i<vpr.size(); i++) LOG_DEBUG (" Profilo nuovo - livello %2d valore %6.2f",i,vpr.val[i]);
901
902 return(0);
903}
904
905int CalcoloVPR::profile_heating(bool has_inst_vpr)
906#include <radarelab/vpr_par.h>
907{
908 LOG_CATEGORY("radar.vpr");
909 //---leggo ultimo file contenente riscaldamento , se non esiste impongo heating=0 (verificare comando)
910 int heating = cum_bac.assets.read_vpr_heating();
911
912 //--una volta letto il file, se il calcolo del vpr è andato bene incremento di uno heating sottraendo però la differenza di date (in quarti d'ora)-1 tra gli ultimi due profili
913 //--lo faccio perchè potrei avere heating più alto del dovuto se ho avuto un interruzione del flusso dei dati
914 //-- heating ha un valore massimo pari WARM dopodichè diventa heating = MEMORY e così resta finchè non sono passati MEMORY istanti non aggiornati ( va bene?)
915 //---se il profil non è stato aggiornato invece decremento la variabile riscaldamento di gap con un minimo pari a 0
916 if (!has_inst_vpr){
917 heating=heating-gap; /*se il profilo non è stato aggiornato, ho raffreddamento, in caso arrivi sotto WARM riparto da 0, cioè serve riscaldamento */
918 }
919 else {
920 heating = heating - (gap-1) + 1; /*se il profilo è stato aggiornato, ho riscaldamento , in caso arrivi sopra WARM riparto da MEMORY */
921 if (heating>=WARM) heating=MEMORY; /* se heating raggiunge WARM allora lo pongo uguale a MEMORY */
922 }
923 if (heating<0) heating=0;
924
925 //----stampo heating su file
926 cum_bac.assets.write_vpr_heating(heating);
927
928 //----stampo log vpr
929 LOG_INFO("gap %li heating %i",gap,heating);
930
931 return(heating);
932}
933
934
936{
937 float vpr_dbz;
938
939 File file(logging_category);
940 file.open_from_env("VPR_ARCH", "wt", "ultimo vpr in dBZ per il plot");
941
942 fprintf(file," QUOTA DBZ AREA PRECI(KM^2/1000)\n" );
943 for (unsigned ilay=0; ilay<NMAXLAYER; ilay++){
944 if (vpr.val[ilay]> 0.001 ) {
945 vpr_dbz=cum_bac.dbz.RtoDBZ(vpr.val[ilay]);
946 fprintf(file," %i %10.3f %li\n", ilay*TCK_VPR+TCK_VPR/2, vpr_dbz, vpr.area[ilay]);
947 }
948 else
949 fprintf(file," %i %10.3f %li\n", ilay*TCK_VPR+TCK_VPR/2, NODATAVPR, vpr.area[ilay]);
950 }
951
952 return 0;
953}
954/*=======================================================================================*/
955/*
956 comstart corr_vpr
957 idx corregge i dati tramite profilo verticale
958 corregge i dati tramite profilo verticale a partire da quelli con valore maggiore di THR_CORR(v vpr_par.h): riporto il dato alla quota del livello 'liquido' tramite "traslazione" cioè aggiungo al valore in dBZ la differenza tra il valore del VPR alla quota 'liquida' e quello alla quota della misura, anche in caso di neve,purchè esista il livello liquido nel profilo. In questo caso pero' flaggo il bin tramte la variabile neve[][]. In caso il profilo intero sia di neve allora riporto al valore di Z al suolo (o al livello rappresentativo) perchè non ho un valore di riferimento 'liquido'.
959 La correzione avviene solo se heating>warm
960
961 int ilref : livello del suolo o della quota rappresentativa di esso nel VPR (dove considero buoni i dati a partire dati da REPR_LEV)
962 int ilray : livello a cui ho il dato
963 int ilay2 : livello sopra o sotto ilray a seconda che il fascio stia sopra o sotto il centro di ilray, serve per interpolare il dato vpr su retta ed evitare correzioni a gradino
964 int heating,warm: indicano quanto è caldo il profilo, e la soglia di riscaldamento
965 int snow : indica se il profilo è di neve (snow=1)
966 int neve[NUM_AZ_X_PPI][MAX_BIN]: 1 se c'è neve, identificata se quota dem> hvprmax+300m
967 float corr: correzione in dB
968 float vpr_liq: valore di R nel VPR alla quota 'liquida' ricavato dalla funzione analyse_VPR
969
970 ilref= (dem[i][k]>REPR_LEV)?(floor(dem[i][k]/TCK_VPR)):( floor(REPR_LEV/TCK_VPR)); in pratica livello dem se > REPR_LEV, livello di REPR_LEV altrimenti.
971 ilray=floor((hbin[i][k])/TCK_VPR)
972 corr=RtoDBZ(vpr_liq)-RtoDBZ(vpr_hray)
973 comend
974*/
976 //* ====correzione profilo====================================*/
977
978#include <radarelab/vpr_par.h>
979
980{
981 LOG_CATEGORY("radar.vpr");
982
983 int ilray,ilref,ilay2,ier_ana,snow,strat;
984 float corr,vpr_liq,vpr_hray,hbin,hliq;
985
986 /*inizializzazione variabili */
987 snow=0;
988 //vpr al livello liquido liv liquido e liv max
989 vpr_liq=NODATAVPR;
990 hliq=NODATAVPR;
991 hvprmax=INODATA;
992
993 // analisi vpr
994
995 ier_max=trovo_hvprmax(&hvprmax);
996 ier_ana=analyse_VPR(&vpr_liq,&snow,&hliq);
997 LOG_INFO("ier_analisi %i",ier_ana) ;
998
999 /* se analisi dice che non è il caso di correggere non correggo (NB in questo caso non riempio la matrice di neve)*/
1000 if (ier_ana) return 1;
1001
1002 LOG_INFO("altezza bright band %i",hvprmax);
1003 LOG_INFO("CORREGGO VPR");
1004
1005 //correzione vpr
1006 for (unsigned i=0; i<NUM_AZ_X_PPI; i++){
1007 for (unsigned k=0; k<cum_bac.volume[0].beam_size; k++){
1008 corr=0.;
1009 /* trovo elevazione reale e quota bin*/
1010 //elevaz=(float)(volume_at_elev_preci(i, k).teta_true)*CONV_RAD;
1011 hbin=(float)cum_bac.anaprop.quota(i, k);
1012
1013 /* se dall'analisi risulta che nevica assegno neve ovunque*/
1014 if (snow) neve(i, k)=1;
1015 strat=1;
1016 if (cum_bac.do_class)
1017 {
1018 if (conv(i,k) >= CONV_VAL){
1019 strat=0;
1020 }
1021 }
1022 //--- impongo una soglia per la correzione pari a 0 dBZ
1023 if (cum_bac.volume[0].get(i, k) > THR_CORR && hbin > hliq && strat){
1024
1025 //---trovo lo strato del pixel, se maggiore o uguale a NMAXLAYER lo retrocedo di 2, se minore di livmn lo pongo uguale a livmin
1026 ilray=(hbin>=livmin)?(floor(hbin/TCK_VPR)):(floor(livmin/TCK_VPR));//discutibile :livello del fascio se minore di livmin posto=livmin
1027 if (ilray>= NMAXLAYER ) ilray=NMAXLAYER-2;//livello del fascio se >= NMAXLAYER posto =NMAXLAYER-2
1028
1029 //---trovo ilay2 strato con cui mediare per calcolare il vpr a una quota intermedia tra 2 livelli, se l'altezza del bin è sopra metà strato prendo quello sopra altrimenti quello sotto
1030 if ((int)hbin%TCK_VPR > TCK_VPR/2) ilay2=ilray+1;
1031 else ilay2=ilray-1;
1032 if (ilay2< floor(livmin/TCK_VPR)) ilay2=floor(livmin/TCK_VPR);
1033
1034 //trovo ilref: livello di riferimento per ricostruire il valore vpr al suolo nel caso di neve.
1035 // in caso di profilo di pioggia mi riporto sempre al valore del livello liquido e questo può essere un punto critico.. vedere come modificarlo.
1036
1037 ilref=(cum_bac.dem(i, k)>livmin)?(floor(cum_bac.dem(i, k)/TCK_VPR)):(floor(livmin/TCK_VPR));//livello di riferimento; se livello dem>livmin = livello dem altrimenti livmin
1038
1039
1040 if (vpr.val[ilref] > 0 && vpr.val[ilray] > 0 ){ /*devo avere dati validi nel VPR alle quote considerate!*/
1041 //-- calcolo il valore del profilo alla quota di interesse
1042 vpr_hray=vpr.val[ilray]+((vpr.val[ilray]-vpr.val[ilay2])/(ilray*TCK_VPR-TCK_VPR/2-ilay2*TCK_VPR))*(hbin-ilray*TCK_VPR-TCK_VPR/2); /*per rendere la correzione continua non a gradini */
1043 //--identifico le aree dove nevica stando alla quota teorica dello zero termico
1044
1045 if (cum_bac.dem(i, k)> hvprmax+HALF_BB-TCK_VPR || snow){ /*classifico neve*/
1046 neve(i, k)=1;
1047
1048 }
1049
1050 //--se nevica la correzione consiste solo nel riportare il valore del vpr al suolo: PROPOSTA: qui si potrebbe generare una mappa di intensità di neve ma deve essere rivisto tutto
1051
1052
1053 //if(snow) //A rimosso, faccio una cosa diversa
1054 if(neve(i, k)){
1055
1056 //faccio la regressione lineare dei punti del profilo sopra il punto del dem
1057 //calcolo il valore al livello del dem e lo sostituisco a vpr.val[ilref] nella correzione
1058 // faccio linearizzazione in maniera becera:
1059 //vpr.val[ilref]=(vpr.val[ilref+7]-vpr.val[ilref+2])/(5)*(ilref-(ilref+2))+vpr.val[ilref+2];
1060
1061 //passaggio=BYTEtoR(volume.vol_pol,aMP_SNOW,bMP_SNOW)
1062
1063 //volpol[0][i][k]=RtoBYTE(passaggio)
1064
1065 corr=cum_bac.dbz.RtoDBZ(vpr.val[ilref])-cum_bac.dbz.RtoDBZ(vpr_hray);
1066
1067 cum_bac.volume[0].set(i, k, cum_bac.dbz.DBZ_snow(cum_bac.volume[0].get(i, k)));
1068 }
1069 else{
1070 // -- altrimenti correggo comunque a livello liquido :
1071 corr = cum_bac.dbz.RtoDBZ_class(vpr_liq) - cum_bac.dbz.RtoDBZ_class(vpr_hray);/*riporto comunque al valore liquido anche se sono sopra la bright band*/
1072 }
1073 // -- controllo qualità su valore correzione
1074 if (corr>MAX_CORR) corr=MAX_CORR; /*soglia sulla massima correzione*/
1075 if (hbin<hvprmax && corr>0.) corr=0; /*evito effetti incrementi non giustificati*/
1076
1077 //controllo qualità su valore corretto e correzione
1078 double corrected = cum_bac.volume[0].get(i, k) + corr;
1079 if (corrected > MAXVAL_DB) // se dato corretto va fuori scala assegno valore massimo
1080 cum_bac.volume[0].set(i, k, MAXVAL_DB);
1081 else if ( corrected < MINVAL_DB) // se dato corretto va a fodoscala assegno valore di fondo scala
1082 cum_bac.volume[0].set(i, k, MINVAL_DB);
1083 else
1084 cum_bac.volume[0].set(i, k, corrected); // correggo
1085
1086 corr_polar(i, k)=(unsigned char)(corr)+128;
1087
1088 //inserisco un ponghino per rifare la neve con aMP e bMP modificati // DA SCOMMENTARE SE DECIDO DI FARLO
1089
1090 //if (neve[i][k]) volume.scan(0).get_raw(i, k)=DBtoBYTE(RtoDBZ( BYTE_to_mp_func(volume.scan(0).get_raw(i, k),aMP_SNOW,bMP_SNOW),aMP_class,bMP_class )) ;
1091
1092
1093 }
1094 }
1095 }
1096 }
1097 return(0);
1098}
1099
1101{
1102 int i,imax,istart,foundlivmax;
1103 float h0start,peak,soglia;
1104
1105
1106 if (t_ground != NODATAVPR)
1107 {
1108 LOG_DEBUG("trovo hvprmax a partire da 400 m sotto lo zero dell'adiabatica secca");
1109 h0start=t_ground/9.8*1000 ;
1110 istart=h0start/TCK_VPR -2;
1111 if (istart< livmin/TCK_VPR) istart=livmin/TCK_VPR;
1112 LOG_DEBUG("t_ground h0start istart %f %f %i",t_ground,h0start,istart);
1113 }
1114 else {
1115 LOG_DEBUG("trovo hvprmax a partire da livmin");
1116 istart=livmin/TCK_VPR+1;
1117 }
1118
1119
1120 /* trovo hvprmax e il suo livello a partire dal livello istart */
1121
1122 //--inizializzazione
1123 foundlivmax=0;
1124 peak=NODATAVPR;
1125 *hmax=INODATA;
1126 // Enrico vprmax=NODATAVPR;
1127 imax=INODATA;
1128 soglia = DBZ::DBZtoR(THR_VPR,200,1.6); // CAMBIATO, ERRORE, PRIMA ERA RtoDBZ!!!!VERIFICARE CHE IL NUMERO PARAMETRI FUNZIONE SIA CORRETTO
1129
1130 //--se vpr al livello corrente e 4 layer sopra> soglia, calcolo picco
1131 LOG_DEBUG(" istart %d low %6.2f up %6.2f soglia %6.2f peak %6.2f imax %d", istart, vpr.val[istart] , vpr.val[istart+4], soglia, peak, imax);
1132 if (vpr.val[istart] >soglia && vpr.val[istart+4] > soglia){
1133 peak=10*log10(vpr.val[istart]/vpr.val[istart+4]);//inizializzo il picco
1134 LOG_DEBUG("peak1 = %f",peak);
1135 }
1136 //----se picco > MINIMO il punto è ok
1137 if(peak> MIN_PEAK_VPR){
1138 imax=istart;
1139 // Enrico vprmax=vpr.val[imax];
1140 LOG_DEBUG("il primo punto soddisfa le condizioni di picco");
1141 }
1142 for (i=istart+1;i<NMAXLAYER-4;i++) //la ricerca è un po' diversa dall'originale.. trovo il picco + alto con valore rispetto a 4 sopra > soglia
1143 {
1144 if (vpr.val[i] <soglia || vpr.val[i+4] < soglia) break;
1145 peak=10*log10(vpr.val[i]/vpr.val[i+4]);
1146 if (vpr.val[i]>vpr.val[i-1] && peak> MIN_PEAK_VPR ) // se vpr(i) maggiore del massimo e picco sufficientemente alto
1147 {
1148 imax=i;
1149 // Enrico vprmax=vpr.val[imax];
1150 }
1151 LOG_DEBUG(" low %6.2f up %6.2f soglia %6.2f peak %6.2f imax %d", vpr.val[i] , vpr.val[i+4], soglia, peak, imax);
1152 }
1153
1154 if ( imax > INODATA ){
1155 foundlivmax=1;
1156 peak=10*log10(vpr.val[imax]/vpr.val[imax+4]);
1157 *hmax=imax*TCK_VPR+TCK_VPR/2;
1158 LOG_DEBUG("trovato ilaymax %i %i",*hmax,imax);
1159 LOG_DEBUG(" picco in dbR %f",peak);
1160 }
1161
1162 LOG_DEBUG("exit status trovo_hvprmax %i",foundlivmax);
1163 return (foundlivmax);
1164}
1165
1166/*
1167 1)hvprmax - semiampiezza BB > liv 100 m => la bright band sta sopra il suolo => interpolo il profilo per trovare il livello liquido
1168 se interpolazione fallisce NON CORREGGO (scelta cautelativa, potrei avere caso convettivo
1169 o pochi punti o molto disomogeneo)
1170 2)hvprmax - semiampiezza BB < liv 100 m => la bright contiene o è sotto il livello 100 m oppure ho un massimo 'spurio'
1171 quindi calcolo la Tground, temperatura media nelle stazioni più vicine al radar, se non trovo T torno al caso 1.
1172 2A) Tground>T_MAX_ML ----> prob. caso convettivo o max preci passa sopra radar, interpolo il profilo e calcolo livello liquido
1173 se interpolazione fallisce NON CORREGGO
1174 2B) T_MIN_ML<T0<T_MAX_ML : prob. bright band al suolo, interpolo il profilo per trovare il livello liquido
1175 se interpolazione fallisce NON CORREGGO
1176 2C) Tground<T_MIN_ML ----> prob. profilo neve NON interpolo e non calcolo vpr al livello liquido perchè non ho livello liquido
1177
1178 comend
1179*/
1180int CalcoloVPR::analyse_VPR(float *vpr_liq,int *snow,float *hliq)
1181 /*=======analisi profilo============ */
1182{
1183 int ier=1,ier_ana=0,liv0;
1184 char date[]="000000000000";
1185 struct tm *tempo;
1186 time_t Time;
1187
1188 // ------------inizializzazione delle variabili ----------
1189
1190 //strcpy(date,"000000000000");
1191
1192 int tipo_profilo=-1;
1193 float v600sottobb=NODATAVPR;
1194 float v1000=NODATAVPR;
1195 float v1500=NODATAVPR;
1196 float vliq=NODATAVPR;
1197 float vhliquid=NODATAVPR;
1198 float vprmax=NODATAVPR;
1199 //*togliere gli ultimi tre*/;
1200
1201 //ier_max=trovo_hvprmax(&hvprmax);
1202
1203
1204 if (t_ground == NODATAVPR) //1 se non ho nè T nè il massimo esco altrimenti tipo_profilo=0
1205 {
1206 LOG_WARN("non ho T,...");
1207
1208 if ( ! ier_max ) {
1209 LOG_ERROR(" non ho trovato hvprmax, nè T, esco");
1210 return 1;
1211 }
1212 tipo_profilo=0;
1213 }
1214 else
1215 {
1216
1217 if (t_ground >= T_MAX_ML+0.65*(float)(livmin+TCK_VPR/2)/100.){ //1 se T > T_MAX_ML e non ho il massimo esco
1218 if ( ! ier_max ) {
1219 LOG_ERROR(" temperatura alta e non ho trovato hvprmax, esco");
1220 return 1;
1221 }
1222 tipo_profilo=0;
1223 }
1224
1225
1226 // if (t_ground >= T_MAX_SN+0.65*(float)(livmin+TCK_VPR/2)/100 && t_ground < T_MAX_ML+0.65*(float)(livmin+TCK_VPR/2)/100. )
1227 if (t_ground >= T_MAX_SN && t_ground < T_MAX_ML+0.65*(float)(livmin+TCK_VPR/2)/100. )
1228 {
1229
1230 if ( ier_max ) {
1231 LOG_INFO(" temperatura da scioglimento e massimo in quota");
1232 tipo_profilo=2;
1233 }
1234 else{
1235 LOG_ERROR(" temperatura da scioglimento ma superiore a temperatura max neve e non ho trovato hvprmax, esco");
1236 return 1;
1237 }
1238 // solo una scritta per descrivere cos'è accaduto
1239 liv0=livmin+HALF_BB;
1240 if (hvprmax > liv0) LOG_INFO(" il livello %i è sotto la Bright band, ma T bassa interpolo",livmin);
1241 else LOG_INFO(" il livello %i potrebbe essere dentro la Bright Band, interpolo",livmin);
1242
1243 }
1244
1245 //if (t_ground >= T_MIN_ML && t_ground < T_MAX_SN+0.65*(float)(livmin+TCK_VPR/2)/100.)
1246 if (t_ground < T_MAX_SN)
1247 {
1248 if ( ier_max ){
1249 LOG_INFO(" temperatura da neve o scioglimento e massimo in quota");
1250 tipo_profilo=2;
1251 }
1252 else {
1253 LOG_INFO(" temperatura da neve o scioglimento e massimo non trovato,neve , non interpolo");
1254 tipo_profilo=3;
1255 hvprmax=0;
1256 }
1257 }
1258
1259 }
1260
1261 // InterpolaVPR_NR iv;
1262 InterpolaVPR_GSL iv;
1263
1264 switch
1265 (tipo_profilo)
1266 {
1267 case 0:
1268 case 1:
1269 case 2:
1270 ier=iv.interpola_VPR(vpr.val.data(), hvprmax, livmin);
1271 if (ier){
1272 LOG_INFO(" interpolazione fallita");
1273 switch (tipo_profilo)
1274 {
1275
1276 /*Questo fallisce se la bright band non è al suolo (per evitare correzioni dannose in casi poco omogenei)*/
1277 case 0:
1278 case 1:
1279 ier_ana=1;
1280 *vpr_liq=NODATAVPR;
1281 *hliq=NODATAVPR;
1282 break;
1283 case 2:
1284 *vpr_liq=vpr.val[(hvprmax+1000)/TCK_VPR]*2.15;/*21 aprile 2008*/
1285 *hliq=0;
1286 break;
1287 }
1288 }
1289 else{
1290 LOG_INFO(" interpolazione eseguita con successo");
1291 //
1292 // stampa del profilo interpolato
1293 const char* vpr_arch = getenv("VPR_ARCH");
1294 if (!vpr_arch) throw runtime_error("VPR_ARCH is not defined");
1295 string fname(vpr_arch);
1296 fname += "_int";
1297 File file(logging_category);
1298 file.open(fname, "wt", "vpr interpolato");
1299 for (unsigned i = 0; i < NMAXLAYER; ++i)
1300 fprintf(file," %f \n", cum_bac.dbz.RtoDBZ(iv.vpr_int[i]));
1301
1302 /*calcolo valore di riferimento di vpr_liq per l'acqua liquida nell'ipotesi che a[2]=quota_bright_band e a[2]-1.5*a[3]=quota acqua liquida*/
1303 if (tipo_profilo == 2 ) {
1304 *hliq=(iv.E-2.1*iv.G)*1000.;
1305 //lineargauss(a[2]-2.1*a[3], a, vpr_liq, dyda, ndata);
1306 if (*hliq<0)
1307 *hliq=0; /*con casi di bright band bassa.. cerco di correggere il più possibile*/
1308 *vpr_liq=vpr.val[(hvprmax+1000)/TCK_VPR]*2.15;
1309 }
1310 else {
1311 *hliq=(iv.E-2.1*iv.G)*1000.;
1312 //lineargauss(a[2]-2.1*a[3], a, vpr_liq, dyda, ndata);
1313 if ( *hliq > livmin) {
1314 *vpr_liq=vpr.val[(int)(*hliq/TCK_VPR)]; // ... SE HO IL VALORE VPR USO QUELLO.
1315 }
1316 else // altrimenti tengo il valore vpr neve + 6 dB* e metto tipo_profilo=2
1317 {
1318 if (*hliq<0) *hliq=0;
1319 tipo_profilo=2;
1320 //*vpr_liq=vpr.val[(hvprmax+1000)/TCK_VPR]*2.15;
1321 // iv.C = risultato dell'interpolazione (interpolated vpr)
1322 *vpr_liq=iv.C;
1323 }
1324 }
1325 }
1326 break;
1327 case 3:
1328 *snow=1;
1329 *vpr_liq=NODATAVPR;
1330 *hliq=NODATAVPR;
1331 break;
1332 }
1333 LOG_INFO("TIPO_PROFILO= %i vpr_liq %f hliq %f", tipo_profilo, *vpr_liq,*hliq );
1334
1335
1336 /* parte di stampa test vpr*/
1337
1338 /* nome data */
1339 //definisco stringa data in modo predefinito
1340 Time = cum_bac.volume.load_info->acq_date;
1341 tempo = gmtime(&Time);
1342 sprintf(date,"%04d%02d%02d%02d%02d",tempo->tm_year+1900, tempo->tm_mon+1,
1343 tempo->tm_mday,tempo->tm_hour, tempo->tm_min);
1344 if (! ier ) {
1345 if(*hliq > livmin +200 )
1346 vhliquid=cum_bac.dbz.RtoDBZ(vpr.val[(int)(*hliq)/TCK_VPR]);
1347 vliq=cum_bac.dbz.RtoDBZ(*vpr_liq);
1348 }
1349 if (ier_max) {
1350 if ( hvprmax-600 >= livmin )
1351 v600sottobb=cum_bac.dbz.RtoDBZ(vpr.val[(hvprmax-600)/TCK_VPR]);
1352 if ((hvprmax+1000)/TCK_VPR < NMAXLAYER )
1353 v1000=cum_bac.dbz.RtoDBZ(vpr.val[(hvprmax+1000)/TCK_VPR]);
1354 if ((hvprmax+1500)/TCK_VPR < NMAXLAYER )
1355 v1500=cum_bac.dbz.RtoDBZ(vpr.val[(hvprmax+1500)/TCK_VPR]);
1356 vprmax=cum_bac.dbz.RtoDBZ(vpr.val[(hvprmax/TCK_VPR)]);
1357 }
1358
1359 if (FILE* test_vpr=fopen(getenv("TEST_VPR"),"a+"))
1360 {
1361 fprintf(test_vpr,"%s %i %i -1 %f %f %f %f %f %f %f %f %f %f %f %f %f %f \n",date,hvprmax,tipo_profilo,iv.chisqfin,*hliq,vliq,vhliquid,v600sottobb,v1000+6,v1500+6,vprmax,iv.rmsefin,iv.B,iv.E,iv.G,iv.C,iv.F);
1362 fclose(test_vpr);
1363 }
1364
1365 // fine parte di stampa test vpr
1366
1367 //---SCRIVO ALTEZZA MASSIMO PER CLASSIFICAZIONE AL RUN SUCCESSIVO
1368
1370 LOG_INFO("fatta scrittura hmax vpr = %d",hvprmax);
1371
1372 return (ier_ana);
1373}
1374
1376{
1377 for (unsigned i=0; i<NUM_AZ_X_PPI; i++)
1378 for (unsigned k=0; k<volume[0].beam_size; k++)
1379 if (calcolo_vpr->conv(i,k) > 0)
1380 volume[0].set(i, k, dbz.DBZ_conv(volume[0].get(i, k)));
1381}
1382
1384{
1385 //--------------se definita la qualita procedo con il calcolo qualita e del VPR (perchè prendo solo i punti con qual > soglia?)-----------------//
1386 if (do_quality)
1387 {
1388 //-------------calcolo qualita' e trovo il top
1389
1391 /* //---------trovo il top (a X dbZ) */
1392 /* printf ("trovo top \n") ; */
1393 /* ier=trovo_top(); */
1394
1395
1396 //--------------se definita CLASS procedo con classificazione -----------------//
1397 if (do_class)
1398 {
1400 }
1401
1402 //--------------se definito VPR procedo con calcolo VPR -----------------//
1403
1404 if (calcolo_vpr)
1406 }
1407
1408 if (do_class)
1410}
1411
1412void CUM_BAC::generate_maps(CartProducts& products)
1413{
1414 // Generate products and write them out
1415 LOG_INFO("Scrittura File Precipitazione 1X1\n");
1416 if (do_zlr_media)
1417 {
1418 std::function<unsigned char(const vector<double>&)> convert = [this](const vector<double>& samples) {
1419 // Samples are in contains dB (logaritmic values), so mediating
1420 // them is not the correct operation, and we need to convert
1421 // them to Z (linear) values to average them.
1422 // TODO: there may be more efficient way to mediate logaritmic
1423 // values.
1424 double sum = 0;
1425 for (const auto& s: samples)
1426 sum += DBZ::DBZtoZ(s);
1427 unsigned char res = DBZ::DBtoBYTE(DBZ::ZtoDBZ(sum / samples.size()));
1428 // il max serve perchè il valore di MISSING è 0
1429 if (res == 0) return (unsigned char)1;
1430 return res;
1431 };
1432 products.scaled.to_cart_average(volume[0], convert, products.z_out);
1433 } else {
1434 std::function<unsigned char(unsigned, unsigned)> assign_cart =
1435 [this](unsigned azimuth, unsigned range) {
1436 // il max serve perchè il valore di MISSING è 0
1437 unsigned char sample = DBZ::DBtoBYTE(volume[0].get(azimuth, range));
1438 return max(sample, (unsigned char)1);
1439 };
1440 products.scaled.to_cart(assign_cart, products.z_out);
1441 }
1442
1443 std::function<unsigned char(unsigned, unsigned)> assign_cart =
1444 [this](unsigned azimuth, unsigned range) {
1445 // il max serve perchè il valore di MISSING è 0
1446 unsigned char sample = DBZ::DBtoBYTE(volume[0].get(azimuth, range));
1447 return max(sample, (unsigned char)1);
1448 };
1449 products.fullres.to_cart(assign_cart, products.z_fr);
1450
1451 products.scaled.to_cart(top, products.top_1x1);
1452 if (do_quality)
1453 {
1454 const auto& elev_fin = anaprop.elev_fin;
1455 const auto& quota = anaprop.quota;
1456
1457 std::function<unsigned char(unsigned, unsigned)> assign_qual =
1458 [this, &elev_fin](unsigned azimuth, unsigned range) {
1459 const auto& el = elev_fin(azimuth, range);
1460 if (range >= volume[el].beam_size)
1461 return (unsigned char)0;
1462 return qual.scan(el).get(azimuth, range);
1463 };
1464 products.scaled.to_cart(assign_qual, products.qual_Z_1x1);
1465
1466 std::function<unsigned char(unsigned, unsigned)> assign_quota =
1467 [&quota](unsigned azimuth, unsigned range) {
1468 return 128 + round(quota(azimuth, range) / 100.0);
1469 };
1470 products.scaled.to_cart(assign_quota, products.quota_1x1);
1471 products.scaled.to_cart(anaprop.dato_corrotto, products.dato_corr_1x1);
1472
1473 std::function<unsigned char(unsigned, unsigned)> assign_elev_fin = [&elev_fin](unsigned azimuth, unsigned range) {
1474 return elev_fin(azimuth, range);
1475 };
1476 products.scaled.to_cart(assign_elev_fin, products.elev_fin_1x1);
1477
1478 products.scaled.to_cart(beam_blocking, products.beam_blocking_1x1);
1479 }
1480
1481 if (calcolo_vpr)
1482 {
1483 const auto& neve = calcolo_vpr->neve;
1484 std::function<unsigned char(unsigned, unsigned)> assign = [&neve](unsigned azimuth, unsigned range) {
1485 return neve(azimuth, range) ? 0 : 1;
1486 };
1487 products.scaled.to_cart(assign, products.neve_1x1);
1488
1489 products.scaled.to_cart(calcolo_vpr->corr_polar, products.corr_1x1);
1490
1491 if (do_class)
1492 products.scaled.to_cart(calcolo_vpr->conv, products.conv_1x1);
1493 }
1494
1495 if (do_devel)
1496 {
1497 SD_Z6 *= 10.;
1498
1499 SingleCart SC_SD(SD_Z6.max_beam_size());
1500 for (unsigned int i=0; i<SD_Z6.size(); i++){
1501 SC_SD.creo_cart(SD_Z6, i);
1502 std::ostringstream oss;
1503 oss<<"SD_"<<i;
1504 SC_SD.write_out(assets,oss.str());
1505 }
1506 }
1507
1508// return true;
1509}
1510
1511
1513 : cum_bac(cum_bac),
1514 inst_vpr(cum_bac.volume, cum_bac.qual, cum_bac.flag_vpr, cum_bac.site.vpr_iaz_min, cum_bac.site.vpr_iaz_max),
1515 conv(NUM_AZ_X_PPI, cum_bac.volume.max_beam_size(),0),
1516 corr_polar(NUM_AZ_X_PPI, cum_bac.volume.max_beam_size()),
1517 neve(NUM_AZ_X_PPI, cum_bac.volume.max_beam_size())
1518{
1519 logging_category = log4c_category_get("radar.vpr");
1521 htbb=-9999.; hbbb=-9999.;
1522 t_ground=NODATAVPR;
1523
1524 /*
1525 for (int k=0; k<NUM_AZ_X_PPI*MyMAX_BIN;k++ ){
1526 lista_conv[k][0]=-999;
1527 lista_conv[k][1]=-999;
1528 }
1529 */
1530
1532}
1533
1534CalcoloVPR::~CalcoloVPR()
1535{
1536}
1537
1539{
1540 LOG_INFO("processo file dati: %s", cum_bac.volume.load_info->filename.c_str());
1541 printf ("calcolo VPR \n") ;
1542
1543 //VPR // ------------inizializzo hvprmax ---------------
1544
1545 hvprmax=INODATA;
1546
1547 //VPR // ------------chiamo combina profili con parametri sito, sito alternativo ---------------
1548
1549 inst_vpr.compute(); // ho fatto func_vpr, il profilo istantaneo
1550 LOG_INFO("fatta func vpr %s", inst_vpr.success ? "ok" : "errore");
1551 for (unsigned i=0; i<inst_vpr.vpr.size(); i++) LOG_DEBUG (" Profilo istantaneo - livello %2d valore %6.2f",i,inst_vpr.vpr.val[i]);
1552
1553 int ier_comb;
1554
1555 // ier_comb=combina_profili(sito,argv[4]);
1556 ier_comb=combina_profili(inst_vpr);
1557 LOG_INFO ("exit status calcolo VPR istantaneo: %s", inst_vpr.success ? "ok": "fallito") ; // debug
1558 LOG_INFO("exit status combinaprofili: (1--fallito 0--ok) %i ",ier_comb) ; // debug
1559
1560
1561 //VPR // ------------chiamo profile_heating che calcola riscaldamento profilo ---------------
1562
1563 heating=profile_heating(inst_vpr.success);
1564 printf ("heating %i \n", heating);
1565 LOG_INFO("ier_comb %i", ier_comb);
1566
1567 if (inst_vpr.success)
1569
1570 //VPR // ------------se combina profili ok e profilo caldo correggo --------------
1571
1572 if (!ier_comb && heating >= WARM){
1573
1574 int ier=corr_vpr();
1575 LOG_INFO("exit status correggo vpr: (1--fallito 0--ok) %i",ier) ; // debug
1576
1577
1578 //VPR // ------------se la correzione è andata bene e il profilo è 'fresco' stampo profilo con data-------
1579
1580 if ( ! ier && inst_vpr.success)
1582 }
1583}
1584
1585namespace {
1586struct CartData
1587{
1588 Image<double> azimut;
1589 Image<double> range;
1590
1591 CartData(int max_bin=512)
1592 : azimut(max_bin), range(max_bin)
1593 {
1594 for(int i=0; i<max_bin; i++)
1595 for(int j=0; j<max_bin; j++)
1596 {
1597 range(i, j) = hypot(i+.5,j+.5) ;
1598 azimut(i, j) = 90. - atan((j+.5)/(i+.5)) * M_1_PI*180.;
1599 }
1600 }
1601};
1602}
1603
1604SingleCart::SingleCart(unsigned max_bin)
1605 : max_bin(max_bin),
1606 cart(max_bin*2)
1607{
1608}
1609
1610void SingleCart::creo_cart(const Volume <double>& volume, unsigned el_index)
1611{
1612 LOG_CATEGORY("radar.singlecart");
1613
1614 //matrici per ricampionamento cartesiano
1615 //int x,y,irange,az,iaz,az_min,az_max,cont;
1616 int x,y,iaz,az_min,az_max;
1617 float az;
1618 CartData cd(max_bin);
1619
1620 for(unsigned i=0; i<max_bin *2; i++)
1621 for(unsigned j=0; j<max_bin *2; j++)
1622 cart(i, j) = MISSING;
1623
1624 LOG_INFO("Creo_cart - %u", max_bin);
1625
1626 for(unsigned quad=0; quad<4; quad++)
1627 for(unsigned i=0; i<max_bin; i++)
1628 for(unsigned j=0; j<max_bin; j++)
1629 {
1630 unsigned irange = (unsigned)round(cd.range(i, j));
1631 if (irange >= max_bin)
1632 continue;
1633 switch(quad)
1634 {
1635 case 0:
1636 x = max_bin + i;
1637 y = max_bin - j;
1638 az = cd.azimut(i, j);
1639 break;
1640 case 1:
1641 x = max_bin + j;
1642 y = max_bin + i;
1643 az = cd.azimut(i, j) + 90.;
1644 break;
1645 case 2:
1646 x = max_bin - i;
1647 y = max_bin + j;
1648 az = cd.azimut(i, j) + 180.;
1649 break;
1650 case 3:
1651 x = max_bin - j;
1652 y = max_bin - i;
1653 az = cd.azimut(i, j)+270.;
1654 break;
1655 }
1656
1657 az_min = (int)((az - .45)/.9);
1658 az_max = ceil((az + .45)/.9);
1659
1660
1661 if(az_min < 0)
1662 {
1663 az_min = az_min + NUM_AZ_X_PPI;
1664 az_max = az_max + NUM_AZ_X_PPI;
1665 }
1666 for(iaz = az_min; iaz<az_max; iaz++){
1667 // Enrico: cerca di non leggere fuori dal volume effettivo
1668 unsigned char sample = 0;
1669 if (irange < volume[el_index].beam_size)
1670 sample = max((unsigned char) (volume[el_index].get(iaz%NUM_AZ_X_PPI, irange)), (unsigned char)1); // il max serve perchè il valore di MISSING è 0
1671 if(cart(y, x) <= sample) cart(y, x) = sample;
1672 }
1673 }
1674}
1675
1676void SingleCart::write_out(Assets& assets, const std::string tagname, const std::string format)
1677{
1678 if (getenv("DIR_DEBUG") == NULL) return;
1679 assets.write_gdal_image(cart, "DIR_DEBUG", tagname.c_str(), format.c_str());
1680}
1681
1682}
1683
1684char *PrendiOra()
1685{
1686 time_t clock;
1687 struct tm *tempo;
1688
1689 clock=time(0);
1690
1691 tempo=gmtime(&clock);
1692
1693 return asctime(tempo);
1694}
1695
1696void prendo_tempo()
1697{
1698 static time_t time_tot = 0,time1 = 0,time2 = 0;
1699 LOG_CATEGORY("radar.timing");
1700
1701 if(time1 == 0){
1702 time1=time(&time1);
1703 time_tot = time1;
1704 }
1705 time2 = time(&time2);
1706
1707 LOG_INFO("tempo parziale %ld ---- totale %ld", time2-time1, time2-time_tot);
1708 time1=time2;
1709 return;
1710}
void load_first_level_bb_bloc(radarelab::Matrix2D< unsigned char > &matrix)
Open the first level elevation BB bloc file.
Definition assets.cpp:135
void write_vpr0(const radarelab::algo::VPR &vpr)
Write in $VPR0_FILE the vpr calculated.
Definition assets.cpp:437
bool read_0term(float &zeroterm)
Read $FILE_ZERO_TERMICO.
Definition assets.cpp:278
void write_vpr_hmax(int hvprmax)
write in $VPR_HMAX the vpr peak's height.
Definition assets.cpp:322
void load_dem(radarelab::Matrix2D< float > &matrix)
Open the dem file.
Definition assets.cpp:115
void configure(const Site &site, time_t acq_time)
Configure asset lookup with the given details.
Definition assets.cpp:41
float read_t_ground() const
fornisce temperatura al suolo, da lettura file esterno
Definition assets.cpp:201
void write_gdal_image(const radarelab::Matrix2D< T > &image, const char *dir_env_var, const char *name, const char *format)
Write a graphic image with gdal.
Definition assets.cpp:635
long int read_profile_gap() const
Read the gap between the time in $LAST_VPR and the current acquisition time.
Definition assets.cpp:231
bool find_vpr0(const radarelab::algo::DBZ &dbz, radarelab::algo::VPR &vpr0, long int &gap)
Read the gap and the vpr0, and if vpr0 is not found, look it up among the archived VPRs.
Definition assets.cpp:391
void load_first_level(radarelab::Matrix2D< unsigned char > &matrix)
Open the first level file.
Definition assets.cpp:120
void load_first_level_bb_el(radarelab::Matrix2D< unsigned char > &matrix)
Open the first level elevation BB el file.
Definition assets.cpp:130
int read_vpr_hmax()
Read in $VPR_HMAX the vpr peak's height.
Definition assets.cpp:305
void write_last_vpr()
Write the acquisition time in $LAST_VPR file.
Definition assets.cpp:294
Finds resources, like data files, used by the program.
Definition assets.h:38
bool do_declutter
use only static declutter map
Definition cum_bac.h:96
bool do_devel
Produce additional output.
Definition cum_bac.h:99
void declutter_anaprop()
funzione che elabora il dato radar rimuovendo anaprop e beam blocking
Definition cum_bac.cpp:389
bool do_readStaticMap
Read Static clutter map.
Definition cum_bac.h:100
radarelab::algo::Anaprop< double > anaprop
Oggetto per correzione ANAPRO.
Definition cum_bac.h:129
radarelab::PolarScan< float > dem
dem in coordinate azimut range
Definition cum_bac.h:132
CalcoloVPR * calcolo_vpr
Oggetto per calcolare e correggere con VPR.
Definition cum_bac.h:113
radarelab::Volume< unsigned char > qual
qualita volume polare
Definition cum_bac.h:135
void want_vpr()
Call this just after creating the CUM_BAC object, to signal that VPR should also be computed.
Definition cum_bac.cpp:147
radarelab::PolarScan< unsigned char > first_level_static
mappa statica
Definition cum_bac.h:124
void leggo_first_level()
funzione che legge la mappa statica e la mappa di elevazioni da beam blocking e le condensa in un uni...
Definition cum_bac.cpp:503
log4c_category_t * logging_category
logging category
Definition cum_bac.h:83
bool do_class
Convective-stratiform classification.
Definition cum_bac.h:97
radarelab::Volume< double > & volume
Set to Z undetect value the Zpixels classified as non-meteo echoes.
Definition cum_bac.h:106
radarelab::PolarScan< unsigned char > bb_first_level
mappa di elevazioni da beam blocking (input)
Definition cum_bac.h:126
static void read_sp20_volume(radarelab::Volume< double > &volume, const Site &site, const char *nome_file, bool do_clean=false, bool do_medium=false)
Read from SP20 data file.
Definition cum_bac.cpp:152
const Site & site
site information object
Definition cum_bac.h:87
radarelab::PolarScan< unsigned char > first_level
mappa dinamica complessiva
Definition cum_bac.h:123
char date[20]
Acquisition date.
Definition cum_bac.h:120
radarelab::Volume< double > SD_Z6
Polar volume of standard deviation of reflectivity over 6 km length.
Definition cum_bac.h:107
CUM_BAC(radarelab::Volume< double > &volume, const Config &cfg, const Site &site, bool medium=false, unsigned max_bin=512)
Constructor.
Definition cum_bac.cpp:119
bool do_anaprop
anaprop correction
Definition cum_bac.h:101
bool do_bloccorr
bloccorrection
Definition cum_bac.h:95
static void read_odim_volume(radarelab::Volume< double > &volume, const Site &site, const char *nome_file, char *fuzzypath, bool do_clean=false, bool do_medium=false, bool set_undetect=false)
Read from ODIM data file.
Definition cum_bac.cpp:190
void vpr_class()
Esegue tutta la catena vpr (e classificazione) se richiesta.
Definition cum_bac.cpp:1383
Assets assets
others
Definition cum_bac.h:88
bool do_medium
medium processing flag
Definition cum_bac.h:90
radarelab::algo::DBZ dbz
????
Definition cum_bac.h:110
bool do_quality
Feature set required for this run.
Definition cum_bac.h:93
radarelab::CylindricalVolume cil
Volume resampled as a cylindrical volume.
Definition cum_bac.h:108
radarelab::PolarScan< unsigned char > top
Echo top a ???? dBZ [hm].
Definition cum_bac.h:137
bool do_beamblocking
beamblocking corretion
Definition cum_bac.h:94
void ScrivoStatistica(const radarelab::algo::anaprop::GridStats &)
funzione scrittura matrici statistica
Definition cum_bac.cpp:570
radarelab::PolarScan< unsigned char > beam_blocking
mappa di beam blocking (input)
Definition cum_bac.h:127
bool do_zlr_media
Compute ZLR map using averaging.
Definition cum_bac.h:98
void caratterizzo_volume()
funzione che caratterizza i volumi polari tramite la qualita'
Definition cum_bac.cpp:636
unsigned MyMAX_BIN
maximum number of beam size
Definition cum_bac.h:85
void conversione_convettiva()
Nei punti convettivi ricalcola la Z come MP classica, partendo dalla stima di R (convettiva)
Definition cum_bac.cpp:1375
Classe principale del programma.
Definition cum_bac.h:62
bool open_from_env(const char *varname, const char *mode, const char *desc=nullptr)
Opens a file taking its name from the environment variable envname.
Definition utils.cpp:37
bool open(const std::string &fname, const char *mode, const char *desc=nullptr)
Opens a file by its pathname.
Definition utils.cpp:49
Open a file taking its name from a given env variable.
Definition utils.h:22
void to_cart(const PolarScan< SRC > &src, Matrix2D< DST > &dst) const
Copy data from the polar scan src to the cartesian map dst.
Definition cart.h:82
void resize_beams_and_propagate_last_bin(unsigned new_beam_size)
Enlarges the PolarScan increasing beam_size and propagating the last bin value.
Definition volume.h:212
PolarScan - structure to describe a polarScan containing a matrix of data and conversion factors.
Definition volume.h:113
const unsigned max_beam_size() const
Return the maximum beam size in all PolarScans.
Definition volume.h:457
const unsigned beam_count
Number of beam_count used ast each elevations.
Definition volume.h:434
Homogeneous volume with a common beam count for all PolarScans.
Definition volume.h:431
static unsigned char DBtoBYTE(double DB, double gain=80./255., double offset=-20.)
funzione che converte dB in valore intero tra 0 e 255
Definition dbz.h:94
static constexpr double ZtoDBZ(double Z)
funzione che converte Z in dBZ
Definition dbz.h:126
static constexpr double beam_blocking_correction(double val_db, double beamblocking)
@function Compute the corrected value (in dB) given the original value (in dB) and the beam blocking ...
Definition dbz.h:70
double attenuation(unsigned char DBZbyte, double PIA)
funzione che calcola l'attenuazione totale
Definition dbz.cpp:51
Class to manage reflectivity functions (simply attenuation correction, conversion between Z,...
Definition dbz.h:23
RadarSite radarSite
RadarSite.
Definition volume.h:274
T offset
Data Offset.
Definition volume.h:276
void normalize_elevations(const std::vector< double > &elevations)
Change the elevations in the PolarScans to match the given elevation vector.
Definition volume.h:397
PolarScan< T > & append_scan(unsigned beam_count, unsigned beam_size, double elevation, double cell_size, const T &default_value=algo::DBZ::BYTEtoDB(1))
Append a scan to this volume.
Definition volume.h:332
std::shared_ptr< LoadInfo > load_info
Polar volume information.
Definition volume.h:272
PolarScan< T > & scan(unsigned idx)
Access a polar scan.
Definition volume.h:313
Sequence of PolarScans which can have a different beam count for each elevation.
Definition volume.h:264
codice principale di elaborazione dei volumi di riflettivita' radar usato per impulso corto
float func_q_Z(unsigned char cl, unsigned char bb, float dst, float dr, float dt, float dh, float dhst, float PIA)
funzione che calcola la qualita' per Z
Definition func_Q3d.cpp:8
funzioni che combinano le componenti semplici di qualita' radar
name space generale del programma
Definition assets.h:28
Namespace per volume dati.
Definition elev_fin.h:12
String functions.
Definition cart.cpp:4
Codice per il caricamento di volumi ODIM in radarelab.
settaggio ambiente lavoro nel caso non sia settato dall'esterno
definisce struttura Site Contiene le informazioni di base che caratterizzano il sito radar
void classifica_rain()
funzione che classifica la precipitazione se stratiforme o convettiva
Definition cum_bac.cpp:748
void esegui_tutto()
Metodo che lancia tutte le operazioni per il calcolo e la correzione del vpr.
Definition cum_bac.cpp:1538
int ier_stampa_vpr
flag d'errore su stampa profilo
Definition cum_bac.h:249
radarelab::PolarScan< unsigned char > neve
matrice az-range che memorizza punti di neve
Definition cum_bac.h:247
unsigned MyMAX_BIN
LUNGHEZZA MASSIMA.
Definition cum_bac.h:251
int heating
variabile di riscaldamento
Definition cum_bac.h:242
int ier_max
flag d'errore su calcolo quota max
Definition cum_bac.h:248
int corr_vpr()
correzione vpr
Definition cum_bac.cpp:975
int combina_profili(const radarelab::algo::InstantaneousVPR &inst_vpr)
funzione che combina il profilo verticale corrente con quello precedente tramite il metodo di Germann
Definition cum_bac.cpp:847
long int gap
distanza temporale dall'ultimo file vpr [numero acquisizioni intercorse dall'ultimo vpr ?...
Definition cum_bac.h:232
void merge_metodi(const radarelab::algo::CalcoloSteiner &steiner, const radarelab::algo::CalcoloVIZ &viz)
fa il merge dei metodi
Definition cum_bac.cpp:817
int stampa_vpr()
stampa profilo combinato
Definition cum_bac.cpp:935
int livmin
quota livello minimo calcolato
Definition cum_bac.h:243
radarelab::algo::VPR vpr
Informa se il pixel è convettivo.
Definition cum_bac.h:236
int analyse_VPR(float *vpr_liq, int *snow, float *hliq)
funzione che analizza il profilo
Definition cum_bac.cpp:1180
int hvprmax
quota picco vpr
Definition cum_bac.h:237
double hbbb
altezza bottom brightband
Definition cum_bac.h:245
int trovo_hvprmax(int *hmax)
trova il massimo del profilo
Definition cum_bac.cpp:1100
CUM_BAC & cum_bac
oggeto CUM_BAC di riferimento
Definition cum_bac.h:230
float t_ground
2m temperature
Definition cum_bac.h:233
int profile_heating(bool has_inst_vpr)
calcola riscaldamento in quarti d'ora
Definition cum_bac.cpp:905
double htbb
altezza top brightband
Definition cum_bac.h:244
CalcoloVPR(CUM_BAC &cum_bac)
Constructor.
Definition cum_bac.cpp:1512
radarelab::PolarScan< unsigned char > corr_polar
correzione vpr in byte 0-128 negativa 128-256 positiva, in coord az-ra
Definition cum_bac.h:246
Struttara per il calcolo del VPR.
Definition cum_bac.h:227
const unsigned max_bin
dimensione matrice
Definition cum_bac.h:337
radarelab::Image< unsigned char > cart
vol_pol interpolated in a cartesian map
Definition cum_bac.h:340
void write_out(Assets &assets, const std::string tagname, const std::string format="PNG")
Produce in output le immagini PNG dei campi in $DIR_DEBUG.
Definition cum_bac.cpp:1676
SingleCart(unsigned max_bin)
Constructor.
Definition cum_bac.cpp:1604
void creo_cart(const radarelab::Volume< double > &volume, unsigned int el_index)
conversione da polare a cartesiano alta risoluzione
Definition cum_bac.cpp:1610
std::string name
Nome sito radar.
Definition site.h:29
virtual unsigned char get_bin_wind_magic_number(time_t when) const =0
Return the magic number for wind to be used in clean procedure.
RadarSite radarSite
Description of radar site.
Definition site.h:41
virtual std::vector< double > get_elev_array(bool medium=false) const =0
return the elev array used
Radar site information.
Definition site.h:24
std::vector< Matrix2D< double > * > slices
Vertical rectangular x,z semi-slices of the cylinder, with one side resting on the cylinder axis.
Definition cylindrical.h:26
Radar volume mapped to cylindrical coordinates.
Definition cylindrical.h:17
Base for all matrices we use, since we rely on row-major data.
Definition matrix.h:37
double sample_height(unsigned cell_idx) const
Return the height (in meters) of the sample at the given cell indexa.
Definition volume.cpp:37
void to_cart_average(const PolarScan< double > &src, std::function< T(const std::vector< double > &)> &convert, Matrix2D< T > &dst) const
Fill the cartesian map dst with the output of the function src(azimuth, range)
Definition cart.h:186
void resample_volume(const Volume< T > &src, Volume< T > &dst, double src_beam_width) const
Merge.
std::map< std::string, Scans< double > * > to_load
Map used to specify quantity to be loaded.
Definition loader.h:26
void request_quantity(const std::string &name, Scans< double > *volume)
Define a request - Fill to_load attribute
Definition odim.cpp:29
void load(const std::string &pathname)
Load method.
Definition odim.cpp:34
Struttura che eredita da Loader e definisce i metodi per accedere ai dati ODIM.
Definition odim.h:23