11 #include <radarelab/algo/vpr.h>
16 #include "cartproducts.h"
21 #include <radarlib/radar.hpp>
49 #define OVERBLOCKING 51
50 #define SOGLIA_TOP 20.0 // soglia per trovare top
54 #define AMPLITUDE 0.9 // ampiezza fascio radar
57 #define LIMITE_ANAP 240
59 #define DTOR M_PI/180. //fattore conversione gradi-radianti
60 #define CONV_RAD 360./4096.*DTOR // fattore conversione unità angolare radar-radianti
63 using namespace radarelab;
64 using namespace radarelab::algo;
66 namespace elaboradar {
74 struct HRay :
public Matrix2D<double>
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;
119 CUM_BAC::CUM_BAC(
Volume<double>& volume,
const Config& cfg,
const Site& site,
bool medium,
unsigned max_bin)
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())
130 assets.configure(site, volume.
load_info->acq_date);
133 time_t Time = volume.
load_info->acq_date;
134 struct tm* tempo = gmtime(&Time);
136 strftime(
date, 20,
"%Y%m%d%H%M", tempo);
139 algo::compute_top(volume, SOGLIA_TOP,
top);
154 using namespace radarelab::volume;
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);
170 z_volume.normalize_elevations(elev_array);
171 w_volume.normalize_elevations(elev_array);
172 v_volume.normalize_elevations(elev_array);
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);
192 using namespace radarelab::volume;
193 LOG_CATEGORY(
"radar.io");
194 namespace odim = OdimH5v21;
195 LOG_INFO(
"Reading %s for site %s", nome_file, site.
name.c_str());
204 loader.request_quantity(odim::PRODUCT_QUANTITY_DBZH, &dbzh_volume);
205 loader.request_quantity(odim::PRODUCT_QUANTITY_TH, &th_volume);
209 loader.request_quantity(odim::PRODUCT_QUANTITY_VRAD, &v_volume);
210 loader.request_quantity(odim::PRODUCT_QUANTITY_WRAD, &w_volume);
211 loader.request_quantity(odim::PRODUCT_QUANTITY_ZDR, &zdr_volume);
213 loader.load(nome_file);
216 if (dbzh_volume.empty() && th_volume.empty())
218 LOG_ERROR(
"neither DBZH nor TH were found in %s", nome_file);
219 throw runtime_error(
"neither DBZH nor TH were found");
224 for (
auto i: loader.to_load){
225 if(!i.second->empty() ) i.second->normalize_elevations(elev_array);
228 if (!dbzh_volume.empty()) {
229 LOG_WARN(
" DBZH found");
230 z_volume = &dbzh_volume;
233 LOG_WARN(
"no DBZH found: using TH");
234 z_volume = &th_volume;
237 if (do_clean && !w_volume.empty() && !v_volume.empty())
239 if (zdr_volume.empty())
243 for (
unsigned i = 0; i < z_volume->size(); ++i){
245 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);
246 radarelab::algo::Cleaner::evaluateCleanID(z_volume->at(i), w_volume.at(i), v_volume.at(i),full_volume_cleanID.at(i),i);
247 for (
unsigned ibeam=0;ibeam<z_volume->at(i).beam_count; ibeam++)
248 for (
unsigned j=0; j<z_volume->at(i).beam_size; j++){
249 if (full_volume_cleanID.at(i)(ibeam,j) != 0) z_volume->at(i)(ibeam,j)=z_volume->at(i).undetect;
253 for (
unsigned i = 0; i < z_volume->size(); ++i){
254 algo::Cleaner::clean(z_volume->at(i), w_volume.at(i), v_volume.at(i),zdr_volume.at(i),i);
255 algo::Cleaner::clean(z_volume->at(i), w_volume.at(i), v_volume.at(i),zdr_volume.at(i),i+100);
260 algo::azimuthresample::MaxOfClosest<double> resampler;
261 resampler.resample_volume(*z_volume, volume, 1.0);
318 assets.load_dem(
dem);
326 for(
unsigned i=0; i<NUM_AZ_X_PPI; i++)
327 for(
unsigned k=0; k<
volume[0].beam_size; k++)
332 for(
unsigned l=0; l<=el_inf; l++)
335 if (k >=
volume[l].beam_size)
continue;
336 if (k <
volume[el_inf].beam_size)
339 volume[l].set(i, k, MISSING_DB);
347 LOG_INFO(
"declutter_anaprop completed with static declutter");
365 std::vector <double> above_0 (4,0);
366 std::vector <double> above_15 (4,0);
367 std::vector <double> above_30 (4,0);
368 std::vector <double> above_40 (4,0);
369 for(
unsigned el =0; el <4; el++)
371 for (
unsigned k=80; k <
volume[el].beam_size; k ++)
372 if (
volume[el].
get(iaz,k) > 40.){
377 }
else if (
volume[el].
get(iaz,k) > 30.){
381 }
else if (
volume[el].
get(iaz,k) > 15.){
384 }
else if (
volume[el].
get(iaz,k) > 0.){
391 if ( above_15[2]/above_15[0] >= 0.025){
392 if (above_0[1]/above_0[0] >= 0.6 && above_30[2]/above_15[2] <0.15 && above_0[1] >=50000){
393 anaprop.conf_texture_threshold = 5.;
394 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] );
399 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] );
403 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] );
406 LOG_INFO(
"declutter_anaprop completed with anaprop");
411 LOG_WARN(
"declutter_anaprop completed without doing anything");
421 LOG_INFO(
"elabora_Dato completata");
436 LOG_INFO(
"Letta mappa statica");
482 for (
unsigned i=NUM_AZ_X_PPI; i<800; i++)
484 for (
unsigned k=i-1; k<i+2; k++)
485 if(
first_level(i%NUM_AZ_X_PPI, j) < first_level_tmp(k%NUM_AZ_X_PPI, j))
486 first_level(i%NUM_AZ_X_PPI, j)=first_level_tmp(k%NUM_AZ_X_PPI, j);
487 LOG_INFO(
" fine patch %u ",MyMAX_BIN);
494 static const int DIM1_ST = 16;
495 static const int DIM2_ST = 13;
498 static const int N_MIN_BIN = 500;
501 unsigned char statistica[DIM1_ST][DIM2_ST];
502 unsigned char statistica_bl[DIM1_ST][DIM2_ST];
503 unsigned char statistica_el[DIM1_ST][DIM2_ST];
505 memset(statistica,255,DIM1_ST*DIM2_ST);
506 memset(statistica_bl,255,DIM1_ST*DIM2_ST);
507 memset(statistica_el,255,DIM1_ST*DIM2_ST);
509 for(az=0; az<DIM1_ST; az++)
510 for(ran=0; ran<DIM2_ST; ran++){
511 if (grid_stats.count(az, ran) >= N_MIN_BIN)
513 statistica[az][ran] = grid_stats.perc_anap(az, ran);
514 statistica_bl[az][ran] = grid_stats.perc_bloc(az, ran);
515 statistica_el[az][ran] = grid_stats.perc_elev(az, ran);
522 fwrite(
date,12,1,f_stat);
523 fwrite(statistica,DIM1_ST*DIM2_ST,1,f_stat);
528 fwrite(
date,12,1,f_stat);
529 fwrite(statistica_bl,DIM1_ST*DIM2_ST,1,f_stat);
534 fwrite(
date,12,1,f_stat);
535 fwrite(statistica_el,DIM1_ST*DIM2_ST,1,f_stat);
559 LOG_DEBUG(
"start caratterizzo_volume");
562 hray_inf.load_hray_inf(assets);
582 for (
unsigned l=0; l<
volume.size(); l++)
584 const auto& scan =
volume[l];
585 for (
int i=0; i<NUM_AZ_X_PPI; i++)
587 const double elevaz = scan.elevations_real(i);
592 for (
unsigned k=0; k<scan.beam_size; k++)
594 double sample = scan.get(i, k);
597 dist = k * scan.cell_size + scan.cell_size / 2.;
613 dhst = PolarScanBase::sample_height(elevaz + 0.45, dist)
614 - PolarScanBase::sample_height(elevaz - 0.45, dist);
618 if (l<
volume.size() -1 ) {
620 dh = hray_inf(l + 1, k) - hray_inf(l, k);
627 if (l <
anaprop.elev_fin(i, k)) {
630 }
else if (l ==
anaprop.elev_fin(i, k)) {
631 cl=
anaprop.dato_corrotto(i, k);
633 }
else if (l >
anaprop.elev_fin(i, k)) {
641 if (l <
anaprop.elev_fin(i, k)) {
649 qual.
scan(l).set(i, k, (
unsigned char)(
func_q_Z(cl,bb,dist,drrs,hray_inf.dtrs,dh,dhst,PIA)*100));
656 if(cl==0 && bb<BBMAX_VPR )
657 flag_vpr.
scan(l).set(i, k, 1);
663 LOG_DEBUG(
"End caratterizzo_volume");
669 LOG_CATEGORY(
"radar.class");
678 LOG_INFO(
"data= %s",cum_bac.date);
680 gap = cum_bac.assets.read_profile_gap();
683 hmax = cum_bac.assets.read_vpr_hmax();
689 hbbb=(hmax-600.)/1000.;
690 htbb=(hmax+600.)/1000.;
694 LOG_INFO(
"leggo 0termico per class da file %s",getenv(
"FILE_ZERO_TERMICO"));
697 if (cum_bac.assets.read_0term(zeroterm))
700 htbb=zeroterm/1000. + 0.4;
701 hbbb=zeroterm/1000. - 1.0;
703 LOG_ERROR(
"non ho trovato il file dello zero termico");
704 LOG_INFO(
"attenzione, non ho trovat zero termico ne da vpr ne da radiosondaggio");
711 if (hbbb<0.) hbbb=0.;
713 LOG_INFO(
"calcolati livelli sopra e sotto bright band hbbb=%f htbb=%f",hbbb,htbb);
718 LOG_DEBUG (
"Matrice cilindrica Naz %3d Nrange %4d Nheight %4d", cil.
slices.size(), cil.x_size, cil.z_size);
721 algo::CalcoloVIZ viz(cil, htbb, hbbb, t_ground);
722 viz.classifico_VIZ();
728 algo::CalcoloSteiner steiner(cum_bac.volume, cum_bac.anaprop.elev_fin, cil.x_size);
729 steiner.calcolo_background();
730 steiner.classifico_STEINER();
732 merge_metodi(steiner, viz);
738 for (
unsigned j = 0; j < NUM_AZ_X_PPI; ++j)
739 for (
unsigned k = 0; k < steiner.max_bin; ++k)
740 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)
741 conv(j,k) = steiner.conv_STEINER(j, k);
743 if (steiner.conv_STEINER(j, k) == viz.conv_VIZ(j, k) && steiner.conv_STEINER(j, k) > 0 && viz.stratiform(j, k) < 1)
744 conv(j,k) = viz.conv_VIZ(j, k);
768 LOG_CATEGORY(
"radar.vpr");
770 LOG_DEBUG (
" modalita %d", MOD_VPR);
776 combinante = cum_bac.assets.find_vpr0(cum_bac.dbz, vpr0, gap);
778 for (
unsigned i=0; i<vpr0.size(); i++)
779 LOG_DEBUG (
" Profilo vecchio - livello %2d valore %6.2f",i,vpr0.val[i]);
781 LOG_INFO(
"gap %li",gap);
789 if (inst_vpr.success)
791 vpr = combine_profiles(vpr0, inst_vpr.vpr, inst_vpr.cv, inst_vpr.ct);
797 if (inst_vpr.success)
809 LOG_INFO(
" livmin %i", livmin.livmin);
811 if (livmin.idx >= vpr.size() - 1 || !livmin.found)
814 this->livmin = livmin.livmin;
818 cum_bac.assets.write_vpr0(vpr);
819 for (
unsigned i=0; i<vpr.size(); i++) LOG_DEBUG (
" Profilo nuovo - livello %2d valore %6.2f",i,vpr.val[i]);
825 #include <radarelab/vpr_par.h>
827 LOG_CATEGORY(
"radar.vpr");
829 int heating = cum_bac.assets.read_vpr_heating();
839 heating = heating - (gap-1) + 1;
840 if (heating>=WARM) heating=MEMORY;
842 if (heating<0) heating=0;
845 cum_bac.assets.write_vpr_heating(heating);
848 LOG_INFO(
"gap %li heating %i",gap,heating);
859 file.
open_from_env(
"VPR_ARCH",
"wt",
"ultimo vpr in dBZ per il plot");
861 fprintf(file,
" QUOTA DBZ AREA PRECI(KM^2/1000)\n" );
862 for (
unsigned ilay=0; ilay<NMAXLAYER; ilay++){
863 if (vpr.val[ilay]> 0.001 ) {
864 vpr_dbz=cum_bac.dbz.RtoDBZ(vpr.val[ilay]);
865 fprintf(file,
" %i %10.3f %li\n", ilay*TCK_VPR+TCK_VPR/2, vpr_dbz, vpr.area[ilay]);
868 fprintf(file,
" %i %10.3f %li\n", ilay*TCK_VPR+TCK_VPR/2, NODATAVPR, vpr.area[ilay]);
897 #include <radarelab/vpr_par.h>
900 LOG_CATEGORY(
"radar.vpr");
902 int ilray,ilref,ilay2,ier_ana,snow,strat;
903 float corr,vpr_liq,vpr_hray,hbin,hliq;
914 ier_max=trovo_hvprmax(&hvprmax);
915 ier_ana=analyse_VPR(&vpr_liq,&snow,&hliq);
916 LOG_INFO(
"ier_analisi %i",ier_ana) ;
919 if (ier_ana)
return 1;
921 LOG_INFO(
"altezza bright band %i",hvprmax);
922 LOG_INFO(
"CORREGGO VPR");
925 for (
unsigned i=0; i<NUM_AZ_X_PPI; i++){
926 for (
unsigned k=0; k<cum_bac.volume[0].beam_size; k++){
930 hbin=(float)cum_bac.anaprop.quota(i, k);
933 if (snow) neve(i, k)=1;
935 if (cum_bac.do_class)
937 if (conv(i,k) >= CONV_VAL){
942 if (cum_bac.volume[0].get(i, k) > THR_CORR && hbin > hliq && strat){
945 ilray=(hbin>=livmin)?(floor(hbin/TCK_VPR)):(floor(livmin/TCK_VPR));
946 if (ilray>= NMAXLAYER ) ilray=NMAXLAYER-2;
949 if ((
int)hbin%TCK_VPR > TCK_VPR/2) ilay2=ilray+1;
951 if (ilay2< floor(livmin/TCK_VPR)) ilay2=floor(livmin/TCK_VPR);
956 ilref=(cum_bac.dem(i, k)>livmin)?(floor(cum_bac.dem(i, k)/TCK_VPR)):(floor(livmin/TCK_VPR));
959 if (vpr.val[ilref] > 0 && vpr.val[ilray] > 0 ){
961 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);
964 if (cum_bac.dem(i, k)> hvprmax+HALF_BB-TCK_VPR || snow){
984 corr=cum_bac.dbz.RtoDBZ(vpr.val[ilref])-cum_bac.dbz.RtoDBZ(vpr_hray);
986 cum_bac.volume[0].set(i, k, cum_bac.dbz.DBZ_snow(cum_bac.volume[0].get(i, k)));
990 corr = cum_bac.dbz.RtoDBZ_class(vpr_liq) - cum_bac.dbz.RtoDBZ_class(vpr_hray);
993 if (corr>MAX_CORR) corr=MAX_CORR;
994 if (hbin<hvprmax && corr>0.) corr=0;
997 double corrected = cum_bac.volume[0].get(i, k) + corr;
998 if (corrected > MAXVAL_DB)
999 cum_bac.volume[0].set(i, k, MAXVAL_DB);
1000 else if ( corrected < MINVAL_DB)
1001 cum_bac.volume[0].set(i, k, MINVAL_DB);
1003 cum_bac.volume[0].set(i, k, corrected);
1005 corr_polar(i, k)=(
unsigned char)(corr)+128;
1021 int i,imax,istart,foundlivmax;
1022 float h0start,peak,soglia;
1025 if (t_ground != NODATAVPR)
1027 LOG_DEBUG(
"trovo hvprmax a partire da 400 m sotto lo zero dell'adiabatica secca");
1028 h0start=t_ground/9.8*1000 ;
1029 istart=h0start/TCK_VPR -2;
1030 if (istart< livmin/TCK_VPR) istart=livmin/TCK_VPR;
1031 LOG_DEBUG(
"t_ground h0start istart %f %f %i",t_ground,h0start,istart);
1034 LOG_DEBUG(
"trovo hvprmax a partire da livmin");
1035 istart=livmin/TCK_VPR+1;
1047 soglia = DBZ::DBZtoR(THR_VPR,200,1.6);
1050 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);
1051 if (vpr.val[istart] >soglia && vpr.val[istart+4] > soglia){
1052 peak=10*log10(vpr.val[istart]/vpr.val[istart+4]);
1053 LOG_DEBUG(
"peak1 = %f",peak);
1056 if(peak> MIN_PEAK_VPR){
1059 LOG_DEBUG(
"il primo punto soddisfa le condizioni di picco");
1061 for (i=istart+1;i<NMAXLAYER-4;i++)
1063 if (vpr.val[i] <soglia || vpr.val[i+4] < soglia)
break;
1064 peak=10*log10(vpr.val[i]/vpr.val[i+4]);
1065 if (vpr.val[i]>vpr.val[i-1] && peak> MIN_PEAK_VPR )
1070 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);
1073 if ( imax > INODATA ){
1075 peak=10*log10(vpr.val[imax]/vpr.val[imax+4]);
1076 *hmax=imax*TCK_VPR+TCK_VPR/2;
1077 LOG_DEBUG(
"trovato ilaymax %i %i",*hmax,imax);
1078 LOG_DEBUG(
" picco in dbR %f",peak);
1081 LOG_DEBUG(
"exit status trovo_hvprmax %i",foundlivmax);
1082 return (foundlivmax);
1102 int ier=1,ier_ana=0,liv0;
1103 char date[]=
"000000000000";
1111 int tipo_profilo=-1;
1112 float v600sottobb=NODATAVPR;
1113 float v1000=NODATAVPR;
1114 float v1500=NODATAVPR;
1115 float vliq=NODATAVPR;
1116 float vhliquid=NODATAVPR;
1117 float vprmax=NODATAVPR;
1123 if (t_ground == NODATAVPR)
1125 LOG_WARN(
"non ho T,...");
1128 LOG_ERROR(
" non ho trovato hvprmax, nè T, esco");
1136 if (t_ground >= T_MAX_ML+0.65*(
float)(livmin+TCK_VPR/2)/100.){
1138 LOG_ERROR(
" temperatura alta e non ho trovato hvprmax, esco");
1146 if (t_ground >= T_MAX_SN && t_ground < T_MAX_ML+0.65*(
float)(livmin+TCK_VPR/2)/100. )
1150 LOG_INFO(
" temperatura da scioglimento e massimo in quota");
1154 LOG_ERROR(
" temperatura da scioglimento ma superiore a temperatura max neve e non ho trovato hvprmax, esco");
1158 liv0=livmin+HALF_BB;
1159 if (hvprmax > liv0) LOG_INFO(
" il livello %i è sotto la Bright band, ma T bassa interpolo",livmin);
1160 else LOG_INFO(
" il livello %i potrebbe essere dentro la Bright Band, interpolo",livmin);
1165 if (t_ground < T_MAX_SN)
1168 LOG_INFO(
" temperatura da neve o scioglimento e massimo in quota");
1172 LOG_INFO(
" temperatura da neve o scioglimento e massimo non trovato,neve , non interpolo");
1181 InterpolaVPR_GSL iv;
1189 ier=iv.interpola_VPR(vpr.val.data(), hvprmax, livmin);
1191 LOG_INFO(
" interpolazione fallita");
1192 switch (tipo_profilo)
1203 *vpr_liq=vpr.val[(hvprmax+1000)/TCK_VPR]*2.15;
1209 LOG_INFO(
" interpolazione eseguita con successo");
1212 const char* vpr_arch = getenv(
"VPR_ARCH");
1213 if (!vpr_arch)
throw runtime_error(
"VPR_ARCH is not defined");
1214 string fname(vpr_arch);
1217 file.
open(fname,
"wt",
"vpr interpolato");
1218 for (
unsigned i = 0; i < NMAXLAYER; ++i)
1219 fprintf(file,
" %f \n", cum_bac.dbz.RtoDBZ(iv.vpr_int[i]));
1222 if (tipo_profilo == 2 ) {
1223 *hliq=(iv.E-2.1*iv.G)*1000.;
1227 *vpr_liq=vpr.val[(hvprmax+1000)/TCK_VPR]*2.15;
1230 *hliq=(iv.E-2.1*iv.G)*1000.;
1232 if ( *hliq > livmin) {
1233 *vpr_liq=vpr.val[(int)(*hliq/TCK_VPR)];
1237 if (*hliq<0) *hliq=0;
1251 LOG_INFO(
"TIPO_PROFILO= %i vpr_liq %f hliq %f", tipo_profilo, *vpr_liq,*hliq );
1258 Time = cum_bac.volume.load_info->acq_date;
1259 tempo = gmtime(&Time);
1260 sprintf(date,
"%04d%02d%02d%02d%02d",tempo->tm_year+1900, tempo->tm_mon+1,
1261 tempo->tm_mday,tempo->tm_hour, tempo->tm_min);
1263 if(*hliq > livmin +200 )
1264 vhliquid=cum_bac.dbz.RtoDBZ(vpr.val[(
int)(*hliq)/TCK_VPR]);
1265 vliq=cum_bac.dbz.RtoDBZ(*vpr_liq);
1268 if ( hvprmax-600 >= livmin )
1269 v600sottobb=cum_bac.dbz.RtoDBZ(vpr.val[(hvprmax-600)/TCK_VPR]);
1270 if ((hvprmax+1000)/TCK_VPR < NMAXLAYER )
1271 v1000=cum_bac.dbz.RtoDBZ(vpr.val[(hvprmax+1000)/TCK_VPR]);
1272 if ((hvprmax+1500)/TCK_VPR < NMAXLAYER )
1273 v1500=cum_bac.dbz.RtoDBZ(vpr.val[(hvprmax+1500)/TCK_VPR]);
1274 vprmax=cum_bac.dbz.RtoDBZ(vpr.val[(hvprmax/TCK_VPR)]);
1277 if (FILE* test_vpr=fopen(getenv(
"TEST_VPR"),
"a+"))
1279 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);
1287 cum_bac.assets.write_vpr_hmax(hvprmax);
1288 LOG_INFO(
"fatta scrittura hmax vpr = %d",hvprmax);
1295 for (
unsigned i=0; i<NUM_AZ_X_PPI; i++)
1296 for (
unsigned k=0; k<
volume[0].beam_size; k++)
1330 void CUM_BAC::generate_maps(CartProducts& products)
1333 LOG_INFO(
"Scrittura File Precipitazione 1X1\n");
1336 std::function<unsigned char(const vector<double>&)> convert = [
this](
const vector<double>& samples) {
1343 for (
const auto& s: samples)
1344 sum += DBZ::DBZtoZ(s);
1345 unsigned char res = DBZ::DBtoBYTE(DBZ::ZtoDBZ(sum / samples.size()));
1347 if (res == 0)
return (
unsigned char)1;
1350 products.scaled.to_cart_average(
volume[0], convert, products.z_out);
1352 std::function<unsigned char(unsigned, unsigned)> assign_cart =
1353 [
this](
unsigned azimuth,
unsigned range) {
1355 unsigned char sample = DBZ::DBtoBYTE(
volume[0].
get(azimuth, range));
1356 return max(sample, (
unsigned char)1);
1358 products.scaled.to_cart(assign_cart, products.z_out);
1361 std::function<unsigned char(unsigned, unsigned)> assign_cart =
1362 [
this](
unsigned azimuth,
unsigned range) {
1364 unsigned char sample = DBZ::DBtoBYTE(
volume[0].
get(azimuth, range));
1365 return max(sample, (
unsigned char)1);
1367 products.fullres.to_cart(assign_cart, products.z_fr);
1369 products.scaled.to_cart(
top, products.top_1x1);
1372 const auto& elev_fin =
anaprop.elev_fin;
1373 const auto& quota =
anaprop.quota;
1375 std::function<unsigned char(unsigned, unsigned)> assign_qual =
1376 [
this, &elev_fin](
unsigned azimuth,
unsigned range) {
1377 const auto& el = elev_fin(azimuth, range);
1378 if (range >=
volume[el].beam_size)
1379 return (
unsigned char)0;
1380 return qual.
scan(el).get(azimuth, range);
1382 products.scaled.to_cart(assign_qual, products.qual_Z_1x1);
1384 std::function<unsigned char(unsigned, unsigned)> assign_quota =
1385 ["a](
unsigned azimuth,
unsigned range) {
1386 return 128 + round(quota(azimuth, range) / 100.0);
1388 products.scaled.to_cart(assign_quota, products.quota_1x1);
1389 products.scaled.to_cart(
anaprop.dato_corrotto, products.dato_corr_1x1);
1391 std::function<unsigned char(unsigned, unsigned)> assign_elev_fin = [&elev_fin](
unsigned azimuth,
unsigned range) {
1392 return elev_fin(azimuth, range);
1394 products.scaled.to_cart(assign_elev_fin, products.elev_fin_1x1);
1396 products.scaled.to_cart(
beam_blocking, products.beam_blocking_1x1);
1402 std::function<unsigned char(unsigned, unsigned)> assign = [&neve](
unsigned azimuth,
unsigned range) {
1403 return neve(azimuth, range) ? 0 : 1;
1405 products.scaled.to_cart(assign, products.neve_1x1);
1410 products.scaled.to_cart(
calcolo_vpr->conv, products.conv_1x1);
1418 for (
unsigned int i=0; i<
SD_Z6.size(); i++){
1419 SC_SD.creo_cart(
SD_Z6, i);
1420 std::ostringstream oss;
1422 SC_SD.write_out(assets,oss.str());
1432 inst_vpr(cum_bac.
volume, cum_bac.
qual, cum_bac.flag_vpr, cum_bac.
site.vpr_iaz_min, cum_bac.
site.vpr_iaz_max),
1433 conv(NUM_AZ_X_PPI, cum_bac.
volume.max_beam_size(),0),
1434 corr_polar(NUM_AZ_X_PPI, cum_bac.
volume.max_beam_size()),
1435 neve(NUM_AZ_X_PPI, cum_bac.
volume.max_beam_size())
1437 logging_category = log4c_category_get(
"radar.vpr");
1452 CalcoloVPR::~CalcoloVPR()
1459 printf (
"calcolo VPR \n") ;
1468 LOG_INFO(
"fatta func vpr %s", inst_vpr.success ?
"ok" :
"errore");
1469 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]);
1475 LOG_INFO (
"exit status calcolo VPR istantaneo: %s", inst_vpr.success ?
"ok":
"fallito") ;
1476 LOG_INFO(
"exit status combinaprofili: (1--fallito 0--ok) %i ",ier_comb) ;
1482 printf (
"heating %i \n",
heating);
1483 LOG_INFO(
"ier_comb %i", ier_comb);
1485 if (inst_vpr.success)
1490 if (!ier_comb &&
heating >= WARM){
1493 LOG_INFO(
"exit status correggo vpr: (1--fallito 0--ok) %i",ier) ;
1498 if ( ! ier && inst_vpr.success)
1506 Image<double> azimut;
1507 Image<double> range;
1509 CartData(
int max_bin=512)
1510 : azimut(max_bin), range(max_bin)
1512 for(
int i=0; i<max_bin; i++)
1513 for(
int j=0; j<max_bin; j++)
1515 range(i, j) = hypot(i+.5,j+.5) ;
1516 azimut(i, j) = 90. - atan((j+.5)/(i+.5)) * M_1_PI*180.;
1530 LOG_CATEGORY(
"radar.singlecart");
1534 int x,y,iaz,az_min,az_max;
1536 CartData cd(max_bin);
1538 for(
unsigned i=0; i<max_bin *2; i++)
1539 for(
unsigned j=0; j<max_bin *2; j++)
1540 cart(i, j) = MISSING;
1542 LOG_INFO(
"Creo_cart - %u", max_bin);
1544 for(
unsigned quad=0; quad<4; quad++)
1545 for(
unsigned i=0; i<
max_bin; i++)
1546 for(
unsigned j=0; j<
max_bin; j++)
1548 unsigned irange = (unsigned)round(cd.range(i, j));
1549 if (irange >= max_bin)
1556 az = cd.azimut(i, j);
1561 az = cd.azimut(i, j) + 90.;
1566 az = cd.azimut(i, j) + 180.;
1571 az = cd.azimut(i, j)+270.;
1575 az_min = (int)((az - .45)/.9);
1576 az_max = ceil((az + .45)/.9);
1581 az_min = az_min + NUM_AZ_X_PPI;
1582 az_max = az_max + NUM_AZ_X_PPI;
1584 for(iaz = az_min; iaz<az_max; iaz++){
1586 unsigned char sample = 0;
1587 if (irange < volume[el_index].beam_size)
1588 sample = max((
unsigned char) (volume[el_index].
get(iaz%NUM_AZ_X_PPI, irange)), (
unsigned char)1);
1589 if(
cart(y, x) <= sample)
cart(y, x) = sample;
1596 if (getenv(
"DIR_DEBUG") == NULL)
return;
1609 tempo=gmtime(&clock);
1611 return asctime(tempo);
1616 static time_t time_tot = 0,time1 = 0,time2 = 0;
1617 LOG_CATEGORY(
"radar.timing");
1623 time2 = time(&time2);
1625 LOG_INFO(
"tempo parziale %ld ---- totale %ld", time2-time1, time2-time_tot);
radarelab::PolarScan< unsigned char > corr_polar
correzione vpr in byte 0-128 negativa 128-256 positiva, in coord az-ra
radarelab::algo::DBZ dbz
????
bool do_quality
Feature set required for this run.
bool do_readStaticMap
Read Static clutter map.
double htbb
altezza top brightband
PolarScan< T > & scan(unsigned idx)
Access a polar scan.
definisce struttura Site Contiene le informazioni di base che caratterizzano il sito radar ...
radarelab::PolarScan< unsigned char > bb_first_level
mappa di elevazioni da beam blocking (input)
Codice per il caricamento di volumi ODIM in radarelab.
float read_t_ground() const
fornisce temperatura al suolo, da lettura file esterno
Struttura che eredita da Loader e definisce i metodi per accedere ai dati ODIM.
CUM_BAC & cum_bac
oggeto CUM_BAC di riferimento
radarelab::PolarScan< unsigned char > neve
matrice az-range che memorizza punti di neve
unsigned MyMAX_BIN
LUNGHEZZA MASSIMA.
int ier_stampa_vpr
flag d'errore su stampa profilo
bool do_medium
medium processing flag
bool do_beamblocking
beamblocking corretion
settaggio ambiente lavoro nel caso non sia settato dall'esterno
void resize_beams_and_propagate_last_bin(unsigned new_beam_size)
Enlarges the PolarScan increasing beam_size and propagating the last bin value.
int corr_vpr()
correzione vpr
bool do_anaprop
anaprop correction
radarelab::PolarScan< unsigned char > first_level_static
mappa statica
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.
int profile_heating(bool has_inst_vpr)
calcola riscaldamento in quarti d'ora
int combina_profili(const radarelab::algo::InstantaneousVPR &inst_vpr)
funzione che combina il profilo verticale corrente con quello precedente tramite il metodo di Germann...
Struttara per il calcolo del VPR.
CalcoloVPR * calcolo_vpr
Oggetto per calcolare e correggere con VPR.
void leggo_first_level()
funzione che legge la mappa statica e la mappa di elevazioni da beam blocking e le condensa in un uni...
RadarSite radarSite
Description of radar site.
void classifica_rain()
funzione che classifica la precipitazione se stratiforme o convettiva
Base for all matrices we use, since we rely on row-major data.
const unsigned beam_count
Number of beam_count used ast each elevations.
radarelab::PolarScan< unsigned char > first_level
mappa dinamica complessiva
virtual std::vector< double > get_elev_array(bool medium=false) const =0
return the elev array used
bool do_zlr_media
Compute ZLR map using averaging.
unsigned MyMAX_BIN
maximum number of beam size
void ScrivoStatistica(const radarelab::algo::anaprop::GridStats &)
funzione scrittura matrici statistica
std::shared_ptr< LoadInfo > load_info
Polar volume information.
codice principale di elaborazione dei volumi di riflettivita' radar usato per impulso corto ...
int trovo_hvprmax(int *hmax)
trova il massimo del profilo
PolarScan< T > & append_scan(unsigned beam_count, unsigned beam_size, double elevation, double cell_size)
Append a scan to this volume.
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.
SingleCart(unsigned max_bin)
Constructor.
funzioni che combinano le componenti semplici di qualita' radar
CalcoloVPR(CUM_BAC &cum_bac)
Constructor.
float t_ground
2m temperature
void esegui_tutto()
Metodo che lancia tutte le operazioni per il calcolo e la correzione del vpr.
char date[20]
Acquisition date.
log4c_category_t * logging_category
logging category
double hbbb
altezza bottom brightband
radarelab::Volume< double > & volume
Polar volume of Reflectivity.
int hvprmax
quota picco vpr
radarelab::PolarScan< float > dem
dem in coordinate azimut range
int stampa_vpr()
stampa profilo combinato
radarelab::CylindricalVolume cil
Volume resampled as a cylindrical volume.
Classe principale del programma.
Homogeneous volume with a common beam count for all PolarScans.
bool do_class
Convective-stratiform classification.
void write_last_vpr()
Write the acquisition time in $LAST_VPR file.
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.
bool do_bloccorr
bloccorrection
static void read_odim_volume(radarelab::Volume< double > &volume, const Site &site, const char *nome_file, bool do_clean=false, bool do_medium=false)
Read from ODIM data file.
bool do_devel
Produce additional output.
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
void vpr_class()
Esegue tutta la catena vpr (e classificazione) se richiesta.
radarelab::PolarScan< unsigned char > beam_blocking
mappa di beam blocking (input)
RadarSite radarSite
RadarSite.
radarelab::Volume< unsigned char > qual
qualita volume polare
radarelab::Image< unsigned char > cart
vol_pol interpolated in a cartesian map
const unsigned max_beam_size() const
Return the maximum beam size in all PolarScans.
Finds resources, like data files, used by the program.
void creo_cart(const radarelab::Volume< double > &volume, unsigned int el_index)
conversione da polare a cartesiano alta risoluzione
std::vector< Matrix2D< double > * > slices
Vertical rectangular x,z semi-slices of the cylinder, with one side resting on the cylinder axis...
double attenuation(unsigned char DBZbyte, double PIA)
funzione che calcola l'attenuazione totale
void caratterizzo_volume()
funzione che caratterizza i volumi polari tramite la qualita'
const unsigned max_bin
dimensione matrice
int analyse_VPR(float *vpr_liq, int *snow, float *hliq)
funzione che analizza il profilo
int heating
variabile di riscaldamento
Open a file taking its name from a given env variable.
radarelab::PolarScan< unsigned char > top
Echo top a ???? dBZ [hm].
radarelab::Volume< double > SD_Z6
Polar volume of standard deviation of reflectivity over 6 km length.
std::string name
Nome sito radar.
Radar volume mapped to cylindrical coordinates.
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.
void merge_metodi(const radarelab::algo::CalcoloSteiner &steiner, const radarelab::algo::CalcoloVIZ &viz)
fa il merge dei metodi
void want_vpr()
Call this just after creating the CUM_BAC object, to signal that VPR should also be computed...
bool open(const std::string &fname, const char *mode, const char *desc=nullptr)
Opens a file by its pathname.
radarelab::algo::Anaprop< double > anaprop
Oggetto per correzione ANAPRO.
const Site & site
site information object
void declutter_anaprop()
funzione che elabora il dato radar rimuovendo anaprop e beam blocking
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.
bool do_declutter
use only static declutter map
void conversione_convettiva()
Nei punti convettivi ricalcola la Z come MP classica, partendo dalla stima di R (convettiva) ...