Elaboradar  0.1
 Tutto Classi Namespace File Funzioni Variabili Tipi enumerati (enum) Gruppi
volume.h
Vai alla documentazione di questo file.
1 
6 #ifndef RADARELAB_VOLUME_H
7 #define RADARELAB_VOLUME_H
8 
9 #include <radarelab/logging.h>
10 #include <radarelab/matrix.h>
11 #include <radarelab/algo/dbz.h>
12 #include <string>
13 #include <vector>
14 #include <cstdio>
15 #include <cmath>
16 #include <algorithm>
17 #include <memory>
18 #include <Eigen/Core>
19 #include <radarelab/RadarSite.h>
20 
21 namespace radarelab {
22 
28 {
30  unsigned beam_count = 0;
31 
33  unsigned beam_size = 0;
34 
36  Eigen::VectorXd azimuths_real;
37 
42  double elevation = 0;
43 
45  Eigen::VectorXd elevations_real;
46 
48  double cell_size = 0;
49 
50  /*
51  * Constructor
52  * Create a PolarScanBase defining beam_count and beam_size
53  * @param [in] beam_count
54  * @param [in] beam_size
55  */
56  PolarScanBase(unsigned beam_count, unsigned beam_size);
57 
58  /*
59  * Constructor
60  * Create a copy of a PolarScanBase
61  * @param [in] s - PolarScanBase to be copied
62  */
63  PolarScanBase(const PolarScanBase& s);
64 
70  double height(unsigned rg, double beam_half_width=0.0);
71 
78  double diff_height(unsigned rg_start, unsigned rg_end);
79 
85  double sample_height(unsigned cell_idx) const;
86 
95  static double sample_height(double elevation, double range, double equiv_earth_radius);
96 
105  static double sample_height(double elevation, double range);
106 };
107 
111 template<typename T>
112 class PolarScan : public PolarScanBase, public Matrix2D<T>
113 {
114 public:
116  double nodata = 0;
118  double undetect = 0;
120  double gain = 1;
122  double offset = 0;
123 
130  PolarScan(unsigned beam_count, unsigned beam_size, const T& default_value=algo::DBZ::BYTEtoDB(1))
131  : PolarScanBase(beam_count, beam_size),
132  Matrix2D<T>(PolarScan::Constant(beam_count, beam_size, default_value))
133  {
134  }
135 
141  // FIXME: Cannot use '= default' or g++ 4.8.3 will not instantiate it on explicit template instantiation
143  : PolarScanBase(s), Matrix2D<T>(s), nodata(s.nodata),
144  undetect(s.undetect), gain(s.gain), offset(s.offset) {}
145 
146  template<class OT>
147  PolarScan(const PolarScan<OT>& s, const T& default_value)
148  : PolarScanBase(s), Matrix2D<T>(PolarScan::Constant(s.beam_count, s.beam_size, default_value)),
150  {
151  }
152 
153  ~PolarScan()
154  {
155  }
156 
161  T get(unsigned az, unsigned beam) const
162  {
163  return (*this)(az, beam);
164  }
165 
170  void set(unsigned az, unsigned beam, T val)
171  {
172  (*this)(az, beam) = val;
173  }
174 
182  double sample_height_real(unsigned az, unsigned cell_idx) const
183  {
184  return PolarScanBase::sample_height(elevations_real(az), (cell_idx + 0.5) * cell_size);
185  }
186 
194  void read_beam(unsigned az, T* out, unsigned out_size, T missing=0) const
195  {
196  using namespace std;
197 
198  // Prima riempio il minimo tra ray.size() e out_size
199  size_t set_count = min(beam_size, out_size);
200 
201  for (unsigned i = 0; i < set_count; ++i)
202  out[i] = get(az, i);
203 
204  for (unsigned i = set_count; i < out_size; ++i)
205  out[i] = missing;
206  }
207 
212  void resize_beams_and_propagate_last_bin(unsigned new_beam_size)
213  {
214  if (new_beam_size <= beam_size) return;
215  this->conservativeResize(Eigen::NoChange, new_beam_size);
216  this->rightCols(new_beam_size - this->beam_size).colwise() = this->col(this->beam_size - 1);
217  this->beam_size = new_beam_size;
218  }
219 };
220 
221 
222 struct VolumeStats
223 {
224  std::vector<unsigned> count_zeros;
225  std::vector<unsigned> count_ones;
226  std::vector<unsigned> count_others;
227  std::vector<unsigned> sum_others;
228 
229  void print(FILE* out);
230 };
231 
235 namespace volume {
236 
240 struct LoadInfo
241 {
243  std::string filename;
244  // std::string date;
245  // std::string time;
247  time_t acq_date;
250 
251  LoadInfo()
252  : declutter_rsp(false)
253  {
254  }
255 };
256 
260 template<typename T>
261 class Scans : public std::vector<PolarScan<T>>
262 {
263 public:
264  typedef T Scalar;
266  std::string quantity;
268  std::string units;
270  std::shared_ptr<LoadInfo> load_info;
272  RadarSite radarSite;
273 
274  // FIXME: Cannot use '= default' or g++ 4.8.3 will not instantiate it on explicit template instantiation
275  Scans() : std::vector<PolarScan<T>>() {}
276 
283  template<typename OT>
284  Scans(const Scans<OT>& v, const T& default_value)
285  {
286  this->quantity = v.quantity;
287  this->units = v.units;
288  this->load_info = v.load_info;
289  this->reserve(v.size());
290  this->radarSite = v.radarSite;
291  for (const auto& src_scan : v)
292  this->push_back(PolarScan<T>(src_scan, default_value));
293  }
294 
298  PolarScan<T>& scan(unsigned idx) { return (*this)[idx]; }
302  const PolarScan<T>& scan(unsigned idx) const { return (*this)[idx]; }
303 
317  PolarScan<T>& append_scan(unsigned beam_count, unsigned beam_size, double elevation, double cell_size)
318  {
319  // Ensure elevations grow as scan indices grow
320  if (!this->empty() && elevation <= this->back().elevation)
321  {
322  LOG_CATEGORY("radar.io");
323  LOG_ERROR("append_scan(beam_count=%u, beam_size=%u, elevation=%f, cell_size=%f) called with an elevation that is not above the last one (%f)", beam_count, beam_size, elevation, cell_size, this->back().elevation);
324  throw std::runtime_error("elevation not greather than the last one");
325  }
326  // Add the new polar scan
327  this->push_back(PolarScan<T>(beam_count, beam_size));
328  this->back().elevation = elevation;
329  this->back().cell_size = cell_size;
330  return this->back();
331  }
332 
342  PolarScan<T>& make_scan(unsigned idx, unsigned beam_count, unsigned beam_size, double elevation, double cell_size)
343  {
344  if (idx < this->size())
345  {
346  if (beam_count != (*this)[idx].beam_count)
347  {
348  LOG_CATEGORY("radar.io");
349  LOG_ERROR("make_scan(idx=%u, beam_count=%u, beam_size=%u) called, but the scan already existed with beam_count=%u", idx, beam_count, beam_size, (*this)[idx].beam_count);
350  throw std::runtime_error("beam_size mismatch");
351  }
352  if (beam_size != (*this)[idx].beam_size)
353  {
354  LOG_CATEGORY("radar.io");
355  LOG_ERROR("make_scan(idx=%u, beam_count=%u, beam_size=%u) called, but the scan already existed with beam_size=%u", idx, beam_count, beam_size, (*this)[idx].beam_size);
356  throw std::runtime_error("beam_size mismatch");
357  }
358  } else {
359  // If some elevation has been skipped, fill in the gap
360  if (idx > this->size())
361  {
362  if (this->empty())
363  this->push_back(PolarScan<T>(beam_count, beam_size));
364  while (this->size() < idx)
365  this->push_back(PolarScan<T>(beam_count, this->back().beam_size));
366  }
367 
368  // Add the new polar scan
369  this->push_back(PolarScan<T>(beam_count, beam_size));
370  this->back().elevation = elevation;
371  this->back().cell_size = cell_size;
372  }
373 
374  // Return it
375  return (*this)[idx];
376  }
377 
382  void normalize_elevations(const std::vector<double>& elevations)
383  {
384  // Ensure that we have enough standard elevations
385  if (elevations.size() < this->size())
386  {
387  LOG_CATEGORY("radar.io");
388  LOG_ERROR("normalize_elevations: standard elevation array has %zd elements, but we have %zd scans",
389  elevations.size(), this->size());
390  throw std::runtime_error("not enough standard elevations");
391  }
392  // Ensure that the nudging that we do do not confuse a scan
393  // with another
394  for (size_t i = 0; i < this->size() - 1; ++i)
395  {
396  if (abs(elevations[i] - this->at(i).elevation) > abs(elevations[i] - this->at(i + 1).elevation))
397  {
398  LOG_CATEGORY("radar.io");
399  LOG_ERROR("normalize_elevations: elevation %zd (%f) should be set to %f but it would make it closer to the next elevation %f", i, this->at(i).elevation, elevations[i], this->at(i + 1).elevation);
400  throw std::runtime_error("real elevation is too different than standard elevations");
401  }
402  }
403  // Assign the new elevations
404  for (size_t i = 0; i < this->size(); ++i)
405  this->at(i).elevation = elevations[i];
406  }
407 };
408 
409 }
410 
414 template<typename T>
415 class Volume : public volume::Scans<T>
416 {
417 public:
419  const unsigned beam_count;
420 
425  Volume(unsigned beam_count)
426  : beam_count(beam_count)
427  {
428  }
429 
435  template<typename OT>
436  Volume(const Volume<OT>& v, const T& default_value)
437  : volume::Scans<T>(v, default_value), beam_count(v.beam_count)
438  {
439  }
440 
442  const unsigned max_beam_size() const
443  {
444  unsigned res = 0;
445  for (const auto& scan : *this)
446  res = std::max(res, scan.beam_size);
447  return res;
448  }
449 
451  bool is_unique_cell_size() const
452  {
453  double cell_size = this->at(0).cell_size;
454  for (size_t i = 1; i < this->size(); ++i)
455  if ( this->at(i).cell_size != cell_size) return false;
456  return true;
457  }
458 
460  double elevation_min() const
461  {
462  return this->front().elevation;
463  }
464 
466  double elevation_max() const
467  {
468  return this->back().elevation;
469  }
470 
477  void read_vertical_slice(unsigned az, Matrix2D<T>& slice, double missing_value) const
478  {
479  unsigned size_z = std::max(this->size(), (size_t)slice.rows());
480  for (unsigned el = 0; el < size_z; ++el)
481  this->scan(el).read_beam(az, slice.row_ptr(el), slice.cols(), missing_value);
482  }
483 
485  void compute_stats(VolumeStats& stats) const
486  {
487  stats.count_zeros.resize(this->size());
488  stats.count_ones.resize(this->size());
489  stats.count_others.resize(this->size());
490  stats.sum_others.resize(this->size());
491 
492  for (unsigned iel = 0; iel < this->size(); ++iel)
493  {
494  stats.count_zeros[iel] = 0;
495  stats.count_ones[iel] = 0;
496  stats.count_others[iel] = 0;
497  stats.sum_others[iel] = 0;
498 
499  for (unsigned iaz = 0; iaz < this->scan(iel).beam_count; ++iaz)
500  {
501  for (size_t i = 0; i < this->scan(iel).beam_size; ++i)
502  {
503  int val = algo::DBZ::DBtoBYTE(this->scan(iel).get(iaz, i));
504  switch (val)
505  {
506  case 0: stats.count_zeros[iel]++; break;
507  case 1: stats.count_ones[iel]++; break;
508  default:
509  stats.count_others[iel]++;
510  stats.sum_others[iel] += val;
511  break;
512  }
513  }
514  }
515  }
516  }
517 
530  {
531  return volume::Scans<T>::append_scan(beam_count, beam_size, elevation, cell_size);
532  }
533 
542  PolarScan<T>& make_scan(unsigned idx, unsigned beam_size, double elevation, double cell_size)
543  {
544  return volume::Scans<T>::make_scan(idx, beam_count, beam_size, elevation, cell_size);
545  }
546 
551  Volume& operator*=(const T coefficient)
552  {
553  for(unsigned el=0;el<this->size();el++)
554  {
555  this->scan(el)*=coefficient;
556  }
557  return *this;
558  }
559 
565  {
566  for (unsigned el=0;el<this->size();el++)
567  {
568  this->scan(el)+=addend[el];
569  }
570  return *this;
571  }
572 
573 protected:
574  void resize_elev_fin();
575 
576 private:
577  Volume& operator=(Volume&);
578  Volume(const Volume&);
579 };
580 
581 /*
582  * Tell the compiler that the implementation of these templates is already
583  * explicitly instantiated in volume.cpp, and it should not spend time
584  * reinstantiating it every time the template is used.
585  */
586 extern template class PolarScan<double>;
587 extern template class Volume<double>;
588 namespace volume {
589 extern template class Scans<double>;
590 }
591 
592 }
593 
594 #endif
double gain
Conversion factor.
Definition: volume.h:120
double undetect
Minimum amount that can be measured.
Definition: volume.h:118
Eigen::VectorXd elevations_real
Vector of actual elevations for each beam.
Definition: volume.h:45
PolarScan< T > & scan(unsigned idx)
Access a polar scan.
Definition: volume.h:298
static constexpr double BYTEtoDB(unsigned char DBZbyte, double gain=80./255., double offset=-20.)
funzione che converte Z unsigned char in DBZ
Definition: dbz.h:82
void read_beam(unsigned az, T *out, unsigned out_size, T missing=0) const
Fill an array with beam data .
Definition: volume.h:194
PolarScan(unsigned beam_count, unsigned beam_size, const T &default_value=algo::DBZ::BYTEtoDB(1))
Definition: volume.h:130
void compute_stats(VolumeStats &stats) const
Compute Volume statistics.
Definition: volume.h:485
Volume & operator*=(const T coefficient)
*= operator defined
Definition: volume.h:551
PolarScan - structure to describe a polarScan containing a matrix of data and conversion factors...
Definition: volume.h:112
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
double cell_size
Size of a beam cell in meters.
Definition: volume.h:48
unsigned beam_count
Count of beams in this scan.
Definition: volume.h:30
PolarScan< T > & append_scan(unsigned beam_size, double elevation, double cell_size)
Append a scan to this volume.
Definition: volume.h:529
std::string filename
Original file name.
Definition: volume.h:243
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:419
double offset
Conversion factor.
Definition: volume.h:122
LoadInfo structure - Contains generic volume information.
Definition: volume.h:240
std::shared_ptr< LoadInfo > load_info
Polar volume information.
Definition: volume.h:270
Volume(const Volume< OT > &v, const T &default_value)
Copy constructor.
Definition: volume.h:436
PolarScan< T > & append_scan(unsigned beam_count, unsigned beam_size, double elevation, double cell_size)
Append a scan to this volume.
Definition: volume.h:317
void normalize_elevations(const std::vector< double > &elevations)
Change the elevations in the PolarScans to match the given elevation vector.
Definition: volume.h:382
unsigned beam_size
Number of samples in each beam.
Definition: volume.h:33
void read_vertical_slice(unsigned az, Matrix2D< T > &slice, double missing_value) const
Fill a matrix elevations x beam_size with the vertical slice at a given azimuth.
Definition: volume.h:477
Basic structure to describe a polar scan, independently of the type of its samples.
Definition: volume.h:27
double diff_height(unsigned rg_start, unsigned rg_end)
Height difference in kilometers (legacy) between two range gates.
Definition: volume.cpp:32
double elevation_min() const
Return the lowest elevation.
Definition: volume.h:460
Volume(unsigned beam_count)
Constructor.
Definition: volume.h:425
std::string units
Data units according to ODIM documentation.
Definition: volume.h:268
bool is_unique_cell_size() const
Test if same cell_size in all PolarScans.
Definition: volume.h:451
Homogeneous volume with a common beam count for all PolarScans.
Definition: volume.h:415
PolarScan(const PolarScan &s)
Constructor Create a copy of a PolarScan.
Definition: volume.h:142
Sequence of PolarScans which can have a different beam count for each elevation.
Definition: volume.h:261
void set(unsigned az, unsigned beam, T val)
Set a beam value.
Definition: volume.h:170
bool declutter_rsp
flag true if data have been decluttered with Doppler at rsp level
Definition: volume.h:249
double sample_height(unsigned cell_idx) const
Return the height (in meters) of the sample at the given cell indexa.
Definition: volume.cpp:37
static unsigned char DBtoBYTE(double DB, double gain=80./255., double offset=-20.)
funzione che converte dB in valore intero tra 0 e 255
Definition: dbz.h:94
RadarSite radarSite
RadarSite.
Definition: volume.h:272
Scans(const Scans< OT > &v, const T &default_value)
Constructor Copy from another Scans.
Definition: volume.h:284
const unsigned max_beam_size() const
Return the maximum beam size in all PolarScans.
Definition: volume.h:442
PolarScan< T > & make_scan(unsigned idx, unsigned beam_count, unsigned beam_size, double elevation, double cell_size)
Create or reuse a scan at position idx, with the given beam size.
Definition: volume.h:342
std::string quantity
Odim quantity name.
Definition: volume.h:266
Eigen::VectorXd azimuths_real
Vector of actual azimuths for each beam.
Definition: volume.h:36
double sample_height_real(unsigned az, unsigned cell_idx) const
Return the height (in meters) of a sample given its azimuth and cell indices use the real beam elevat...
Definition: volume.h:182
double nodata
Value used as &#39;no data&#39; value.
Definition: volume.h:116
double height(unsigned rg, double beam_half_width=0.0)
Height in kilometers (legacy) at range gate for beam elevation + beam_half_width. ...
Definition: volume.cpp:27
const PolarScan< T > & scan(unsigned idx) const
Access a polar scan (const)
Definition: volume.h:302
double elevation_max() const
Return the highest elevation.
Definition: volume.h:466
PolarScan< T > & make_scan(unsigned idx, unsigned beam_size, double elevation, double cell_size)
Create or reuse a scan at position idx, with the given beam size.
Definition: volume.h:542
time_t acq_date
Acquisition date.
Definition: volume.h:247
Volume & operator+=(Volume &addend)
+= operator defined
Definition: volume.h:564
double elevation
Nominal elevation of this PolarScan, which may be different from the effective elevation of each sing...
Definition: volume.h:42
Radar Site description.