4#include <radarlib/radar.hpp>
15using namespace radarelab::algo;
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;
37 unsigned char beam_header[40];
38 if (!in.
fread(beam_header,
sizeof(beam_header)))
41 decode_header_sp20_date_from_name(beam_header, &beam_info, (
char*)in.name());
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);
54 beam_size = beam_info.cell_num;
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);
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]; }
74struct Beams :
public std::vector<sp20::Beam*>
79 unsigned beam_size = 0;
81 Beams(
double elevation) : elevation(elevation) {}
82 Beams(
const Beams&) =
delete;
83 Beams(
const Beams&&) =
delete;
89 Beams& operator=(
const Beams&) =
delete;
92 unsigned min_beam_size()
const
94 if (empty())
return 0;
95 unsigned res = (*this)[0]->beam_size;
97 res = min(res, i->beam_size);
102struct Elevations :
public std::vector<unique_ptr<Beams>>
104 Elevations(
float* elevations,
unsigned el_count)
108 for (
unsigned i = 0; i < el_count; ++i)
109 els.push_back(elevations[i]);
110 std::sort(els.begin(), els.end());
114 emplace_back(
new Beams(el));
117 Beams* get_closest(
double elevation)
119 for (
unsigned i=0; i < size(); ++i)
120 if (elevation >= (at(i)->elevation-0.5) && elevation < (at(i)->elevation+0.5))
128void SP20Loader::load(
const std::string& pathname)
130 namespace odim = OdimH5v21;
131 LOG_CATEGORY(
"Volume");
133 static const float size_cell_by_resolution[]={62.5,125.,250.,500.,1000.,2000.};
135 HD_DBP_SP20_RAW hd_char;
137 shared_ptr<LoadInfo> load_info = make_shared<LoadInfo>();
138 load_info->filename = pathname;
142 File sp20_in(logging_category);
143 sp20_in.open(pathname,
"rb",
"input sp20 file");
145 throw std::runtime_error(
"failed to open sp20 input file");
148 sp20_in.fread(&hd_char,
sizeof(hd_char));
150 HD_DBP_SP20_DECOD hd_file;
151 decode_header_DBP_SP20(&hd_char, &hd_file);
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;
164 Elevations elevations(hd_file.ele, hd_file.num_ele);
168 unique_ptr<sp20::Beam> beam(
new sp20::Beam);
170 if (!beam->read(sp20_in))
173 Beams* beams = elevations.get_closest(beam->beam_info.elevation);
174 if (!beams)
continue;
176 beams->push_back(beam.release());
179 if (elevations.empty())
180 throw std::runtime_error(
"sp20 file has no beams");
183 for (
auto& beams: elevations)
185 unsigned beam_size = beams->min_beam_size();
186 beams->beam_size = beam_size;
189 if (elevations.back()->beam_size == 0)
190 throw std::runtime_error(
"last elevation beam size is 0");
194 for (
int i = elevations.size() - 2; i >= 0; --i)
196 if (elevations[i]->beam_size == 0)
197 elevations[i]->beam_size = elevations[i + 1]->beam_size;
201 for (
unsigned idx = 0; idx < elevations.size(); ++idx)
203 const Beams& beams = *elevations[idx];
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);
212 for (
unsigned i = 0; i < beams.size(); ++i)
213 beam_to_volumes(*beams.at(i), i, beams.beam_size, idx);
218 vol_z->load_info = load_info;
219 if (load_info->declutter_rsp)
220 vol_z->quantity = odim::PRODUCT_QUANTITY_DBZH;
222 vol_z->quantity = odim::PRODUCT_QUANTITY_TH;
223 for (
auto& scan : *vol_z)
227 scan.gain = 80.0 / 255.0;
233 vol_d->load_info = load_info;
234 vol_v->quantity = odim::PRODUCT_QUANTITY_ZDR;
235 for (
auto& scan : *vol_z)
239 scan.gain = 16.0 / 255.0;
245 vol_v->load_info = load_info;
246 vol_v->quantity = odim::PRODUCT_QUANTITY_VRAD;
247 for (
auto& scan : *vol_v)
252 scan.undetect = -49.5;
253 scan.gain = 99.0 / 254;
257 scan.undetect = -16.5;
258 scan.gain = 33.0 / 254;
265 vol_w->load_info = load_info;
266 vol_w->quantity = odim::PRODUCT_QUANTITY_WRAD;
267 for (
auto& scan : *vol_w)
271 scan.gain = 10.0 / 255.0;
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);
281void SP20Loader::beam_to_volumes(
const sp20::Beam& beam,
unsigned az_idx,
unsigned beam_size,
unsigned el_num)
283 const unsigned max_range = min(beam_size, beam.beam_size);
285 if (vol_z && beam.has_z())
288 Eigen::VectorXd dbs(max_range);
289 for (
unsigned i = 0; i < max_range; ++i)
291 PolarScan<double>& scan = vol_z->at(el_num);
292 scan.row(az_idx) = dbs;
297 if (vol_d && beam.has_d())
300 const double range_zdr = 16.;
301 const double min_zdr = -6.;
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;
308 PolarScan<double>& scan = vol_d->at(el_num);
309 scan.row(az_idx) = dbs;
314 if (vol_v && beam.has_v())
317 Eigen::VectorXd ms(max_range);
318 if (beam.beam_info.PRF ==
'S')
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;
326 ms(i) = beam.data_v[i] * range_v / 254.;
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;
334 ms(i) = beam.data_v[i] * range_v / 254.;
336 PolarScan<double>& scan = vol_v->at(el_num);
337 scan.row(az_idx) = ms;
342 if (vol_w && beam.has_w())
345 const double range_sig_v = 10.;
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;
351 PolarScan<double>& scan = vol_w->at(el_num);
352 scan.row(az_idx) = ms;
bool fread(void *buf, size_t size)
Performs a fread on the file, throwing an exception if anything goes wrong.
static unsigned char DBtoBYTE(double DB, double gain=80./255., double offset=-20.)
funzione che converte dB in valore intero tra 0 e 255
static constexpr double BYTEtoDB(unsigned char DBZbyte, double gain=80./255., double offset=-20.)
funzione che converte Z unsigned char in DBZ
Namespace per volume dati.
Eigen::VectorXd elevations_real
Vector of actual elevations for each beam.
Eigen::VectorXd azimuths_real
Vector of actual azimuths for each beam.