Elaboradar  0.1
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))
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;
251  std::string source_name;
252 
253  LoadInfo()
254  : declutter_rsp(false)
255  {
256  }
257 };
258 
262 template<typename T>
263 class Scans : public std::vector<PolarScan<T>>
264 {
265 public:
266  typedef T Scalar;
268  std::string quantity;
270  std::string units;
272  std::shared_ptr<LoadInfo> load_info;
274  RadarSite radarSite;
276  T offset = 0;
277 
278  // FIXME: Cannot use '= default' or g++ 4.8.3 will not instantiate it on explicit template instantiation
279  Scans() : std::vector<PolarScan<T>>() {}
280 
287  template<typename OT>
288  Scans(const Scans<OT>& v, const T& default_value)
289  {
290  this->quantity = v.quantity;
291  this->units = v.units;
292  this->load_info = v.load_info;
293  this->reserve(v.size());
294  this->radarSite = v.radarSite;
295  for (const auto& src_scan : v)
296  this->push_back(PolarScan<T>(src_scan, default_value));
297  this->offset=v.offset;
298  }
299 
302  //
303  void SetOffset (T offset)
304  {
305  this->offset = offset;
306  }
307 
308 
309 
313  PolarScan<T>& scan(unsigned idx) { return (*this)[idx]; }
317  const PolarScan<T>& scan(unsigned idx) const { return (*this)[idx]; }
318 
332  PolarScan<T>& append_scan(unsigned beam_count, unsigned beam_size, double elevation, double cell_size, const T& default_value=algo::DBZ::BYTEtoDB(1))
333  {
334  // Ensure elevations grow as scan indices grow
335  if (!this->empty() && elevation <= this->back().elevation)
336  {
337  LOG_CATEGORY("radar.io");
338  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);
339  throw std::runtime_error("elevation not greather than the last one");
340  }
341  // Add the new polar scan
342  this->push_back(PolarScan<T>(beam_count, beam_size, default_value));
343  this->back().elevation = elevation;
344  this->back().cell_size = cell_size;
345  return this->back();
346  }
347 
357  PolarScan<T>& make_scan(unsigned idx, unsigned beam_count, unsigned beam_size, double elevation, double cell_size)
358  {
359  if (idx < this->size())
360  {
361  if (beam_count != (*this)[idx].beam_count)
362  {
363  LOG_CATEGORY("radar.io");
364  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);
365  throw std::runtime_error("beam_size mismatch");
366  }
367  if (beam_size != (*this)[idx].beam_size)
368  {
369  LOG_CATEGORY("radar.io");
370  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);
371  throw std::runtime_error("beam_size mismatch");
372  }
373  } else {
374  // If some elevation has been skipped, fill in the gap
375  if (idx > this->size())
376  {
377  if (this->empty())
378  this->push_back(PolarScan<T>(beam_count, beam_size));
379  while (this->size() < idx)
380  this->push_back(PolarScan<T>(beam_count, this->back().beam_size));
381  }
382 
383  // Add the new polar scan
384  this->push_back(PolarScan<T>(beam_count, beam_size));
385  this->back().elevation = elevation;
386  this->back().cell_size = cell_size;
387  }
388 
389  // Return it
390  return (*this)[idx];
391  }
392 
397  void normalize_elevations(const std::vector<double>& elevations)
398  {
399  // Ensure that we have enough standard elevations
400  if (elevations.size() < this->size())
401  {
402  LOG_CATEGORY("radar.io");
403  LOG_ERROR("normalize_elevations: standard elevation array has %zd elements, but we have %zd scans",
404  elevations.size(), this->size());
405  throw std::runtime_error("not enough standard elevations");
406  }
407  // Ensure that the nudging that we do do not confuse a scan
408  // with another
409  for (size_t i = 0; i < this->size() - 1; ++i)
410  {
411  if (abs(elevations[i] - this->at(i).elevation) > abs(elevations[i] - this->at(i + 1).elevation))
412  {
413  LOG_CATEGORY("radar.io");
414  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);
415  throw std::runtime_error("real elevation is too different than standard elevations");
416  }
417  }
418  // Assign the new elevations
419  for (size_t i = 0; i < this->size(); ++i)
420  this->at(i).elevation = elevations[i];
421  }
422 };
423 
424 }
425 
429 template<typename T>
430 class Volume : public volume::Scans<T>
431 {
432 public:
434  const unsigned beam_count;
435 
440  Volume(unsigned beam_count)
442  {
443  }
444 
450  template<typename OT>
451  Volume(const Volume<OT>& v, const T& default_value)
452  : volume::Scans<T>(v, default_value), beam_count(v.beam_count)
453  {
454  }
455 
457  const unsigned max_beam_size() const
458  {
459  unsigned res = 0;
460  for (const auto& scan : *this)
461  res = std::max(res, scan.beam_size);
462  return res;
463  }
464 
466  bool is_unique_cell_size() const
467  {
468  double cell_size = this->at(0).cell_size;
469  for (size_t i = 1; i < this->size(); ++i)
470  if ( this->at(i).cell_size != cell_size) return false;
471  return true;
472  }
473 
475  double elevation_min() const
476  {
477  return this->front().elevation;
478  }
479 
481  double elevation_max() const
482  {
483  return this->back().elevation;
484  }
485 
492  void read_vertical_slice(unsigned az, Matrix2D<T>& slice, double missing_value) const
493  {
494  unsigned size_z = std::max(this->size(), (size_t)slice.rows());
495  for (unsigned el = 0; el < size_z; ++el)
496  this->scan(el).read_beam(az, slice.row_ptr(el), slice.cols(), missing_value);
497  }
498 
500  void compute_stats(VolumeStats& stats) const
501  {
502  stats.count_zeros.resize(this->size());
503  stats.count_ones.resize(this->size());
504  stats.count_others.resize(this->size());
505  stats.sum_others.resize(this->size());
506 
507  for (unsigned iel = 0; iel < this->size(); ++iel)
508  {
509  stats.count_zeros[iel] = 0;
510  stats.count_ones[iel] = 0;
511  stats.count_others[iel] = 0;
512  stats.sum_others[iel] = 0;
513 
514  for (unsigned iaz = 0; iaz < this->scan(iel).beam_count; ++iaz)
515  {
516  for (size_t i = 0; i < this->scan(iel).beam_size; ++i)
517  {
518  int val = algo::DBZ::DBtoBYTE(this->scan(iel).get(iaz, i));
519  switch (val)
520  {
521  case 0: stats.count_zeros[iel]++; break;
522  case 1: stats.count_ones[iel]++; break;
523  default:
524  stats.count_others[iel]++;
525  stats.sum_others[iel] += val;
526  break;
527  }
528  }
529  }
530  }
531  }
532 
545  {
547  }
548 
557  PolarScan<T>& make_scan(unsigned idx, unsigned beam_size, double elevation, double cell_size)
558  {
560  }
561 
566  Volume& operator*=(const T coefficient)
567  {
568  for(unsigned el=0;el<this->size();el++)
569  {
570  this->scan(el)*=coefficient;
571  }
572  return *this;
573  }
574 
580  {
581  for (unsigned el=0;el<this->size();el++)
582  {
583  this->scan(el)+=addend[el];
584  }
585  return *this;
586  }
587 
588 protected:
589  void resize_elev_fin();
590 
591 private:
592  Volume& operator=(Volume&);
593  Volume(const Volume&);
594 };
595 
596 /*
597  * Tell the compiler that the implementation of these templates is already
598  * explicitly instantiated in volume.cpp, and it should not spend time
599  * reinstantiating it every time the template is used.
600  */
601 extern template class PolarScan<double>;
602 extern template class Volume<double>;
603 namespace volume {
604 extern template class Scans<double>;
605 }
606 
607 }
608 
609 #endif
Radar Site description.
PolarScan(const PolarScan &s)
Constructor Create a copy of a PolarScan.
Definition: volume.h:142
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 nodata
Value used as 'no data' value.
Definition: volume.h:116
double undetect
Minimum amount that can be measured.
Definition: volume.h:118
double offset
Conversion factor.
Definition: volume.h:122
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
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
T get(unsigned az, unsigned beam) const
Get a beam value.
Definition: volume.h:161
double gain
Conversion factor.
Definition: volume.h:120
void set(unsigned az, unsigned beam, T val)
Set a beam value.
Definition: volume.h:170
PolarScan - structure to describe a polarScan containing a matrix of data and conversion factors.
Definition: volume.h:113
void compute_stats(VolumeStats &stats) const
Compute Volume statistics.
Definition: volume.h:500
Volume(unsigned beam_count)
Constructor.
Definition: volume.h:440
double elevation_max() const
Return the highest elevation.
Definition: volume.h:481
PolarScan< T > & append_scan(unsigned beam_size, double elevation, double cell_size)
Append a scan to this volume.
Definition: volume.h:544
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:492
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:557
Volume & operator*=(const T coefficient)
*= operator defined
Definition: volume.h:566
const unsigned max_beam_size() const
Return the maximum beam size in all PolarScans.
Definition: volume.h:457
Volume(const Volume< OT > &v, const T &default_value)
Copy constructor.
Definition: volume.h:451
double elevation_min() const
Return the lowest elevation.
Definition: volume.h:475
Volume & operator+=(Volume &addend)
+= operator defined
Definition: volume.h:579
const unsigned beam_count
Number of beam_count used ast each elevations.
Definition: volume.h:434
bool is_unique_cell_size() const
Test if same cell_size in all PolarScans.
Definition: volume.h:466
Homogeneous volume with a common beam count for all PolarScans.
Definition: volume.h:431
static unsigned char DBtoBYTE(double DB, double gain=80./255., double offset=-20.)
funzione che converte dB in valore intero tra 0 e 255
Definition: dbz.h:94
static constexpr double BYTEtoDB(unsigned char DBZbyte, double gain=80./255., double offset=-20.)
funzione che converte Z unsigned char in DBZ
Definition: dbz.h:82
const PolarScan< T > & scan(unsigned idx) const
Access a polar scan (const)
Definition: volume.h:317
void SetOffset(T offset)
set offset value
Definition: volume.h:303
RadarSite radarSite
RadarSite.
Definition: volume.h:274
Scans(const Scans< OT > &v, const T &default_value)
Constructor Copy from another Scans.
Definition: volume.h:288
std::string units
Data units according to ODIM documentation.
Definition: volume.h:270
PolarScan< T > & append_scan(unsigned beam_count, unsigned beam_size, double elevation, double cell_size, const T &default_value=algo::DBZ::BYTEtoDB(1))
Append a scan to this volume.
Definition: volume.h:332
T offset
Data Offset.
Definition: volume.h:276
PolarScan< T > & scan(unsigned idx)
Access a polar scan.
Definition: volume.h:313
std::string quantity
Odim quantity name.
Definition: volume.h:268
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:357
void normalize_elevations(const std::vector< double > &elevations)
Change the elevations in the PolarScans to match the given elevation vector.
Definition: volume.h:397
std::shared_ptr< LoadInfo > load_info
Polar volume information.
Definition: volume.h:272
Sequence of PolarScans which can have a different beam count for each elevation.
Definition: volume.h:264
String functions.
Definition: cart.cpp:4
Base for all matrices we use, since we rely on row-major data.
Definition: matrix.h:37
double diff_height(unsigned rg_start, unsigned rg_end)
Height difference in kilometers (legacy) between two range gates.
Definition: volume.cpp:32
unsigned beam_size
Number of samples in each beam.
Definition: volume.h:33
double elevation
Nominal elevation of this PolarScan, which may be different from the effective elevation of each sing...
Definition: volume.h:42
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
Eigen::VectorXd elevations_real
Vector of actual elevations for each beam.
Definition: volume.h:45
unsigned beam_count
Count of beams in this scan.
Definition: volume.h:30
Eigen::VectorXd azimuths_real
Vector of actual azimuths for each beam.
Definition: volume.h:36
double cell_size
Size of a beam cell in meters.
Definition: volume.h:48
double sample_height(unsigned cell_idx) const
Return the height (in meters) of the sample at the given cell indexa.
Definition: volume.cpp:37
Basic structure to describe a polar scan, independently of the type of its samples.
Definition: volume.h:28
time_t acq_date
Acquisition date.
Definition: volume.h:247
std::string source_name
radar source name in capital letters (SPC or GAT)
Definition: volume.h:251
std::string filename
Original file name.
Definition: volume.h:243
bool declutter_rsp
flag true if data have been decluttered with Doppler at rsp level
Definition: volume.h:249
LoadInfo structure - Contains generic volume information.
Definition: volume.h:241