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;
274  T offset = 0;
275 
276  // FIXME: Cannot use '= default' or g++ 4.8.3 will not instantiate it on explicit template instantiation
277  Scans() : std::vector<PolarScan<T>>() {}
278 
285  template<typename OT>
286  Scans(const Scans<OT>& v, const T& default_value)
287  {
288  this->quantity = v.quantity;
289  this->units = v.units;
290  this->load_info = v.load_info;
291  this->reserve(v.size());
292  this->radarSite = v.radarSite;
293  for (const auto& src_scan : v)
294  this->push_back(PolarScan<T>(src_scan, default_value));
295  this->offset=v.offset;
296  }
297 
300  //
301  void SetOffset (T offset)
302  {
303  this->offset = offset;
304  }
305 
306 
307 
311  PolarScan<T>& scan(unsigned idx) { return (*this)[idx]; }
315  const PolarScan<T>& scan(unsigned idx) const { return (*this)[idx]; }
316 
330  PolarScan<T>& append_scan(unsigned beam_count, unsigned beam_size, double elevation, double cell_size)
331  {
332  // Ensure elevations grow as scan indices grow
333  if (!this->empty() && elevation <= this->back().elevation)
334  {
335  LOG_CATEGORY("radar.io");
336  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);
337  throw std::runtime_error("elevation not greather than the last one");
338  }
339  // Add the new polar scan
340  this->push_back(PolarScan<T>(beam_count, beam_size));
341  this->back().elevation = elevation;
342  this->back().cell_size = cell_size;
343  return this->back();
344  }
345 
355  PolarScan<T>& make_scan(unsigned idx, unsigned beam_count, unsigned beam_size, double elevation, double cell_size)
356  {
357  if (idx < this->size())
358  {
359  if (beam_count != (*this)[idx].beam_count)
360  {
361  LOG_CATEGORY("radar.io");
362  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);
363  throw std::runtime_error("beam_size mismatch");
364  }
365  if (beam_size != (*this)[idx].beam_size)
366  {
367  LOG_CATEGORY("radar.io");
368  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);
369  throw std::runtime_error("beam_size mismatch");
370  }
371  } else {
372  // If some elevation has been skipped, fill in the gap
373  if (idx > this->size())
374  {
375  if (this->empty())
376  this->push_back(PolarScan<T>(beam_count, beam_size));
377  while (this->size() < idx)
378  this->push_back(PolarScan<T>(beam_count, this->back().beam_size));
379  }
380 
381  // Add the new polar scan
382  this->push_back(PolarScan<T>(beam_count, beam_size));
383  this->back().elevation = elevation;
384  this->back().cell_size = cell_size;
385  }
386 
387  // Return it
388  return (*this)[idx];
389  }
390 
395  void normalize_elevations(const std::vector<double>& elevations)
396  {
397  // Ensure that we have enough standard elevations
398  if (elevations.size() < this->size())
399  {
400  LOG_CATEGORY("radar.io");
401  LOG_ERROR("normalize_elevations: standard elevation array has %zd elements, but we have %zd scans",
402  elevations.size(), this->size());
403  throw std::runtime_error("not enough standard elevations");
404  }
405  // Ensure that the nudging that we do do not confuse a scan
406  // with another
407  for (size_t i = 0; i < this->size() - 1; ++i)
408  {
409  if (abs(elevations[i] - this->at(i).elevation) > abs(elevations[i] - this->at(i + 1).elevation))
410  {
411  LOG_CATEGORY("radar.io");
412  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);
413  throw std::runtime_error("real elevation is too different than standard elevations");
414  }
415  }
416  // Assign the new elevations
417  for (size_t i = 0; i < this->size(); ++i)
418  this->at(i).elevation = elevations[i];
419  }
420 };
421 
422 }
423 
427 template<typename T>
428 class Volume : public volume::Scans<T>
429 {
430 public:
432  const unsigned beam_count;
433 
438  Volume(unsigned beam_count)
439  : beam_count(beam_count)
440  {
441  }
442 
448  template<typename OT>
449  Volume(const Volume<OT>& v, const T& default_value)
450  : volume::Scans<T>(v, default_value), beam_count(v.beam_count)
451  {
452  }
453 
455  const unsigned max_beam_size() const
456  {
457  unsigned res = 0;
458  for (const auto& scan : *this)
459  res = std::max(res, scan.beam_size);
460  return res;
461  }
462 
464  bool is_unique_cell_size() const
465  {
466  double cell_size = this->at(0).cell_size;
467  for (size_t i = 1; i < this->size(); ++i)
468  if ( this->at(i).cell_size != cell_size) return false;
469  return true;
470  }
471 
473  double elevation_min() const
474  {
475  return this->front().elevation;
476  }
477 
479  double elevation_max() const
480  {
481  return this->back().elevation;
482  }
483 
490  void read_vertical_slice(unsigned az, Matrix2D<T>& slice, double missing_value) const
491  {
492  unsigned size_z = std::max(this->size(), (size_t)slice.rows());
493  for (unsigned el = 0; el < size_z; ++el)
494  this->scan(el).read_beam(az, slice.row_ptr(el), slice.cols(), missing_value);
495  }
496 
498  void compute_stats(VolumeStats& stats) const
499  {
500  stats.count_zeros.resize(this->size());
501  stats.count_ones.resize(this->size());
502  stats.count_others.resize(this->size());
503  stats.sum_others.resize(this->size());
504 
505  for (unsigned iel = 0; iel < this->size(); ++iel)
506  {
507  stats.count_zeros[iel] = 0;
508  stats.count_ones[iel] = 0;
509  stats.count_others[iel] = 0;
510  stats.sum_others[iel] = 0;
511 
512  for (unsigned iaz = 0; iaz < this->scan(iel).beam_count; ++iaz)
513  {
514  for (size_t i = 0; i < this->scan(iel).beam_size; ++i)
515  {
516  int val = algo::DBZ::DBtoBYTE(this->scan(iel).get(iaz, i));
517  switch (val)
518  {
519  case 0: stats.count_zeros[iel]++; break;
520  case 1: stats.count_ones[iel]++; break;
521  default:
522  stats.count_others[iel]++;
523  stats.sum_others[iel] += val;
524  break;
525  }
526  }
527  }
528  }
529  }
530 
543  {
544  return volume::Scans<T>::append_scan(beam_count, beam_size, elevation, cell_size);
545  }
546 
555  PolarScan<T>& make_scan(unsigned idx, unsigned beam_size, double elevation, double cell_size)
556  {
557  return volume::Scans<T>::make_scan(idx, beam_count, beam_size, elevation, cell_size);
558  }
559 
564  Volume& operator*=(const T coefficient)
565  {
566  for(unsigned el=0;el<this->size();el++)
567  {
568  this->scan(el)*=coefficient;
569  }
570  return *this;
571  }
572 
578  {
579  for (unsigned el=0;el<this->size();el++)
580  {
581  this->scan(el)+=addend[el];
582  }
583  return *this;
584  }
585 
586 protected:
587  void resize_elev_fin();
588 
589 private:
590  Volume& operator=(Volume&);
591  Volume(const Volume&);
592 };
593 
594 /*
595  * Tell the compiler that the implementation of these templates is already
596  * explicitly instantiated in volume.cpp, and it should not spend time
597  * reinstantiating it every time the template is used.
598  */
599 extern template class PolarScan<double>;
600 extern template class Volume<double>;
601 namespace volume {
602 extern template class Scans<double>;
603 }
604 
605 }
606 
607 #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:311
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:498
Volume & operator*=(const T coefficient)
*= operator defined
Definition: volume.h:564
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:542
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:432
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:449
PolarScan< T > & append_scan(unsigned beam_count, unsigned beam_size, double elevation, double cell_size)
Append a scan to this volume.
Definition: volume.h:330
void normalize_elevations(const std::vector< double > &elevations)
Change the elevations in the PolarScans to match the given elevation vector.
Definition: volume.h:395
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:490
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:473
Volume(unsigned beam_count)
Constructor.
Definition: volume.h:438
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:464
Homogeneous volume with a common beam count for all PolarScans.
Definition: volume.h:428
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:286
const unsigned max_beam_size() const
Return the maximum beam size in all PolarScans.
Definition: volume.h:455
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:355
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:315
double elevation_max() const
Return the highest elevation.
Definition: volume.h:479
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:555
time_t acq_date
Acquisition date.
Definition: volume.h:247
Volume & operator+=(Volume &addend)
+= operator defined
Definition: volume.h:577
void SetOffset(T offset)
set offset value
Definition: volume.h:301
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.