Elaboradar  0.1
sp20.cpp
1 #include "sp20.h"
2 #include "logging.h"
3 #include "utils.h"
4 #include <radarlib/radar.hpp>
5 #include <memory>
6 #include <cerrno>
7 #include <cstring>
8 #include "sp20/func_SP20read.h"
9 
10 #ifdef NEL
11 #undef NEL
12 #endif
13 
14 using namespace std;
15 using namespace radarelab::algo;
16 
18 int elev_array[15];
19 
20 namespace radarelab {
21 namespace volume {
22 
23 namespace sp20 {
24 
25 struct Beam
26 {
27  BEAM_HD_SP20_INFO beam_info;
28  unsigned char data_z[1024];
29  unsigned char data_d[1024];
30  unsigned char data_v[1024];
31  unsigned char data_w[1024];
32  unsigned beam_size = 0;
33 
34  // Read the next beam in the file. Returns false on end of file.
35  bool read(File& in)
36  {
37  unsigned char beam_header[40];
38  if (!in.fread(beam_header, sizeof(beam_header)))
39  return false;
40 
41  decode_header_sp20_date_from_name(beam_header, &beam_info, (char*)in.name());
42 
43  // Salta beam non validi
44  if (!beam_info.valid_data){
45  unsigned to_be_skipped=0;
46  if (has_z()) to_be_skipped=to_be_skipped+beam_info.cell_num;
47  if (has_d()) to_be_skipped=to_be_skipped+beam_info.cell_num;
48  if (has_v()) to_be_skipped=to_be_skipped+beam_info.cell_num;
49  if (has_w()) to_be_skipped=to_be_skipped+beam_info.cell_num;
50  in.fseek(to_be_skipped, SEEK_CUR);
51  return true;
52  }
53 
54  beam_size = beam_info.cell_num;
55 
56  if (has_z()) in.fread(data_z, beam_info.cell_num);
57  if (has_d()) in.fread(data_d, beam_info.cell_num);
58  if (has_v()) in.fread(data_v, beam_info.cell_num);
59  if (has_w()) in.fread(data_w, beam_info.cell_num);
60 
61  return true;
62  }
63 
64  bool has_z() const { return beam_info.flag_quantities[0]; }
65  bool has_d() const { return beam_info.flag_quantities[1]; }
66  bool has_v() const { return beam_info.flag_quantities[2]; }
67  bool has_w() const { return beam_info.flag_quantities[3]; }
68 };
69 
70 }
71 
72 namespace {
73 
74 struct Beams : public std::vector<sp20::Beam*>
75 {
76  // Elevation
77  double elevation;
78  // Beam size for the polar scan
79  unsigned beam_size = 0;
80 
81  Beams(double elevation) : elevation(elevation) {}
82  Beams(const Beams&) = delete;
83  Beams(const Beams&&) = delete;
84  ~Beams()
85  {
86  for (auto i: *this)
87  delete i;
88  }
89  Beams& operator=(const Beams&) = delete;
90 
91  // Compute the minimum beam size
92  unsigned min_beam_size() const
93  {
94  if (empty()) return 0;
95  unsigned res = (*this)[0]->beam_size;
96  for (auto i: *this)
97  res = min(res, i->beam_size);
98  return res;
99  }
100 };
101 
102 struct Elevations : public std::vector<unique_ptr<Beams>>
103 {
104  Elevations(float* elevations, unsigned el_count)
105  {
106  // Create a sorted vector of elevations
107  vector<double> els;
108  for (unsigned i = 0; i < el_count; ++i)
109  els.push_back(elevations[i]);
110  std::sort(els.begin(), els.end());
111 
112  // Use it to initialize our Beams structures
113  for (double el: els)
114  emplace_back(new Beams(el));
115  }
116 
117  Beams* get_closest(double elevation)
118  {
119  for (unsigned i=0; i < size(); ++i)
120  if (elevation >= (at(i)->elevation-0.5) && elevation < (at(i)->elevation+0.5))
121  return at(i).get();
122  return nullptr;
123  }
124 };
125 
126 }
127 
128 void SP20Loader::load(const std::string& pathname)
129 {
130  namespace odim = OdimH5v21;
131  LOG_CATEGORY("Volume");
132  // dimensioni cella a seconda del tipo di acquisizione
133  static const float size_cell_by_resolution[]={62.5,125.,250.,500.,1000.,2000.};
134  //LOG_CATEGORY("radar.io");
135  HD_DBP_SP20_RAW hd_char;
136 
137  shared_ptr<LoadInfo> load_info = make_shared<LoadInfo>();
138  load_info->filename = pathname;
139 
140  // Replicato qui la read_dbp_SP20, per poi metterci mano e condividere codice con la lettura di ODIM
141 
142  File sp20_in(logging_category);
143  sp20_in.open(pathname, "rb", "input sp20 file");
144  if (!sp20_in)
145  throw std::runtime_error("failed to open sp20 input file");
146 
147  // Read and decode file header
148  sp20_in.fread(&hd_char, sizeof(hd_char));
149 
150  HD_DBP_SP20_DECOD hd_file;
151  decode_header_DBP_SP20(&hd_char, &hd_file);
152 
153  /*--------
154  ATTENZIONE PRENDO LA DATA DAL NOME DEL FILE
155  -------*/
156  struct tm data_nome;
157  time_t acq_date = get_date_from_name(0, &data_nome, pathname.c_str());
158  load_info->acq_date = acq_date;
159  load_info->declutter_rsp = (bool)hd_file.filtro_clutter;
160  double size_cell = size_cell_by_resolution[(int)hd_file.cell_size];
161  bool has_dual_prf = hd_file.Dual_PRF;
162 
163  // Read all beams from the file
164  Elevations elevations(hd_file.ele, hd_file.num_ele);
165 
166  while (true)
167  {
168  unique_ptr<sp20::Beam> beam(new sp20::Beam);
169 
170  if (!beam->read(sp20_in))
171  break;
172 
173  Beams* beams = elevations.get_closest(beam->beam_info.elevation);
174  if (!beams) continue;
175 
176  beams->push_back(beam.release());
177  }
178 
179  if (elevations.empty())
180  throw std::runtime_error("sp20 file has no beams");
181 
182  // Set the beam size on each elevation to the minimum beam size
183  for (auto& beams: elevations)
184  {
185  unsigned beam_size = beams->min_beam_size();
186  beams->beam_size = beam_size;
187  }
188 
189  if (elevations.back()->beam_size == 0)
190  throw std::runtime_error("last elevation beam size is 0");
191 
192  // For elevations that have beam_size = 0, use the beam size from the
193  // elevation above
194  for (int i = elevations.size() - 2; i >= 0; --i)
195  {
196  if (elevations[i]->beam_size == 0)
197  elevations[i]->beam_size = elevations[i + 1]->beam_size;
198  }
199 
200  // Create the polarscans and fill them with data
201  for (unsigned idx = 0; idx < elevations.size(); ++idx)
202  {
203  const Beams& beams = *elevations[idx];
204 
205  // Create structures for a new PolarScan
206  if (vol_z) vol_z->append_scan(beams.size(), beams.beam_size, beams.elevation, size_cell);
207  if (vol_d) vol_d->append_scan(beams.size(), beams.beam_size, beams.elevation, size_cell);
208  if (vol_v) vol_v->append_scan(beams.size(), beams.beam_size, beams.elevation, size_cell);
209  if (vol_w) vol_w->append_scan(beams.size(), beams.beam_size, beams.elevation, size_cell);
210 
211  // Fill the new scan with data
212  for (unsigned i = 0; i < beams.size(); ++i)
213  beam_to_volumes(*beams.at(i), i, beams.beam_size, idx);
214  }
215 
216  if (vol_z)
217  {
218  vol_z->load_info = load_info;
219  if (load_info->declutter_rsp)
220  vol_z->quantity = odim::PRODUCT_QUANTITY_DBZH;
221  else
222  vol_z->quantity = odim::PRODUCT_QUANTITY_TH;
223  for (auto& scan : *vol_z)
224  {
225  scan.nodata = DBZ::DBtoBYTE(0);
226  scan.undetect = DBZ::DBtoBYTE(1);
227  scan.gain = 80.0 / 255.0;
228  scan.offset = -20;
229  }
230  }
231  if (vol_d)
232  {
233  vol_d->load_info = load_info;
234  vol_v->quantity = odim::PRODUCT_QUANTITY_ZDR;
235  for (auto& scan : *vol_z)
236  {
237  scan.nodata = -7;
238  scan.undetect = -6;
239  scan.gain = 16.0 / 255.0;
240  scan.offset = -6;
241  }
242  }
243  if (vol_v)
244  {
245  vol_v->load_info = load_info;
246  vol_v->quantity = odim::PRODUCT_QUANTITY_VRAD;
247  for (auto& scan : *vol_v)
248  {
249  if (has_dual_prf)
250  {
251  scan.nodata = -49.5;
252  scan.undetect = -49.5;
253  scan.gain = 99.0 / 254; // TODO: check why it's 254 and not 255
254  scan.offset = -49.5;
255  } else {
256  scan.nodata = -16.5;
257  scan.undetect = -16.5;
258  scan.gain = 33.0 / 254; // TODO: check why it's 254 and not 255
259  scan.offset = -16.5;
260  }
261  }
262  }
263  if (vol_w)
264  {
265  vol_w->load_info = load_info;
266  vol_w->quantity = odim::PRODUCT_QUANTITY_WRAD;
267  for (auto& scan : *vol_w)
268  {
269  scan.nodata = -1;
270  scan.undetect = 0;
271  scan.gain = 10.0 / 255.0;
272  scan.offset = 0;
273  }
274  }
275 
276  LOG_DEBUG ("Nel volume ci sono %zd scan", vol_z->size());
277  for (size_t i = 0; i < vol_z->size(); ++i)
278  LOG_DEBUG (" Scan %2zd - dimensione beam %5d - numero raggi %3d", i, vol_z->at(i).beam_size, vol_z->at(i).beam_count);
279 }
280 
281 void SP20Loader::beam_to_volumes(const sp20::Beam& beam, unsigned az_idx, unsigned beam_size, unsigned el_num)
282 {
283  const unsigned max_range = min(beam_size, beam.beam_size);
284 
285  if (vol_z && beam.has_z()) // Riflettività Z
286  {
287  // Convert to DB
288  Eigen::VectorXd dbs(max_range);
289  for (unsigned i = 0; i < max_range; ++i)
290  dbs(i) = DBZ::BYTEtoDB(beam.data_z[i]);
291  PolarScan<double>& scan = vol_z->at(el_num);
292  scan.row(az_idx) = dbs;
293  scan.elevations_real(az_idx) = beam.beam_info.elevation;
294  scan.azimuths_real(az_idx) = beam.beam_info.azimuth;
295  }
296 
297  if (vol_d && beam.has_d()) // Riflettività differenziale ZDR
298  {
299  // range variabilita ZDR
300  const double range_zdr = 16.;
301  const double min_zdr = -6.;
302 
303  // Convert to DB
304  Eigen::VectorXd dbs(max_range);
305  for (unsigned i = 0; i < max_range; ++i)
306  dbs(i) = beam.data_d[i] * range_zdr / 255. + min_zdr;
307 
308  PolarScan<double>& scan = vol_d->at(el_num);
309  scan.row(az_idx) = dbs;
310  scan.elevations_real(az_idx) = beam.beam_info.elevation;
311  scan.azimuths_real(az_idx) = beam.beam_info.azimuth;
312  }
313 
314  if (vol_v && beam.has_v()) // Velocità V
315  {
316  // Convert to m/s
317  Eigen::VectorXd ms(max_range);
318  if (beam.beam_info.PRF == 'S')
319  {
320  // range variabilita V - velocità radiale
321  const double range_v = 33.;
322  for (unsigned i = 0; i < max_range; ++i)
323  if (beam.data_v[i] == -128)
324  ms(i) = -range_v / 2;
325  else
326  ms(i) = beam.data_v[i] * range_v / 254.;
327  } else {
328  // range variabilita V - velocità radiale
329  const double range_v = 99.;
330  for (unsigned i = 0; i < max_range; ++i)
331  if (beam.data_v[i] == -128)
332  ms(i) = -range_v / 2;
333  else
334  ms(i) = beam.data_v[i] * range_v / 254.;
335  }
336  PolarScan<double>& scan = vol_v->at(el_num);
337  scan.row(az_idx) = ms;
338  scan.elevations_real(az_idx) = beam.beam_info.elevation;
339  scan.azimuths_real(az_idx) = beam.beam_info.azimuth;
340  }
341 
342  if (vol_w && beam.has_w()) // Spread - Sigma V
343  {
344  // range variabilita Sigma V - Spread velocità
345  const double range_sig_v = 10.;
346  // Convert to m/s
347  Eigen::VectorXd ms(max_range);
348  for (unsigned i = 0; i < max_range; ++i)
349  ms[i] = beam.data_w[i] * range_sig_v / 255.0;
350 
351  PolarScan<double>& scan = vol_w->at(el_num);
352  scan.row(az_idx) = ms;
353  scan.elevations_real(az_idx) = beam.beam_info.elevation;
354  scan.azimuths_real(az_idx) = beam.beam_info.azimuth;
355  }
356 }
357 
358 }
359 }
bool fread(void *buf, size_t size)
Performs a fread on the file, throwing an exception if anything goes wrong.
Definition: utils.cpp:109
Open a file taking its name from a given env variable.
Definition: utils.h:22
PolarScan - structure to describe a polarScan containing a matrix of data and conversion factors.
Definition: volume.h:113
String functions.
Definition: cart.cpp:4
Eigen::VectorXd elevations_real
Vector of actual elevations for each beam.
Definition: volume.h:45
Eigen::VectorXd azimuths_real
Vector of actual azimuths for each beam.
Definition: volume.h:36