Elaboradar 0.1
Caricamento in corso...
Ricerca in corso...
Nessun risultato
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>
9
10#ifdef NEL
11#undef NEL
12#endif
13
14using namespace std;
15using namespace radarelab::algo;
16
18int elev_array[15];
19
20namespace radarelab {
21namespace volume {
22
23namespace sp20 {
24
25struct 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
72namespace {
73
74struct 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
102struct 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
128void 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
281void 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
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
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