Elaboradar 0.1
Caricamento in corso...
Ricerca in corso...
Nessun risultato
odim.cpp
1#include "odim.h"
2#include "logging.h"
3#include "utils.h"
4#include <radarlib/radar.hpp>
5#include <memory>
6
7using namespace std;
8using namespace OdimH5v21;
9
10namespace {
11
12unsigned 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
26namespace radarelab {
27namespace volume {
28
29void ODIMLoader::request_quantity(const std::string& name, Scans<double>* volume)
30{
31 to_load.insert(make_pair(name, volume));
32}
33
34void 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
163void 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
278void 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
std::string quantity
Odim quantity name.
Definition volume.h:268
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
Sequence of PolarScans which can have a different beam count for each elevation.
Definition volume.h:264
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
std::map< std::string, Scans< double > * > to_load
Map used to specify quantity to be loaded.
Definition loader.h:26
void request_quantity(const std::string &name, Scans< double > *volume)
Define a request - Fill to_load attribute
Definition odim.cpp:29
void load(const std::string &pathname)
Load method.
Definition odim.cpp:34