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