Elaboradar  0.1
odim.cpp
1 #include "odim.h"
2 #include "logging.h"
3 #include "utils.h"
4 #include <radarlib/radar.hpp>
5 #include <memory>
6 
7 using namespace std;
8 using namespace OdimH5v21;
9 
10 namespace {
11 
12 unsigned int_to_unsigned(int val, const char* desc)
13 {
14  if (val < 0)
15  {
16  string errmsg(desc);
17  errmsg += " has a negative value";
18  throw std::runtime_error(errmsg);
19  }
20  return (unsigned)val;
21 }
22 
23 }
24 
25 
26 namespace radarelab {
27 namespace volume {
28 
29 void ODIMLoader::request_quantity(const std::string& name, Scans<double>* volume)
30 {
31  to_load.insert(make_pair(name, volume));
32 }
33 
34 void ODIMLoader::load(const std::string& pathname)
35 {
36  LOG_CATEGORY("radar.io");
37 
38  namespace odim = OdimH5v21;
39  using namespace Radar;
40  using namespace std;
41 
42  shared_ptr<LoadInfo> load_info = make_shared<LoadInfo>();
43  load_info->filename = pathname;
44 
45  unique_ptr<odim::OdimFactory> factory(new odim::OdimFactory());
46  unique_ptr<odim::PolarVolume> volume(factory->openPolarVolume(pathname));
47 
48  load_info->acq_date = volume->getDateTime();
49 
50  try
51  {
52  //string sourcename(volume->getSource().Place.substr(2).c_str());
53  //transform(sourcename.begin(), sourcename.end(), sourcename.begin(),::toupper);
54  //cout<<"volume getsource : "<<sourcename<<endl;
55 // string sourcename(volume->getSource().Place.c_str());
56 // cout<<"volume getsource : "<<sourcename<<endl;
57  load_info->source_name = volume->getSource().Place.c_str() ; //sourcename;
58  }
59  catch (std::exception& stde)
60  {
61  std::cerr << "Errore durante l'esecuzione: " << stde.what() << std::endl;
62  }
63  catch (...)
64  {
65  std::cerr << "Errore sconosciuto" << std::endl;
66  }
67 
68  double range_scale = 0;
69 
70  // Iterate all scans
71  unsigned scan_count = int_to_unsigned(volume->getScanCount(), "scan count");
72  double old_elevation = -1000.;
73  for (unsigned src_elev = 0; src_elev < scan_count; ++src_elev)
74  {
75  unique_ptr<odim::PolarScan> scan(volume->getScan(src_elev));
76  double elevation = scan->getEAngle();
77  Available_Elevations.push_back(elevation);
78  // FIXME: please add a comment on why this is needed: are there faulty
79  // ODIM files out there which repeat elevations? Is it correct that
80  // only the first elevation is used and not, say, the last one?
81  if( elevation == old_elevation ) continue;
82  old_elevation=elevation;
83  // Read and and validate resolution information
84  if (src_elev == 0)
85  range_scale = scan->getRangeScale();
86  else {
87  double rs = scan->getRangeScale();
88  if (rs != range_scale)
89  {
90  LOG_ERROR("scan %d (elevation %f) has rangeScale %f that is different from %f in the previous scans",
91  src_elev, elevation, rs, range_scale);
92  throw runtime_error("rangeScale mismatch");
93  }
94  }
95 
96  // Get and validate the azimuth angles for this scan
97  std::vector<odim::AZAngles> azangles = scan->getAzimuthAngles();
98  int rpm_sign = scan->getDirection();
99 
100  std::vector<double> elevation_angles = scan->getElevationAngles();
101 
102  unsigned beam_count = int_to_unsigned(scan->getNumRays(), "number of rays");
103  if (azangles.size() != beam_count)
104  {
105  LOG_ERROR("elevation %f has %zd azimuth angles and %d rays", elevation, azangles.size(), beam_count);
106  throw runtime_error("mismatch between number of azumuth angles and number of rays");
107  }
108 
109  unsigned beam_size = int_to_unsigned(scan->getNumBins(), "beam size");
110 
111  // Read all quantities that have been requested
112  for (auto& todo : to_load)
113  {
114  // Create azimuth maps and PolarScan objects for this elevation
115  const string& name = todo.first;
116  Scans<double>& target = *todo.second;
117 
118  // Pick the best quantity among the ones available
119  if (!scan->hasQuantityData(name))
120  {
121  LOG_WARN("no %s found for elevation angle %f: skipping", name.c_str(), elevation);
122  continue;
123  }
124  PolarScan<double>& vol_pol_scan = target.append_scan(beam_count, beam_size, elevation, range_scale);
125 
126  unique_ptr<odim::PolarScanData> data(scan->getQuantityData(name));
127 
128  // Fill/overwrite variable metadata at volume level
129  target.quantity = name;
130  target.radarSite.height_r = volume->getAltitude()/1000.;
131  target.radarSite.antennaTowerHeight = 0.;
132  // Fill variable metadata at scan level
133  vol_pol_scan.nodata = data->getNodata() * data->getGain() + data->getOffset();
134  vol_pol_scan.undetect = data->getUndetect() * data->getGain() + data->getOffset();
135  vol_pol_scan.gain = data->getGain();
136  vol_pol_scan.offset = data->getOffset();
137  vol_pol_scan.cell_size = scan->getRangeScale();
138 
139  // Read actual data from ODIM
140  odim::RayMatrix<double> matrix;
141  matrix.resize(beam_count, beam_size);
142  data->readTranslatedData(matrix);
143 
144  for (unsigned src_az = 0; src_az < beam_count; ++src_az)
145  {
146  Eigen::VectorXd beam(beam_size);
147 
148  for (unsigned i = 0; i < beam_size; ++i)
149  beam(i) = matrix.elem(src_az, i);
150 
151  vol_pol_scan.row(src_az) = beam;
152  vol_pol_scan.elevations_real(src_az) = elevation_angles[src_az];
153  vol_pol_scan.azimuths_real(src_az) = azangles[src_az].averagedAngle(rpm_sign);
154  }
155  }
156  }
157  for (auto& i: to_load)
158  {
159  i.second->load_info = load_info;
160  }
161 }
162 
163 void ODIMStorer::store(const std::string& pathname)
164 {
165  namespace odim = OdimH5v21;
166  using namespace Radar;
167  using namespace std;
168 
169  shared_ptr<LoadInfo> load_info = make_shared<LoadInfo>();
170  load_info->filename = pathname;
171 
172  unique_ptr<odim::OdimFactory> factory(new odim::OdimFactory());
173  unique_ptr<odim::PolarVolume> volume(factory->openPolarVolume(pathname));
174 
175  for(unsigned i=0;i<to_store_int.size();i++)
176  for(unsigned j=0;j<to_store_int[i]->size();j++)
177  {
178  vector<odim::PolarScan*> scans;
179  scans = volume->getScans(to_store_int[i]->scan(j).elevation,0.1);
180  shared_ptr<odim::PolarScan> scan;
181  if(scans.size()) scan.reset(scans[0]);
182  else
183  {
184  scan.reset(volume->createScan());
185  scan->setEAngle(to_store_int[i]->scan(j).elevation);
186  }
187  unique_ptr<odim::PolarScanData> data(scan->createQuantityData(to_store_int[i]->quantity));
188  data->setGain((double)to_store_int[i]->scan(j).gain);
189  data->setOffset((double)to_store_int[i]->scan(j).offset);
190  data->setNodata((double)to_store_int[i]->scan(j).nodata);
191  data->setUndetect((double)to_store_int[i]->scan(j).undetect);
192  odim::RayMatrix<unsigned short> matrix;
193  matrix.resize(to_store_int[i]->scan(j).beam_count,to_store_int[i]->scan(j).beam_size);
194  for(unsigned ii=0;ii<to_store_int[i]->scan(j).beam_count;ii++)
195  for(unsigned jj=0;jj<to_store_int[i]->scan(j).beam_size;jj++)
196  matrix.elem(ii,jj) = to_store_int[i]->scan(j)(ii,jj);
197  data->writeData(matrix);
198  }
199  for(unsigned i=0;i<to_store_fp.size();i++)
200  for(unsigned j=0;j<to_store_fp[i]->size();j++)
201  {
202  vector<odim::PolarScan*> scans;
203  scans = volume->getScans(to_store_fp[i]->scan(j).elevation,0.1);
204  shared_ptr<odim::PolarScan> scan;
205  if(scans.size()) scan.reset(scans[0]);
206  else
207  {
208  scan.reset(volume->createScan());
209  scan->setEAngle(to_store_fp[i]->scan(j).elevation);
210  }
211  unique_ptr<odim::PolarScanData> data(scan->createQuantityData(to_store_fp[i]->quantity));
212  data->setGain((double)to_store_fp[i]->scan(j).gain);
213  data->setOffset((double)to_store_fp[i]->scan(j).offset);
214  data->setNodata((double)to_store_fp[i]->scan(j).nodata);
215  data->setUndetect((double)to_store_fp[i]->scan(j).undetect);
216  odim::RayMatrix<float> matrix;
217  matrix.resize(to_store_fp[i]->scan(j).beam_count,to_store_fp[i]->scan(j).beam_size);
218  for(unsigned ii=0;ii<to_store_fp[i]->scan(j).beam_count;ii++)
219  for(unsigned jj=0;jj<to_store_fp[i]->scan(j).beam_size;jj++)
220  matrix.elem(ii,jj) = to_store_fp[i]->scan(j)(ii,jj);
221  data->writeAndTranslate(matrix,(float)data->getOffset(),(float)data->getGain(),H5::PredType::NATIVE_UINT16);
222  }
223 
224  for(unsigned i=0;i<to_store_uchar.size();i++)
225  for(unsigned j=0;j<to_store_uchar[i]->size();j++)
226  {
227  vector<odim::PolarScan*> scans;
228  scans = volume->getScans(to_store_uchar[i]->scan(j).elevation,0.1);
229  shared_ptr<odim::PolarScan> scan;
230  if(scans.size()) scan.reset(scans[0]);
231  else
232  {
233  scan.reset(volume->createScan());
234  scan->setEAngle(to_store_uchar[i]->scan(j).elevation);
235  }
236  unique_ptr<odim::PolarScanData> data(scan->createQuantityData(to_store_uchar[i]->quantity));
237  data->setGain((double)to_store_uchar[i]->scan(j).gain);
238  data->setOffset((double)to_store_uchar[i]->scan(j).offset);
239  data->setNodata((double)to_store_uchar[i]->scan(j).nodata);
240  data->setUndetect((double)to_store_uchar[i]->scan(j).undetect);
241  odim::RayMatrix<float> matrix;
242  matrix.resize(to_store_uchar[i]->scan(j).beam_count,to_store_uchar[i]->scan(j).beam_size);
243  for(unsigned ii=0;ii<to_store_uchar[i]->scan(j).beam_count;ii++)
244  for(unsigned jj=0;jj<to_store_uchar[i]->scan(j).beam_size;jj++)
245  matrix.elem(ii,jj) = to_store_uchar[i]->scan(j)(ii,jj);
246  data->writeAndTranslate(matrix,(float)data->getOffset(),(float)data->getGain(),H5::PredType::NATIVE_UINT8);
247  }
248 
249  for(unsigned i=0;i<to_replace.size();i++)
250  for(unsigned j=0;j<to_replace[i]->size();j++)
251  {
252  vector<odim::PolarScan*> scans;
253  scans = volume->getScans(to_replace[i]->scan(j).elevation,0.1);
254 
255  shared_ptr<odim::PolarScan> scan;
256  if(scans.size()) scan.reset(scans[0]);
257  else
258  {
259 // scan.reset(volume->createScan());
260 // scan->setEAngle(to_store_uchar[i]->scan(j).elevation);
261  ; // da sistemare. in questo caso bisogna far uscire una segnalazione perchè qui non dovremmo proprio finirci.
262  }
263 
264  unique_ptr<odim::PolarScanData> data(scan->getQuantityData(to_replace[i]->quantity));
265  H5::AtomType type = data->getDataType();
266 
267  odim::RayMatrix<double> matrix;
268  matrix.resize(to_replace[i]->scan(j).beam_count,to_replace[i]->scan(j).beam_size);
269  for(unsigned ii=0;ii<to_replace[i]->scan(j).beam_count;ii++)
270  for(unsigned jj=0;jj<to_replace[i]->scan(j).beam_size;jj++)
271  matrix.elem(ii,jj) = to_replace[i]->scan(j)(ii,jj);
272  data->writeAndTranslate(matrix,(float)data->getOffset(),(float)data->getGain(),type);
273  }
274 
275 }
276 
277 
278 void ODIMStorer::storeQuality(const std::string& pathname, const std::string & task, bool RemoveQualityFields)
279 {
280  namespace odim = OdimH5v21;
281  using namespace Radar;
282  using namespace std;
283 
284  shared_ptr<LoadInfo> load_info = make_shared<LoadInfo>();
285  load_info->filename = pathname;
286 
287  unique_ptr<odim::OdimFactory> factory(new odim::OdimFactory());
288  unique_ptr<odim::PolarVolume> volume(factory->openPolarVolume(pathname));
289 
290  for(unsigned i=0;i<to_store_uchar.size();i++)
291  for(unsigned j=0;j<to_store_uchar[i]->size();j++)
292  {
293  vector<odim::PolarScan*> scans;
294  scans = volume->getScans(to_store_uchar[i]->scan(j).elevation,0.1);
295  shared_ptr<odim::PolarScan> scan;
296  if(scans.size()) scan.reset(scans[0]);
297  else
298  {
299  scan.reset(volume->createScan());
300  scan->setEAngle(to_store_uchar[i]->scan(j).elevation);
301  }
302  if (RemoveQualityFields) {
303  int iqc =scan->getQualityCount();
304  for (int iq=iqc-1; iq >=0 ; iq--){
305  scan->removeQuality(iq);
306  }
307  }
308  unique_ptr<odim::OdimQuality> quality(scan->createQuality());
309  quality->getWhat()->set(ATTRIBUTE_WHAT_OFFSET, 0.);
310  quality->getWhat()->set(ATTRIBUTE_WHAT_GAIN, 1.);
311  quality->getHow() ->set(ATTRIBUTE_HOW_TASK, task);
312  odim::RayMatrix<unsigned short> matrix;
313  matrix.resize(to_store_uchar[i]->scan(j).beam_count,to_store_uchar[i]->scan(j).beam_size);
314  for(unsigned ii=0;ii<to_store_uchar[i]->scan(j).beam_count;ii++)
315  for(unsigned jj=0;jj<to_store_uchar[i]->scan(j).beam_size;jj++)
316  matrix.elem(ii,jj) = to_store_uchar[i]->scan(j)(ii,jj);
317  quality->writeQuality(matrix);
318  }
319 
320 }
321 
322 
323 } // namespace volume
324 } // namespace radarelab
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
double gain
Conversion factor.
Definition: volume.h:120
PolarScan - structure to describe a polarScan containing a matrix of data and conversion factors.
Definition: volume.h:113
RadarSite radarSite
RadarSite.
Definition: volume.h:274
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
std::string quantity
Odim quantity name.
Definition: volume.h:268
String functions.
Definition: cart.cpp:4
Codice per il caricamento di volumi ODIM in radarelab.
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
double cell_size
Size of a beam cell in meters.
Definition: volume.h:48