Elaboradar  0.1
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
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 {
505  if (do_readStaticMap)
506  {
507  // Leggo mappa statica
509  // Allargo per coprire la dimensione del volume
510  if (first_level_static.cols() < volume.max_beam_size())
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;
544  first_level(i, j)=first_level_static(i, j);
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 
570 void 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
634 comend
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
761  gap = cum_bac.assets.read_profile_gap();
762  //-- se gap < memory leggo hmax da VPR
763  if (gap<=MEMORY){
764  hmax = cum_bac.assets.read_vpr_hmax();
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 
817 void 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 */
847 int 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-----
899  cum_bac.assets.write_vpr0(vpr);
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 
905 int 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 */
1180 int 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 
1369  cum_bac.assets.write_vpr_hmax(hvprmax);
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 
1412 void 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 
1534 CalcoloVPR::~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 
1585 namespace {
1586 struct 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 
1604 SingleCart::SingleCart(unsigned max_bin)
1605  : max_bin(max_bin),
1606  cart(max_bin*2)
1607 {
1608 }
1609 
1610 void 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 
1676 void 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 
1684 char *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 
1696 void 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 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
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
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
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 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
const unsigned max_beam_size() const
Return the maximum beam size in all PolarScans.
Definition: volume.h:455
const unsigned beam_count
Number of beam_count used ast each elevations.
Definition: volume.h:432
Homogeneous volume with a common beam count for all PolarScans.
Definition: volume.h:429
double attenuation(unsigned char DBZbyte, double PIA)
funzione che calcola l'attenuazione totale
Definition: dbz.cpp:51
RadarSite radarSite
RadarSite.
Definition: volume.h:272
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:330
T offset
Data Offset.
Definition: volume.h:274
PolarScan< T > & scan(unsigned idx)
Access a polar scan.
Definition: volume.h:311
void normalize_elevations(const std::vector< double > &elevations)
Change the elevations in the PolarScans to match the given elevation vector.
Definition: volume.h:395
std::shared_ptr< LoadInfo > load_info
Polar volume information.
Definition: volume.h:270
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 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
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 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 std::vector< double > get_elev_array(bool medium=false) const =0
return the elev array used
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
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
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