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  double range_scale = 0;
51 
52  // Iterate all scans
53  unsigned scan_count = int_to_unsigned(volume->getScanCount(), "scan count");
54  double old_elevation = -1000.;
55  for (unsigned src_elev = 0; src_elev < scan_count; ++src_elev)
56  {
57  unique_ptr<odim::PolarScan> scan(volume->getScan(src_elev));
58  double elevation = scan->getEAngle();
59  Available_Elevations.push_back(elevation);
60  // FIXME: please add a comment on why this is needed: are there faulty
61  // ODIM files out there which repeat elevations? Is it correct that
62  // only the first elevation is used and not, say, the last one?
63  if( elevation == old_elevation ) continue;
64  old_elevation=elevation;
65  // Read and and validate resolution information
66  if (src_elev == 0)
67  range_scale = scan->getRangeScale();
68  else {
69  double rs = scan->getRangeScale();
70  if (rs != range_scale)
71  {
72  LOG_ERROR("scan %d (elevation %f) has rangeScale %f that is different from %f in the previous scans",
73  src_elev, elevation, rs, range_scale);
74  throw runtime_error("rangeScale mismatch");
75  }
76  }
77 
78  // Get and validate the azimuth angles for this scan
79  std::vector<odim::AZAngles> azangles = scan->getAzimuthAngles();
80  int rpm_sign = scan->getDirection();
81 
82  std::vector<double> elevation_angles = scan->getElevationAngles();
83 
84  unsigned beam_count = int_to_unsigned(scan->getNumRays(), "number of rays");
85  if (azangles.size() != beam_count)
86  {
87  LOG_ERROR("elevation %f has %zd azimuth angles and %d rays", elevation, azangles.size(), beam_count);
88  throw runtime_error("mismatch between number of azumuth angles and number of rays");
89  }
90 
91  unsigned beam_size = int_to_unsigned(scan->getNumBins(), "beam size");
92 
93  // Read all quantities that have been requested
94  for (auto& todo : to_load)
95  {
96  // Create azimuth maps and PolarScan objects for this elevation
97  const string& name = todo.first;
98  Scans<double>& target = *todo.second;
99 
100  // Pick the best quantity among the ones available
101  if (!scan->hasQuantityData(name))
102  {
103  LOG_WARN("no %s found for elevation angle %f: skipping", name.c_str(), elevation);
104  continue;
105  }
106  PolarScan<double>& vol_pol_scan = target.append_scan(beam_count, beam_size, elevation, range_scale);
107 
108  unique_ptr<odim::PolarScanData> data(scan->getQuantityData(name));
109 
110  // Fill/overwrite variable metadata at volume level
111  target.quantity = name;
112  target.radarSite.height_r = volume->getAltitude()/1000.;
113  target.radarSite.antennaTowerHeight = 0.;
114  // Fill variable metadata at scan level
115  vol_pol_scan.nodata = data->getNodata() * data->getGain() + data->getOffset();
116  vol_pol_scan.undetect = data->getUndetect() * data->getGain() + data->getOffset();
117  vol_pol_scan.gain = data->getGain();
118  vol_pol_scan.offset = data->getOffset();
119  vol_pol_scan.cell_size = scan->getRangeScale();
120 
121  // Read actual data from ODIM
122  odim::RayMatrix<double> matrix;
123  matrix.resize(beam_count, beam_size);
124  data->readTranslatedData(matrix);
125 
126  for (unsigned src_az = 0; src_az < beam_count; ++src_az)
127  {
128  Eigen::VectorXd beam(beam_size);
129 
130  for (unsigned i = 0; i < beam_size; ++i)
131  beam(i) = matrix.elem(src_az, i);
132 
133  vol_pol_scan.row(src_az) = beam;
134  vol_pol_scan.elevations_real(src_az) = elevation_angles[src_az];
135  vol_pol_scan.azimuths_real(src_az) = azangles[src_az].averagedAngle(rpm_sign);
136  }
137  }
138  }
139  for (auto& i: to_load)
140  {
141  i.second->load_info = load_info;
142  }
143 }
144 
145 void ODIMStorer::store(const std::string& pathname)
146 {
147  namespace odim = OdimH5v21;
148  using namespace Radar;
149  using namespace std;
150 
151  shared_ptr<LoadInfo> load_info = make_shared<LoadInfo>();
152  load_info->filename = pathname;
153 
154  unique_ptr<odim::OdimFactory> factory(new odim::OdimFactory());
155  unique_ptr<odim::PolarVolume> volume(factory->openPolarVolume(pathname));
156 
157  for(unsigned i=0;i<to_store_int.size();i++)
158  for(unsigned j=0;j<to_store_int[i]->size();j++)
159  {
160  vector<odim::PolarScan*> scans;
161  scans = volume->getScans(to_store_int[i]->scan(j).elevation,0.1);
162  shared_ptr<odim::PolarScan> scan;
163  if(scans.size()) scan.reset(scans[0]);
164  else
165  {
166  scan.reset(volume->createScan());
167  scan->setEAngle(to_store_int[i]->scan(j).elevation);
168  }
169  unique_ptr<odim::PolarScanData> data(scan->createQuantityData(to_store_int[i]->quantity));
170  data->setGain((double)to_store_int[i]->scan(j).gain);
171  data->setOffset((double)to_store_int[i]->scan(j).offset);
172  data->setNodata((double)to_store_int[i]->scan(j).nodata);
173  data->setUndetect((double)to_store_int[i]->scan(j).undetect);
174  odim::RayMatrix<unsigned short> matrix;
175  matrix.resize(to_store_int[i]->scan(j).beam_count,to_store_int[i]->scan(j).beam_size);
176  for(unsigned ii=0;ii<to_store_int[i]->scan(j).beam_count;ii++)
177  for(unsigned jj=0;jj<to_store_int[i]->scan(j).beam_size;jj++)
178  matrix.elem(ii,jj) = to_store_int[i]->scan(j)(ii,jj);
179  data->writeData(matrix);
180  }
181  for(unsigned i=0;i<to_store_fp.size();i++)
182  for(unsigned j=0;j<to_store_fp[i]->size();j++)
183  {
184  vector<odim::PolarScan*> scans;
185  scans = volume->getScans(to_store_fp[i]->scan(j).elevation,0.1);
186  shared_ptr<odim::PolarScan> scan;
187  if(scans.size()) scan.reset(scans[0]);
188  else
189  {
190  scan.reset(volume->createScan());
191  scan->setEAngle(to_store_fp[i]->scan(j).elevation);
192  }
193  unique_ptr<odim::PolarScanData> data(scan->createQuantityData(to_store_fp[i]->quantity));
194  data->setGain((double)to_store_fp[i]->scan(j).gain);
195  data->setOffset((double)to_store_fp[i]->scan(j).offset);
196  data->setNodata((double)to_store_fp[i]->scan(j).nodata);
197  data->setUndetect((double)to_store_fp[i]->scan(j).undetect);
198  odim::RayMatrix<float> matrix;
199  matrix.resize(to_store_fp[i]->scan(j).beam_count,to_store_fp[i]->scan(j).beam_size);
200  for(unsigned ii=0;ii<to_store_fp[i]->scan(j).beam_count;ii++)
201  for(unsigned jj=0;jj<to_store_fp[i]->scan(j).beam_size;jj++)
202  matrix.elem(ii,jj) = to_store_fp[i]->scan(j)(ii,jj);
203  data->writeAndTranslate(matrix,(float)data->getOffset(),(float)data->getGain(),H5::PredType::NATIVE_UINT16);
204  }
205 
206  for(unsigned i=0;i<to_store_uchar.size();i++)
207  for(unsigned j=0;j<to_store_uchar[i]->size();j++)
208  {
209  vector<odim::PolarScan*> scans;
210  scans = volume->getScans(to_store_uchar[i]->scan(j).elevation,0.1);
211  shared_ptr<odim::PolarScan> scan;
212  if(scans.size()) scan.reset(scans[0]);
213  else
214  {
215  scan.reset(volume->createScan());
216  scan->setEAngle(to_store_uchar[i]->scan(j).elevation);
217  }
218  unique_ptr<odim::PolarScanData> data(scan->createQuantityData(to_store_uchar[i]->quantity));
219  data->setGain((double)to_store_uchar[i]->scan(j).gain);
220  data->setOffset((double)to_store_uchar[i]->scan(j).offset);
221  data->setNodata((double)to_store_uchar[i]->scan(j).nodata);
222  data->setUndetect((double)to_store_uchar[i]->scan(j).undetect);
223  odim::RayMatrix<float> matrix;
224  matrix.resize(to_store_uchar[i]->scan(j).beam_count,to_store_uchar[i]->scan(j).beam_size);
225  for(unsigned ii=0;ii<to_store_uchar[i]->scan(j).beam_count;ii++)
226  for(unsigned jj=0;jj<to_store_uchar[i]->scan(j).beam_size;jj++)
227  matrix.elem(ii,jj) = to_store_uchar[i]->scan(j)(ii,jj);
228  data->writeAndTranslate(matrix,(float)data->getOffset(),(float)data->getGain(),H5::PredType::NATIVE_UINT8);
229  }
230 
231  for(unsigned i=0;i<to_replace.size();i++)
232  for(unsigned j=0;j<to_replace[i]->size();j++)
233  {
234  vector<odim::PolarScan*> scans;
235  scans = volume->getScans(to_replace[i]->scan(j).elevation,0.1);
236 
237  shared_ptr<odim::PolarScan> scan;
238  if(scans.size()) scan.reset(scans[0]);
239  else
240  {
241 // scan.reset(volume->createScan());
242 // scan->setEAngle(to_store_uchar[i]->scan(j).elevation);
243  ; // da sistemare. in questo caso bisogna far uscire una segnalazione perchè qui non dovremmo proprio finirci.
244  }
245 
246  unique_ptr<odim::PolarScanData> data(scan->getQuantityData(to_replace[i]->quantity));
247  H5::AtomType type = data->getDataType();
248 
249  odim::RayMatrix<double> matrix;
250  matrix.resize(to_replace[i]->scan(j).beam_count,to_replace[i]->scan(j).beam_size);
251  for(unsigned ii=0;ii<to_replace[i]->scan(j).beam_count;ii++)
252  for(unsigned jj=0;jj<to_replace[i]->scan(j).beam_size;jj++)
253  matrix.elem(ii,jj) = to_replace[i]->scan(j)(ii,jj);
254  data->writeAndTranslate(matrix,(float)data->getOffset(),(float)data->getGain(),type);
255  }
256 
257 }
258 
259 
260 void ODIMStorer::storeQuality(const std::string& pathname, const std::string & task, bool RemoveQualityFields)
261 {
262  namespace odim = OdimH5v21;
263  using namespace Radar;
264  using namespace std;
265 
266  shared_ptr<LoadInfo> load_info = make_shared<LoadInfo>();
267  load_info->filename = pathname;
268 
269  unique_ptr<odim::OdimFactory> factory(new odim::OdimFactory());
270  unique_ptr<odim::PolarVolume> volume(factory->openPolarVolume(pathname));
271 
272  for(unsigned i=0;i<to_store_uchar.size();i++)
273  for(unsigned j=0;j<to_store_uchar[i]->size();j++)
274  {
275  vector<odim::PolarScan*> scans;
276  scans = volume->getScans(to_store_uchar[i]->scan(j).elevation,0.1);
277  shared_ptr<odim::PolarScan> scan;
278  if(scans.size()) scan.reset(scans[0]);
279  else
280  {
281  scan.reset(volume->createScan());
282  scan->setEAngle(to_store_uchar[i]->scan(j).elevation);
283  }
284  if (RemoveQualityFields) {
285  int iqc =scan->getQualityCount();
286  for (int iq=iqc-1; iq >=0 ; iq--){
287  scan->removeQuality(iq);
288  }
289  }
290  unique_ptr<odim::OdimQuality> quality(scan->createQuality());
291  quality->getWhat()->set(ATTRIBUTE_WHAT_OFFSET, 0.);
292  quality->getWhat()->set(ATTRIBUTE_WHAT_GAIN, 1.);
293  quality->getHow() ->set(ATTRIBUTE_HOW_TASK, task);
294  odim::RayMatrix<unsigned short> matrix;
295  matrix.resize(to_store_uchar[i]->scan(j).beam_count,to_store_uchar[i]->scan(j).beam_size);
296  for(unsigned ii=0;ii<to_store_uchar[i]->scan(j).beam_count;ii++)
297  for(unsigned jj=0;jj<to_store_uchar[i]->scan(j).beam_size;jj++)
298  matrix.elem(ii,jj) = to_store_uchar[i]->scan(j)(ii,jj);
299  quality->writeQuality(matrix);
300  }
301 
302 }
303 
304 
305 } // namespace volume
306 } // 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:272
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:330
std::string quantity
Odim quantity name.
Definition: volume.h:266
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