11#include <radarelab/algo/vpr.h>
16#include "cartproducts.h"
21#include <radarlib/radar.hpp>
49#define OVERBLOCKING 51
50#define SOGLIA_TOP 45.0
57#define LIMITE_ANAP 240
60#define CONV_RAD 360./4096.*DTOR
64using namespace radarelab::algo;
76 static const int NSCAN = 6;
84 const double radius = 6378137.0;
85 const double kea = 4. / 3. * radius;
87 for (
unsigned iel = 0; iel < vol.size(); ++iel)
89 const double elev = vol.
scan(iel).elevation;
90 const double cell_size = vol.
scan(iel).cell_size;
92 for (
unsigned ibin = 0; ibin < cols(); ++ibin)
94 double range = (ibin + 0.5) * cell_size;
95 (*this)(iel, ibin) = ::sqrt(::pow(range, 2.) + ::pow(kea, 2.) + 2 * kea * range * ::sin(DTOR * elev)) - kea;
100 void load_hray(Assets& assets)
103 dtrs = assets.read_file_hray([
this](
unsigned el,
unsigned bin,
double value) {
104 if (el >= rows() || bin >= cols())
return;
105 (*this)(el, bin) = value;
108 void load_hray_inf(Assets& assets)
111 dtrs = assets.read_file_hray_inf([
this](
unsigned el,
unsigned bin,
double value) {
112 if (el >= rows() || bin >= cols())
return;
113 (*this)(el, bin) = value;
120 : MyMAX_BIN(max_bin), cfg(cfg), site(site), assets(cfg),
121 do_medium(medium), volume(volume), SD_Z6(volume.beam_count), cil(volume, 0, RES_HOR_CIL, RES_VERT_CIL),
122 dbz(volume), flag_vpr(volume, 0),
123 first_level(NUM_AZ_X_PPI, MyMAX_BIN), first_level_static(NUM_AZ_X_PPI, MyMAX_BIN),
124 bb_first_level(NUM_AZ_X_PPI, 1024), beam_blocking(NUM_AZ_X_PPI, 1024),
125 anaprop(volume), dem(NUM_AZ_X_PPI, 1024),
127 top(volume.beam_count, volume.max_beam_size())
134 struct tm* tempo = gmtime(&Time);
136 strftime(
date, 20,
"%Y%m%d%H%M", tempo);
139 algo::compute_top(
volume, SOGLIA_TOP,
top);
155 LOG_CATEGORY(
"radar.io");
156 LOG_INFO(
"Reading %s for site %s", nome_file,
site.
name.c_str());
163 loader.vol_z = &z_volume;
164 loader.vol_w = &w_volume;
165 loader.vol_v = &v_volume;
166 loader.load(nome_file);
176 for (
unsigned i = 0; i < z_volume.size(); ++i) {
178 * v_volume.at(i).gain + v_volume.at(i).
offset;
179 algo::Cleaner::clean(z_volume.at(i), w_volume.at(i), v_volume.at(i), bin_wind_magic_number);
183 algo::azimuthresample::MaxOfClosest<double> resampler;
184 resampler.resample_volume(z_volume,
volume, 1.0);
193 LOG_CATEGORY(
"radar.io");
194 namespace odim = OdimH5v21;
195 LOG_INFO(
"Reading %s for site %s", nome_file,
site.
name.c_str());
211 string radar_name =
site.
name.c_str();
212 bool init_sqi =
false;
226 loader.
load(nome_file);
229 if (dbzh_volume.empty() && th_volume.empty())
231 LOG_ERROR(
"neither DBZH nor TH were found in %s", nome_file);
232 throw runtime_error(
"neither DBZH nor TH were found");
238 if(!i.second->empty() ) i.second->normalize_elevations(elev_array);
241 if (!dbzh_volume.empty()) {
242 LOG_WARN(
" DBZH found");
243 z_volume = &dbzh_volume;
246 LOG_WARN(
"no DBZH found: using TH");
247 z_volume = &th_volume;
251 if (do_clean && !w_volume.empty() && !v_volume.empty())
253 if(sqi_volume.empty()) init_sqi =
true;
255 if (zdr_volume.empty())
260 for (
unsigned i = 0; i < z_volume->size(); ++i){
262 full_volume_cleanID.
append_scan(z_volume->at(i).beam_count,z_volume->at(i).beam_size,z_volume->at(i).elevation, z_volume->at(i).cell_size, 0);
263 radarelab::algo::Cleaner::evaluateCleanID(z_volume->at(i), w_volume.at(i), v_volume.at(i),full_volume_cleanID.at(i),i);
264 for (
unsigned ibeam=0;ibeam<z_volume->at(i).beam_count; ibeam++)
265 for (
unsigned j=0; j<z_volume->at(i).beam_size; j++){
266 if (full_volume_cleanID.at(i)(ibeam,j) != 0) z_volume->at(i)(ibeam,j)=z_volume->at(i).undetect;
270 cout<<
"applico logica fuzzy"<<endl;
276 unsigned last = z_volume->size() -1;
277 for (
unsigned i = 0; i < z_volume->size(); ++i){
281 full_volume_cleanID.
append_scan(z_volume->at(i).beam_count,z_volume->at(i).beam_size,z_volume->at(i).elevation, z_volume->at(i).cell_size);
282 full_volume_diffprob.
append_scan(z_volume->at(i).beam_count,z_volume->at(i).beam_size,z_volume->at(i).elevation, z_volume->at(i).cell_size);
288 sqi_volume.
append_scan(z_volume->at(i).beam_count,z_volume->at(i).beam_size,z_volume->at(i).elevation, z_volume->at(i).cell_size);
289 sqi_volume.at(i).setZero();
295 Input.push_back(z_volume->at(i));
296 Input2.push_back(z_volume->at(i+1));
297 radarelab::volume::textureVD(Input, Input2, Texture,
true);
298 Texture.at(0).nodata=65535.;
299 Texture.at(0).undetect=0.;
305 Texture.
append_scan(z_volume->at(i).beam_count,z_volume->at(i).beam_size,z_volume->at(i).elevation, z_volume->at(i).cell_size);
306 Texture.at(0).setZero();
307 Texture.at(0).nodata=65535.;
308 Texture.at(0).undetect=0.;
315 algo::Cleaner::evaluateClassID(z_volume->at(i), w_volume.at(i), v_volume.at(i), zdr_volume.at(i), rhohv_volume.at(i), sqi_volume.at(i), snr_volume.at(i), Texture.at(0), full_volume_cleanID.at(i), full_volume_diffprob.at(0), v_volume.at(i).undetect , radar_name, fuzzypath, i,
true);
317 double new_value=z_volume->at(0).nodata;
319 if (set_undetect) new_value=z_volume->at(i).undetect;
321 for (
unsigned ii = 0; ii < z_volume->at(i).beam_count; ++ii)
322 for (
unsigned ib = 0; ib < z_volume->at(i).beam_size; ++ib) {
324 if(full_volume_cleanID.at(i)(ii,ib) ) z_volume->at(i)(ii,ib)= new_value;
337 algo::azimuthresample::MaxOfClosest<double> resampler;
338 resampler.resample_volume(*z_volume,
volume, 1.0);
339 cout<<
"resampler fatto!!!"<<endl;
405 for(
unsigned i=0; i<NUM_AZ_X_PPI; i++)
406 for(
unsigned k=0; k<
volume[0].beam_size; k++)
411 for(
unsigned l=0; l<=el_inf; l++)
414 if (k >=
volume[l].beam_size)
continue;
415 if (k <
volume[el_inf].beam_size)
418 volume[l].set(i, k, MISSING_DB);
426 LOG_INFO(
"declutter_anaprop completed with static declutter");
444 std::vector <double> above_0 (4,0);
445 std::vector <double> above_15 (4,0);
446 std::vector <double> above_30 (4,0);
447 std::vector <double> above_40 (4,0);
448 for(
unsigned el =0; el <4; el++)
450 for (
unsigned k=80; k <
volume[el].beam_size; k ++)
451 if (
volume[el].get(iaz,k) > 40.){
456 }
else if (
volume[el].get(iaz,k) > 30.){
460 }
else if (
volume[el].get(iaz,k) > 15.){
463 }
else if (
volume[el].get(iaz,k) > 0.){
470 if ( above_15[2]/above_15[0] >= 0.025){
471 if (above_0[1]/above_0[0] >= 0.6 && above_30[2]/above_15[2] <0.15 && above_0[1] >=50000){
472 anaprop.conf_texture_threshold = 5.;
473 LOG_WARN(
"TEXTURE THRESHOLD USED %4.1f -- 0. %6d %6d %6d %6d -- 15. %6d %6d %6d %6d -- 30. %6d %6d %6d %6d -- 40. %6d %6d %6d %6d",
anaprop.conf_texture_threshold, (
int)above_0[0], (
int)above_0[1], (
int)above_0[2], (
int)above_0[3], (
int)above_15[0], (
int)above_15[1], (
int)above_15[2], (
int)above_15[3], (
int)above_30[0], (
int)above_30[1], (
int)above_30[2], (
int)above_30[3], (
int)above_40[0], (
int)above_40[1], (
int)above_40[2], (
int)above_40[3] );
478 LOG_WARN(
"THUNDERSTORM %4.1f -- 0. %6d %6d %6d %6d -- 15. %6d %6d %6d %6d -- 30. %6d %6d %6d %6d -- 40. %6d %6d %6d %6d", -9.9, (
int)above_0[0], (
int)above_0[1], (
int)above_0[2], (
int)above_0[3], (
int)above_15[0], (
int)above_15[1], (
int)above_15[2], (
int)above_15[3], (
int)above_30[0], (
int)above_30[1], (
int)above_30[2], (
int)above_30[3], (
int)above_40[0], (
int)above_40[1], (
int)above_40[2], (
int)above_40[3] );
482 LOG_WARN(
"TEXTURE THRESHOLD USED %4.1f -- 0. %6d %6d %6d %6d -- 15. %6d %6d %6d %6d -- 30. %6d %6d %6d %6d -- 40. %6d %6d %6d %6d",
anaprop.conf_texture_threshold, (
int)above_0[0], (
int)above_0[1], (
int)above_0[2], (
int)above_0[3], (
int)above_15[0], (
int)above_15[1], (
int)above_15[2], (
int)above_15[3], (
int)above_30[0], (
int)above_30[1], (
int)above_30[2], (
int)above_30[3], (
int)above_40[0], (
int)above_40[1], (
int)above_40[2], (
int)above_40[3] );
485 LOG_INFO(
"declutter_anaprop completed with anaprop");
490 LOG_WARN(
"declutter_anaprop completed without doing anything");
500 LOG_INFO(
"elabora_Dato completata");
515 LOG_INFO(
"Letta mappa statica");
561 for (
unsigned i=NUM_AZ_X_PPI; i<800; i++)
563 for (
unsigned k=i-1; k<i+2; k++)
564 if(
first_level(i%NUM_AZ_X_PPI, j) < first_level_tmp(k%NUM_AZ_X_PPI, j))
565 first_level(i%NUM_AZ_X_PPI, j)=first_level_tmp(k%NUM_AZ_X_PPI, j);
573 static const int DIM1_ST = 16;
574 static const int DIM2_ST = 13;
577 static const int N_MIN_BIN = 500;
580 unsigned char statistica[DIM1_ST][DIM2_ST];
581 unsigned char statistica_bl[DIM1_ST][DIM2_ST];
582 unsigned char statistica_el[DIM1_ST][DIM2_ST];
584 memset(statistica,255,DIM1_ST*DIM2_ST);
585 memset(statistica_bl,255,DIM1_ST*DIM2_ST);
586 memset(statistica_el,255,DIM1_ST*DIM2_ST);
588 for(az=0; az<DIM1_ST; az++)
589 for(ran=0; ran<DIM2_ST; ran++){
590 if (grid_stats.count(az, ran) >= N_MIN_BIN)
592 statistica[az][ran] = grid_stats.perc_anap(az, ran);
593 statistica_bl[az][ran] = grid_stats.perc_bloc(az, ran);
594 statistica_el[az][ran] = grid_stats.perc_elev(az, ran);
601 fwrite(
date,12,1,f_stat);
602 fwrite(statistica,DIM1_ST*DIM2_ST,1,f_stat);
607 fwrite(
date,12,1,f_stat);
608 fwrite(statistica_bl,DIM1_ST*DIM2_ST,1,f_stat);
613 fwrite(
date,12,1,f_stat);
614 fwrite(statistica_el,DIM1_ST*DIM2_ST,1,f_stat);
638 LOG_DEBUG(
"start caratterizzo_volume");
641 hray_inf.load_hray_inf(
assets);
643 cout<<
"sono in caratterizzo volume"<<endl;
663 for (
unsigned l=0; l<
volume.size(); l++)
665 const auto& scan =
volume[l];
666 for (
int i=0; i<NUM_AZ_X_PPI; i++)
668 const double elevaz = scan.elevations_real(i);
673 for (
unsigned k=0; k<scan.beam_size; k++)
675 double sample = scan.get(i, k);
678 dist = k * scan.cell_size + scan.cell_size / 2.;
699 if (l<
volume.size() -1 ) {
701 dh = hray_inf(l + 1, k) - hray_inf(l, k);
708 if (l <
anaprop.elev_fin(i, k)) {
711 }
else if (l ==
anaprop.elev_fin(i, k)) {
712 cl=
anaprop.dato_corrotto(i, k);
714 }
else if (l >
anaprop.elev_fin(i, k)) {
722 if (l <
anaprop.elev_fin(i, k)) {
730 qual.
scan(l).set(i, k, (
unsigned char)(
func_q_Z(cl,bb,dist,drrs,hray_inf.dtrs,dh,dhst,PIA)*100));
737 if(cl==0 && bb<BBMAX_VPR )
738 flag_vpr.
scan(l).set(i, k, 1);
744 LOG_DEBUG(
"End caratterizzo_volume");
750 LOG_CATEGORY(
"radar.class");
770 hbbb=(hmax-600.)/1000.;
771 htbb=(hmax+600.)/1000.;
775 LOG_INFO(
"leggo 0termico per class da file %s",getenv(
"FILE_ZERO_TERMICO"));
781 htbb=zeroterm/1000. + 0.4;
782 hbbb=zeroterm/1000. - 1.0;
784 LOG_ERROR(
"non ho trovato il file dello zero termico");
785 LOG_INFO(
"attenzione, non ho trovat zero termico ne da vpr ne da radiosondaggio");
794 LOG_INFO(
"calcolati livelli sopra e sotto bright band hbbb=%f htbb=%f",
hbbb,
htbb);
799 LOG_DEBUG (
"Matrice cilindrica Naz %3d Nrange %4d Nheight %4d", cil.
slices.size(), cil.x_size, cil.z_size);
803 viz.classifico_VIZ();
810 steiner.calcolo_background();
811 steiner.classifico_STEINER();
819 for (
unsigned j = 0; j < NUM_AZ_X_PPI; ++j)
820 for (
unsigned k = 0; k < steiner.max_bin; ++k)
821 if (
cum_bac.
anaprop.quota(j, k) <
hbbb*1000. && steiner.conv_STEINER(j, k) == viz.conv_VIZ(j, k) && steiner.conv_STEINER(j, k) > 0)
822 conv(j,k) = steiner.conv_STEINER(j, k);
824 if (steiner.conv_STEINER(j, k) == viz.conv_VIZ(j, k) && steiner.conv_STEINER(j, k) > 0 && viz.stratiform(j, k) < 1)
825 conv(j,k) = viz.conv_VIZ(j, k);
849 LOG_CATEGORY(
"radar.vpr");
851 LOG_DEBUG (
" modalita %d", MOD_VPR);
859 for (
unsigned i=0; i<vpr0.size(); i++)
860 LOG_DEBUG (
" Profilo vecchio - livello %2d valore %6.2f",i,vpr0.val[i]);
862 LOG_INFO(
"gap %li",
gap);
870 if (inst_vpr.success)
872 vpr = combine_profiles(vpr0, inst_vpr.vpr, inst_vpr.cv, inst_vpr.ct);
878 if (inst_vpr.success)
890 LOG_INFO(
" livmin %i",
livmin.livmin);
895 this->livmin =
livmin.livmin;
900 for (
unsigned i=0; i<
vpr.size(); i++) LOG_DEBUG (
" Profilo nuovo - livello %2d valore %6.2f",i,
vpr.val[i]);
906#include <radarelab/vpr_par.h>
908 LOG_CATEGORY(
"radar.vpr");
910 int heating = cum_bac.assets.read_vpr_heating();
920 heating = heating - (gap-1) + 1;
921 if (heating>=WARM) heating=MEMORY;
923 if (heating<0) heating=0;
926 cum_bac.assets.write_vpr_heating(heating);
929 LOG_INFO(
"gap %li heating %i",gap,heating);
939 File file(logging_category);
940 file.
open_from_env(
"VPR_ARCH",
"wt",
"ultimo vpr in dBZ per il plot");
942 fprintf(file,
" QUOTA DBZ AREA PRECI(KM^2/1000)\n" );
943 for (
unsigned ilay=0; ilay<NMAXLAYER; ilay++){
944 if (
vpr.val[ilay]> 0.001 ) {
946 fprintf(file,
" %i %10.3f %li\n", ilay*TCK_VPR+TCK_VPR/2, vpr_dbz,
vpr.area[ilay]);
949 fprintf(file,
" %i %10.3f %li\n", ilay*TCK_VPR+TCK_VPR/2, NODATAVPR,
vpr.area[ilay]);
978#include <radarelab/vpr_par.h>
981 LOG_CATEGORY(
"radar.vpr");
983 int ilray,ilref,ilay2,ier_ana,snow,strat;
984 float corr,vpr_liq,vpr_hray,hbin,hliq;
995 ier_max=trovo_hvprmax(&hvprmax);
996 ier_ana=analyse_VPR(&vpr_liq,&snow,&hliq);
997 LOG_INFO(
"ier_analisi %i",ier_ana) ;
1000 if (ier_ana)
return 1;
1002 LOG_INFO(
"altezza bright band %i",hvprmax);
1003 LOG_INFO(
"CORREGGO VPR");
1006 for (
unsigned i=0; i<NUM_AZ_X_PPI; i++){
1007 for (
unsigned k=0; k<cum_bac.volume[0].beam_size; k++){
1011 hbin=(float)cum_bac.anaprop.quota(i, k);
1014 if (snow) neve(i, k)=1;
1016 if (cum_bac.do_class)
1018 if (conv(i,k) >= CONV_VAL){
1023 if (cum_bac.volume[0].get(i, k) > THR_CORR && hbin > hliq && strat){
1026 ilray=(hbin>=livmin)?(floor(hbin/TCK_VPR)):(floor(livmin/TCK_VPR));
1027 if (ilray>= NMAXLAYER ) ilray=NMAXLAYER-2;
1030 if ((
int)hbin%TCK_VPR > TCK_VPR/2) ilay2=ilray+1;
1032 if (ilay2< floor(livmin/TCK_VPR)) ilay2=floor(livmin/TCK_VPR);
1037 ilref=(cum_bac.dem(i, k)>livmin)?(floor(cum_bac.dem(i, k)/TCK_VPR)):(floor(livmin/TCK_VPR));
1040 if (vpr.val[ilref] > 0 && vpr.val[ilray] > 0 ){
1042 vpr_hray=vpr.val[ilray]+((vpr.val[ilray]-vpr.val[ilay2])/(ilray*TCK_VPR-TCK_VPR/2-ilay2*TCK_VPR))*(hbin-ilray*TCK_VPR-TCK_VPR/2);
1045 if (cum_bac.dem(i, k)> hvprmax+HALF_BB-TCK_VPR || snow){
1065 corr=cum_bac.dbz.RtoDBZ(vpr.val[ilref])-cum_bac.dbz.RtoDBZ(vpr_hray);
1067 cum_bac.volume[0].set(i, k, cum_bac.dbz.DBZ_snow(cum_bac.volume[0].get(i, k)));
1071 corr = cum_bac.dbz.RtoDBZ_class(vpr_liq) - cum_bac.dbz.RtoDBZ_class(vpr_hray);
1074 if (corr>MAX_CORR) corr=MAX_CORR;
1075 if (hbin<hvprmax && corr>0.) corr=0;
1078 double corrected = cum_bac.volume[0].get(i, k) + corr;
1079 if (corrected > MAXVAL_DB)
1080 cum_bac.volume[0].set(i, k, MAXVAL_DB);
1081 else if ( corrected < MINVAL_DB)
1082 cum_bac.volume[0].set(i, k, MINVAL_DB);
1084 cum_bac.volume[0].set(i, k, corrected);
1086 corr_polar(i, k)=(
unsigned char)(corr)+128;
1102 int i,imax,istart,foundlivmax;
1103 float h0start,peak,soglia;
1108 LOG_DEBUG(
"trovo hvprmax a partire da 400 m sotto lo zero dell'adiabatica secca");
1110 istart=h0start/TCK_VPR -2;
1112 LOG_DEBUG(
"t_ground h0start istart %f %f %i",
t_ground,h0start,istart);
1115 LOG_DEBUG(
"trovo hvprmax a partire da livmin");
1128 soglia = DBZ::DBZtoR(THR_VPR,200,1.6);
1131 LOG_DEBUG(
" istart %d low %6.2f up %6.2f soglia %6.2f peak %6.2f imax %d", istart,
vpr.val[istart] ,
vpr.val[istart+4], soglia, peak, imax);
1132 if (
vpr.val[istart] >soglia &&
vpr.val[istart+4] > soglia){
1133 peak=10*log10(
vpr.val[istart]/
vpr.val[istart+4]);
1134 LOG_DEBUG(
"peak1 = %f",peak);
1137 if(peak> MIN_PEAK_VPR){
1140 LOG_DEBUG(
"il primo punto soddisfa le condizioni di picco");
1142 for (i=istart+1;i<NMAXLAYER-4;i++)
1144 if (
vpr.val[i] <soglia ||
vpr.val[i+4] < soglia)
break;
1145 peak=10*log10(
vpr.val[i]/
vpr.val[i+4]);
1146 if (
vpr.val[i]>
vpr.val[i-1] && peak> MIN_PEAK_VPR )
1151 LOG_DEBUG(
" low %6.2f up %6.2f soglia %6.2f peak %6.2f imax %d",
vpr.val[i] ,
vpr.val[i+4], soglia, peak, imax);
1154 if ( imax > INODATA ){
1156 peak=10*log10(
vpr.val[imax]/
vpr.val[imax+4]);
1157 *hmax=imax*TCK_VPR+TCK_VPR/2;
1158 LOG_DEBUG(
"trovato ilaymax %i %i",*hmax,imax);
1159 LOG_DEBUG(
" picco in dbR %f",peak);
1162 LOG_DEBUG(
"exit status trovo_hvprmax %i",foundlivmax);
1163 return (foundlivmax);
1183 int ier=1,ier_ana=0,liv0;
1184 char date[]=
"000000000000";
1192 int tipo_profilo=-1;
1193 float v600sottobb=NODATAVPR;
1194 float v1000=NODATAVPR;
1195 float v1500=NODATAVPR;
1196 float vliq=NODATAVPR;
1197 float vhliquid=NODATAVPR;
1198 float vprmax=NODATAVPR;
1206 LOG_WARN(
"non ho T,...");
1209 LOG_ERROR(
" non ho trovato hvprmax, nè T, esco");
1219 LOG_ERROR(
" temperatura alta e non ho trovato hvprmax, esco");
1231 LOG_INFO(
" temperatura da scioglimento e massimo in quota");
1235 LOG_ERROR(
" temperatura da scioglimento ma superiore a temperatura max neve e non ho trovato hvprmax, esco");
1240 if (
hvprmax > liv0) LOG_INFO(
" il livello %i è sotto la Bright band, ma T bassa interpolo",
livmin);
1241 else LOG_INFO(
" il livello %i potrebbe essere dentro la Bright Band, interpolo",
livmin);
1249 LOG_INFO(
" temperatura da neve o scioglimento e massimo in quota");
1253 LOG_INFO(
" temperatura da neve o scioglimento e massimo non trovato,neve , non interpolo");
1262 InterpolaVPR_GSL iv;
1272 LOG_INFO(
" interpolazione fallita");
1273 switch (tipo_profilo)
1284 *vpr_liq=
vpr.val[(
hvprmax+1000)/TCK_VPR]*2.15;
1290 LOG_INFO(
" interpolazione eseguita con successo");
1293 const char* vpr_arch = getenv(
"VPR_ARCH");
1294 if (!vpr_arch)
throw runtime_error(
"VPR_ARCH is not defined");
1295 string fname(vpr_arch);
1297 File file(logging_category);
1298 file.
open(fname,
"wt",
"vpr interpolato");
1299 for (
unsigned i = 0; i < NMAXLAYER; ++i)
1300 fprintf(file,
" %f \n",
cum_bac.
dbz.RtoDBZ(iv.vpr_int[i]));
1303 if (tipo_profilo == 2 ) {
1304 *hliq=(iv.E-2.1*iv.G)*1000.;
1308 *vpr_liq=
vpr.val[(
hvprmax+1000)/TCK_VPR]*2.15;
1311 *hliq=(iv.E-2.1*iv.G)*1000.;
1314 *vpr_liq=
vpr.val[(int)(*hliq/TCK_VPR)];
1318 if (*hliq<0) *hliq=0;
1333 LOG_INFO(
"TIPO_PROFILO= %i vpr_liq %f hliq %f", tipo_profilo, *vpr_liq,*hliq );
1341 tempo = gmtime(&Time);
1342 sprintf(date,
"%04d%02d%02d%02d%02d",tempo->tm_year+1900, tempo->tm_mon+1,
1343 tempo->tm_mday,tempo->tm_hour, tempo->tm_min);
1352 if ((
hvprmax+1000)/TCK_VPR < NMAXLAYER )
1354 if ((
hvprmax+1500)/TCK_VPR < NMAXLAYER )
1359 if (FILE* test_vpr=fopen(getenv(
"TEST_VPR"),
"a+"))
1361 fprintf(test_vpr,
"%s %i %i -1 %f %f %f %f %f %f %f %f %f %f %f %f %f %f \n",date,
hvprmax,tipo_profilo,iv.chisqfin,*hliq,vliq,vhliquid,v600sottobb,v1000+6,v1500+6,vprmax,iv.rmsefin,iv.B,iv.E,iv.G,iv.C,iv.F);
1370 LOG_INFO(
"fatta scrittura hmax vpr = %d",
hvprmax);
1377 for (
unsigned i=0; i<NUM_AZ_X_PPI; i++)
1378 for (
unsigned k=0; k<
volume[0].beam_size; k++)
1412void CUM_BAC::generate_maps(CartProducts& products)
1415 LOG_INFO(
"Scrittura File Precipitazione 1X1\n");
1418 std::function<
unsigned char(
const vector<double>&)> convert = [
this](
const vector<double>& samples) {
1425 for (
const auto& s: samples)
1429 if (res == 0)
return (
unsigned char)1;
1432 products.scaled.to_cart_average(
volume[0], convert, products.z_out);
1434 std::function<
unsigned char(
unsigned,
unsigned)> assign_cart =
1435 [
this](
unsigned azimuth,
unsigned range) {
1438 return max(sample, (
unsigned char)1);
1440 products.scaled.to_cart(assign_cart, products.z_out);
1443 std::function<
unsigned char(
unsigned,
unsigned)> assign_cart =
1444 [
this](
unsigned azimuth,
unsigned range) {
1447 return max(sample, (
unsigned char)1);
1449 products.fullres.to_cart(assign_cart, products.z_fr);
1451 products.scaled.to_cart(
top, products.top_1x1);
1454 const auto& elev_fin =
anaprop.elev_fin;
1455 const auto& quota =
anaprop.quota;
1457 std::function<
unsigned char(
unsigned,
unsigned)> assign_qual =
1458 [
this, &elev_fin](
unsigned azimuth,
unsigned range) {
1459 const auto& el = elev_fin(azimuth, range);
1460 if (range >=
volume[el].beam_size)
1461 return (
unsigned char)0;
1462 return qual.
scan(el).get(azimuth, range);
1464 products.scaled.to_cart(assign_qual, products.qual_Z_1x1);
1466 std::function<
unsigned char(
unsigned,
unsigned)> assign_quota =
1467 ["a](
unsigned azimuth,
unsigned range) {
1468 return 128 + round(quota(azimuth, range) / 100.0);
1470 products.scaled.to_cart(assign_quota, products.quota_1x1);
1471 products.scaled.to_cart(
anaprop.dato_corrotto, products.dato_corr_1x1);
1473 std::function<
unsigned char(
unsigned,
unsigned)> assign_elev_fin = [&elev_fin](
unsigned azimuth,
unsigned range) {
1474 return elev_fin(azimuth, range);
1476 products.scaled.to_cart(assign_elev_fin, products.elev_fin_1x1);
1478 products.scaled.to_cart(
beam_blocking, products.beam_blocking_1x1);
1484 std::function<
unsigned char(
unsigned,
unsigned)> assign = [&neve](
unsigned azimuth,
unsigned range) {
1485 return neve(azimuth, range) ? 0 : 1;
1487 products.scaled.to_cart(assign, products.neve_1x1);
1492 products.scaled.to_cart(
calcolo_vpr->conv, products.conv_1x1);
1500 for (
unsigned int i=0; i<
SD_Z6.size(); i++){
1501 SC_SD.creo_cart(
SD_Z6, i);
1502 std::ostringstream oss;
1504 SC_SD.write_out(
assets,oss.str());
1514 inst_vpr(cum_bac.volume, cum_bac.qual, cum_bac.flag_vpr, cum_bac.site.vpr_iaz_min, cum_bac.site.vpr_iaz_max),
1515 conv(NUM_AZ_X_PPI, cum_bac.volume.max_beam_size(),0),
1516 corr_polar(NUM_AZ_X_PPI, cum_bac.volume.max_beam_size()),
1517 neve(NUM_AZ_X_PPI, cum_bac.volume.max_beam_size())
1519 logging_category = log4c_category_get(
"radar.vpr");
1534CalcoloVPR::~CalcoloVPR()
1541 printf (
"calcolo VPR \n") ;
1550 LOG_INFO(
"fatta func vpr %s", inst_vpr.success ?
"ok" :
"errore");
1551 for (
unsigned i=0; i<inst_vpr.vpr.size(); i++) LOG_DEBUG (
" Profilo istantaneo - livello %2d valore %6.2f",i,inst_vpr.vpr.val[i]);
1557 LOG_INFO (
"exit status calcolo VPR istantaneo: %s", inst_vpr.success ?
"ok":
"fallito") ;
1558 LOG_INFO(
"exit status combinaprofili: (1--fallito 0--ok) %i ",ier_comb) ;
1564 printf (
"heating %i \n",
heating);
1565 LOG_INFO(
"ier_comb %i", ier_comb);
1567 if (inst_vpr.success)
1572 if (!ier_comb &&
heating >= WARM){
1575 LOG_INFO(
"exit status correggo vpr: (1--fallito 0--ok) %i",ier) ;
1580 if ( ! ier && inst_vpr.success)
1588 Image<double> azimut;
1589 Image<double> range;
1591 CartData(
int max_bin=512)
1592 : azimut(max_bin), range(max_bin)
1594 for(
int i=0; i<max_bin; i++)
1595 for(
int j=0; j<max_bin; j++)
1597 range(i, j) = hypot(i+.5,j+.5) ;
1598 azimut(i, j) = 90. - atan((j+.5)/(i+.5)) * M_1_PI*180.;
1612 LOG_CATEGORY(
"radar.singlecart");
1616 int x,y,iaz,az_min,az_max;
1620 for(
unsigned i=0; i<
max_bin *2; i++)
1621 for(
unsigned j=0; j<
max_bin *2; j++)
1622 cart(i, j) = MISSING;
1624 LOG_INFO(
"Creo_cart - %u",
max_bin);
1626 for(
unsigned quad=0; quad<4; quad++)
1627 for(
unsigned i=0; i<
max_bin; i++)
1628 for(
unsigned j=0; j<
max_bin; j++)
1630 unsigned irange = (unsigned)round(cd.range(i, j));
1638 az = cd.azimut(i, j);
1643 az = cd.azimut(i, j) + 90.;
1648 az = cd.azimut(i, j) + 180.;
1653 az = cd.azimut(i, j)+270.;
1657 az_min = (int)((az - .45)/.9);
1658 az_max = ceil((az + .45)/.9);
1663 az_min = az_min + NUM_AZ_X_PPI;
1664 az_max = az_max + NUM_AZ_X_PPI;
1666 for(iaz = az_min; iaz<az_max; iaz++){
1668 unsigned char sample = 0;
1669 if (irange < volume[el_index].beam_size)
1670 sample = max((
unsigned char) (volume[el_index].get(iaz%NUM_AZ_X_PPI, irange)), (
unsigned char)1);
1671 if(
cart(y, x) <= sample)
cart(y, x) = sample;
1678 if (getenv(
"DIR_DEBUG") == NULL)
return;
1691 tempo=gmtime(&clock);
1693 return asctime(tempo);
1698 static time_t time_tot = 0,time1 = 0,time2 = 0;
1699 LOG_CATEGORY(
"radar.timing");
1705 time2 = time(&time2);
1707 LOG_INFO(
"tempo parziale %ld ---- totale %ld", time2-time1, time2-time_tot);
void load_first_level_bb_bloc(radarelab::Matrix2D< unsigned char > &matrix)
Open the first level elevation BB bloc file.
void write_vpr0(const radarelab::algo::VPR &vpr)
Write in $VPR0_FILE the vpr calculated.
bool read_0term(float &zeroterm)
Read $FILE_ZERO_TERMICO.
void write_vpr_hmax(int hvprmax)
write in $VPR_HMAX the vpr peak's height.
void load_dem(radarelab::Matrix2D< float > &matrix)
Open the dem file.
void configure(const Site &site, time_t acq_time)
Configure asset lookup with the given details.
float read_t_ground() const
fornisce temperatura al suolo, da lettura file esterno
void write_gdal_image(const radarelab::Matrix2D< T > &image, const char *dir_env_var, const char *name, const char *format)
Write a graphic image with gdal.
long int read_profile_gap() const
Read the gap between the time in $LAST_VPR and the current acquisition time.
bool find_vpr0(const radarelab::algo::DBZ &dbz, radarelab::algo::VPR &vpr0, long int &gap)
Read the gap and the vpr0, and if vpr0 is not found, look it up among the archived VPRs.
void load_first_level(radarelab::Matrix2D< unsigned char > &matrix)
Open the first level file.
void load_first_level_bb_el(radarelab::Matrix2D< unsigned char > &matrix)
Open the first level elevation BB el file.
int read_vpr_hmax()
Read in $VPR_HMAX the vpr peak's height.
void write_last_vpr()
Write the acquisition time in $LAST_VPR file.
Finds resources, like data files, used by the program.
bool do_declutter
use only static declutter map
bool do_devel
Produce additional output.
void declutter_anaprop()
funzione che elabora il dato radar rimuovendo anaprop e beam blocking
bool do_readStaticMap
Read Static clutter map.
radarelab::algo::Anaprop< double > anaprop
Oggetto per correzione ANAPRO.
radarelab::PolarScan< float > dem
dem in coordinate azimut range
CalcoloVPR * calcolo_vpr
Oggetto per calcolare e correggere con VPR.
radarelab::Volume< unsigned char > qual
qualita volume polare
void want_vpr()
Call this just after creating the CUM_BAC object, to signal that VPR should also be computed.
radarelab::PolarScan< unsigned char > first_level_static
mappa statica
void leggo_first_level()
funzione che legge la mappa statica e la mappa di elevazioni da beam blocking e le condensa in un uni...
log4c_category_t * logging_category
logging category
bool do_class
Convective-stratiform classification.
radarelab::Volume< double > & volume
Set to Z undetect value the Zpixels classified as non-meteo echoes.
radarelab::PolarScan< unsigned char > bb_first_level
mappa di elevazioni da beam blocking (input)
static void read_sp20_volume(radarelab::Volume< double > &volume, const Site &site, const char *nome_file, bool do_clean=false, bool do_medium=false)
Read from SP20 data file.
const Site & site
site information object
radarelab::PolarScan< unsigned char > first_level
mappa dinamica complessiva
char date[20]
Acquisition date.
radarelab::Volume< double > SD_Z6
Polar volume of standard deviation of reflectivity over 6 km length.
CUM_BAC(radarelab::Volume< double > &volume, const Config &cfg, const Site &site, bool medium=false, unsigned max_bin=512)
Constructor.
bool do_anaprop
anaprop correction
bool do_bloccorr
bloccorrection
static void read_odim_volume(radarelab::Volume< double > &volume, const Site &site, const char *nome_file, char *fuzzypath, bool do_clean=false, bool do_medium=false, bool set_undetect=false)
Read from ODIM data file.
void vpr_class()
Esegue tutta la catena vpr (e classificazione) se richiesta.
bool do_medium
medium processing flag
radarelab::algo::DBZ dbz
????
bool do_quality
Feature set required for this run.
radarelab::CylindricalVolume cil
Volume resampled as a cylindrical volume.
radarelab::PolarScan< unsigned char > top
Echo top a ???? dBZ [hm].
bool do_beamblocking
beamblocking corretion
void ScrivoStatistica(const radarelab::algo::anaprop::GridStats &)
funzione scrittura matrici statistica
radarelab::PolarScan< unsigned char > beam_blocking
mappa di beam blocking (input)
bool do_zlr_media
Compute ZLR map using averaging.
void caratterizzo_volume()
funzione che caratterizza i volumi polari tramite la qualita'
unsigned MyMAX_BIN
maximum number of beam size
void conversione_convettiva()
Nei punti convettivi ricalcola la Z come MP classica, partendo dalla stima di R (convettiva)
Classe principale del programma.
bool open_from_env(const char *varname, const char *mode, const char *desc=nullptr)
Opens a file taking its name from the environment variable envname.
bool open(const std::string &fname, const char *mode, const char *desc=nullptr)
Opens a file by its pathname.
Open a file taking its name from a given env variable.
void resize_beams_and_propagate_last_bin(unsigned new_beam_size)
Enlarges the PolarScan increasing beam_size and propagating the last bin value.
PolarScan - structure to describe a polarScan containing a matrix of data and conversion factors.
const unsigned max_beam_size() const
Return the maximum beam size in all PolarScans.
const unsigned beam_count
Number of beam_count used ast each elevations.
Homogeneous volume with a common beam count for all PolarScans.
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 ZtoDBZ(double Z)
funzione che converte Z in dBZ
static constexpr double beam_blocking_correction(double val_db, double beamblocking)
@function Compute the corrected value (in dB) given the original value (in dB) and the beam blocking ...
double attenuation(unsigned char DBZbyte, double PIA)
funzione che calcola l'attenuazione totale
static constexpr double DBZtoZ(double DBZ)
funzione che converte dBZ in Z
RadarSite radarSite
RadarSite.
void normalize_elevations(const std::vector< double > &elevations)
Change the elevations in the PolarScans to match the given elevation vector.
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.
std::shared_ptr< LoadInfo > load_info
Polar volume information.
PolarScan< T > & scan(unsigned idx)
Access a polar scan.
Sequence of PolarScans which can have a different beam count for each elevation.
codice principale di elaborazione dei volumi di riflettivita' radar usato per impulso corto
float func_q_Z(unsigned char cl, unsigned char bb, float dst, float dr, float dt, float dh, float dhst, float PIA)
funzione che calcola la qualita' per Z
funzioni che combinano le componenti semplici di qualita' radar
name space generale del programma
Namespace per volume dati.
Codice per il caricamento di volumi ODIM in radarelab.
settaggio ambiente lavoro nel caso non sia settato dall'esterno
definisce struttura Site Contiene le informazioni di base che caratterizzano il sito radar
void classifica_rain()
funzione che classifica la precipitazione se stratiforme o convettiva
void esegui_tutto()
Metodo che lancia tutte le operazioni per il calcolo e la correzione del vpr.
int ier_stampa_vpr
flag d'errore su stampa profilo
radarelab::PolarScan< unsigned char > neve
matrice az-range che memorizza punti di neve
unsigned MyMAX_BIN
LUNGHEZZA MASSIMA.
int heating
variabile di riscaldamento
int ier_max
flag d'errore su calcolo quota max
int corr_vpr()
correzione vpr
int combina_profili(const radarelab::algo::InstantaneousVPR &inst_vpr)
funzione che combina il profilo verticale corrente con quello precedente tramite il metodo di Germann
long int gap
distanza temporale dall'ultimo file vpr [numero acquisizioni intercorse dall'ultimo vpr ?...
void merge_metodi(const radarelab::algo::CalcoloSteiner &steiner, const radarelab::algo::CalcoloVIZ &viz)
fa il merge dei metodi
int stampa_vpr()
stampa profilo combinato
int livmin
quota livello minimo calcolato
radarelab::algo::VPR vpr
Informa se il pixel è convettivo.
int analyse_VPR(float *vpr_liq, int *snow, float *hliq)
funzione che analizza il profilo
int hvprmax
quota picco vpr
double hbbb
altezza bottom brightband
int trovo_hvprmax(int *hmax)
trova il massimo del profilo
CUM_BAC & cum_bac
oggeto CUM_BAC di riferimento
float t_ground
2m temperature
int profile_heating(bool has_inst_vpr)
calcola riscaldamento in quarti d'ora
double htbb
altezza top brightband
CalcoloVPR(CUM_BAC &cum_bac)
Constructor.
radarelab::PolarScan< unsigned char > corr_polar
correzione vpr in byte 0-128 negativa 128-256 positiva, in coord az-ra
Struttara per il calcolo del VPR.
const unsigned max_bin
dimensione matrice
radarelab::Image< unsigned char > cart
vol_pol interpolated in a cartesian map
void write_out(Assets &assets, const std::string tagname, const std::string format="PNG")
Produce in output le immagini PNG dei campi in $DIR_DEBUG.
SingleCart(unsigned max_bin)
Constructor.
void creo_cart(const radarelab::Volume< double > &volume, unsigned int el_index)
conversione da polare a cartesiano alta risoluzione
std::string name
Nome sito radar.
virtual unsigned char get_bin_wind_magic_number(time_t when) const =0
Return the magic number for wind to be used in clean procedure.
RadarSite radarSite
Description of radar site.
virtual std::vector< double > get_elev_array(bool medium=false) const =0
return the elev array used
std::vector< Matrix2D< double > * > slices
Vertical rectangular x,z semi-slices of the cylinder, with one side resting on the cylinder axis.
Radar volume mapped to cylindrical coordinates.
Base for all matrices we use, since we rely on row-major data.
double sample_height(unsigned cell_idx) const
Return the height (in meters) of the sample at the given cell indexa.
std::map< std::string, Scans< double > * > to_load
Map used to specify quantity to be loaded.
void request_quantity(const std::string &name, Scans< double > *volume)
Define a request - Fill to_load attribute
void load(const std::string &pathname)
Load method.
Struttura che eredita da Loader e definisce i metodi per accedere ai dati ODIM.