Elaboradar  0.1
 Tutto Classi Namespace File Funzioni Variabili Tipi enumerati (enum) Gruppi
melting_layer.cpp
1 /*
2  * =====================================================================================
3  *
4  * Filename: melting_layer.cpp
5  *
6  * Description: implementation of melting_layer.h methods
7  *
8  * Version: 1.0
9  * Created: 9/07/2014 15:30
10  * Revision:
11  * Compiler: gcc
12  *
13  * Author: Davide Ori
14  * Organization: Dept.Physics University of Bologna
15  *
16  * =====================================================================================
17  */
18 
19 #include "classifier.h"
21 
22 using namespace radarelab;
23 using namespace volume;
24 using namespace std;
25 
26 void increment(MLpoints& matrix,PolarScan<double>& scan, unsigned az_idx, unsigned rg_idx)
27 {
28  unsigned m_h_idx=matrix.h_idx(scan.height(rg_idx));
29  unsigned m_az_idx=matrix.deg2idx((double)az_idx*360./scan.beam_count);
30  matrix(m_h_idx,m_az_idx)++;
31  matrix.count++;
32 }
33 
34 void MLpoints::box_top_bottom(double box_width_deg, double bot_th, double top_th, std::vector<double>& ML_b, std::vector<double>& ML_t)
35 {
36  if(bot_th<0.||bot_th>1.) cout<<"ERROR bot_th must be 0<%<1 "<<endl;
37  if(top_th<0.||top_th>1.) cout<<"ERROR top_th must be 0<%<1 "<<endl;
38  if(top_th<bot_th) cout<<"ERROR top_th must be > than bot_th"<<endl;
39  int width=1+2*std::floor(0.5*box_width_deg*this->cols()/360.);
40  int half=0.5*(width-1);
41  unsigned box_count=0;
42  double bottom_lim;
43  double top_lim;
44  for(unsigned az=0;az<this->cols();az++)
45  {
46  ML_b[az]= -99.;
47  ML_t[az]= -99.;
48  box_count=0;
49  unsigned round_bm;
50  width=az+half+1;
51  //cout<<az<<"\t";
52  for(int bm=az-half;bm<width;bm++)
53  {
54  if(bm<0) round_bm=this->cols()+bm;
55  else if(bm>=this->cols()) round_bm=bm-this->cols();
56  else round_bm=bm;
57  for(unsigned h=0;h<this->rows();h++)box_count+=(*this)(h,round_bm);
58  }
59  bottom_lim=bot_th*box_count;
60  top_lim=top_th*box_count;
61  //cout<<box_count<<endl;
62  if(box_count>=100)
63  {
64  box_count=0;
65  for(unsigned h=0;h<this->rows();h++)
66  {
67  for(int bm=az-half;bm<width;bm++)
68  {
69  if(bm<0) round_bm=this->cols()+bm;
70  else if(bm>=this->cols()) round_bm=bm-this->cols();
71  else round_bm=bm;
72  box_count+=(*this)(h,round_bm);
73  }
74  if(ML_b[az]<0 && box_count>bottom_lim)ML_b[az]=this->Hmin+h*(this->Hmax-this->Hmin)/this->rows();
75  if(ML_t[az]<0 && box_count>top_lim) ML_t[az]=this->Hmin+h*(this->Hmax-this->Hmin)/this->rows();
76  }
77  }
78  }
79 }
80 
81 void MeltingLayer::seek4mlfile(time_t now, MLpoints& mlp)
82 {
83  unsigned max_time=16*60;
84  ostringstream filename;
85  fstream file;
86  unsigned m_h_idx;
87  unsigned az, app;
88  double h;
89  while(max_time)
90  {
91  max_time-=60;
92  now-=60;
93  filename<<now<<".ml";
94  file.open(filename.str());
95  if(file.is_open())
96  {
97  cout<<endl<<"Apro MLfile "<<filename.str()<<endl;
98  while(file)
99  {
100  file>>app;
101  if(file.eof()) break;
102  file>>az>>h;
103  m_h_idx=mlp.h_idx(h);
104  mlp(m_h_idx,az)++;
105  mlp.count++;
106  }
107  file.close();
108  }
109  filename.seekp(0);
110  }
111 }
112 
113 void MeltingLayer::fill_empty_azimuths()
114 {
115  // X FACILE riempo con la media
116  // DIFFICILE interpolo fra settori buoni
117  // TODO se tutti i punti sono -99 la media è 0 di default!! verifica a posteriori la numerosità del campione della media
118  Statistic<double> mean_bot;
119  Statistic<double> mean_top;
120  for(unsigned i=0;i<bot.size();i++)
121  {
122  if(bot[i]!=-99) mean_bot.feed(bot[i]);
123  if(top[i]!=-99) mean_top.feed(top[i]);
124  }
125  mean_bot.compute_mean();
126  mean_top.compute_mean();
127 
128  if(mean_bot.N==0)
129  {
130  mean_bot.mean=3.0; //TODO: i valori di default dovrebbero essere noti da metodi scientifici (modello, T al suolo, ecc...)
131  mean_top.mean=3.5;
132  }
133  for(unsigned i=0;i<bot.size();i++)
134  {
135  if(bot[i]==-99) bot[i]=mean_bot.mean;
136  if(top[i]==-99) top[i]=mean_top.mean;
137  }
138 }
139 
140 MeltingLayer::MeltingLayer(Volume<double>& vol_z,Volume<double>& vol_zdr,Volume<double>& vol_rhohv,
141  vector< vector< vector< HCA_Park> > >& HCA)
142  : vol_z_0_5km(NUM_AZ_X_PPI), vol_zdr_1km(NUM_AZ_X_PPI), vol_rhohv_1km(NUM_AZ_X_PPI)
143 {
144  cout<<"\tInizio melting Layer"<<endl;
145 // filter(vol_z,vol_z_0_5km,1000.,0.,false); // TODO: se tengo questo range di filtro, semplificare la struttura e riusare i vol_1km vol_2km già filtrati
146 // filter(vol_zdr,vol_zdr_1km,2000.,0.,false);
147 // filter(vol_rhohv,vol_rhohv_1km,2000.,0.,false);
148 
149  Volume<double> dummy(NUM_AZ_X_PPI);
150  filter(vol_z,vol_z_0_5km,500.,3.,false);
151  filter(vol_zdr,dummy,1000.,3.,false);
152  filter(dummy,vol_zdr_1km,1000.,3.,false);
153  filter(vol_rhohv,dummy,1000.,3.,false);
154  filter(dummy,vol_rhohv_1km,1000.,3.,false);
155 
156  cout<<"filtrati"<<endl;
157  //correzione attenuazione con phidp fatta a priori da chi ha invocato
158  //altro preprocessing Ryzhkov 2005b ??? sull'articolo non c'è nulla
159 
160  double MAX_ML_H=4.5; //TODO: check for climatological boundaries in ML height
161  double MIN_ML_H=0.;
162 
163  MLpoints melting_points(MIN_ML_H,MAX_ML_H,vol_z.beam_count,100);
164  top.resize(vol_z.beam_count);
165  bot.resize(vol_z.beam_count);
166  unsigned curr_rg=0;
167  bool confirmed=false;
168 
169  ostringstream ML;
170  ostringstream DATE;
171  DATE<<vol_z.load_info->acq_date<<".ml";
172  ofstream OUT;
173  OUT.open(DATE.str());
174 
175  for(unsigned el=0;el<vol_rhohv_1km.size();el++)
176  {
177  cout<<"\t\t ML el ";
178  PolarScan<double>& rho=vol_rhohv_1km.scan(el);
179  if(rho.elevation>3.&&rho.elevation<10.) //TODO abbasso la soglia minima di un'elevazione
180  {
181  cout<<el<<endl;
182  PolarScan<double>& z=vol_z_0_5km.scan(el);
183  PolarScan<double>& zdr=vol_zdr_1km.scan(el);
184  for(unsigned az=0;az<rho.beam_count;az++)
185  {
186  Statistic<double> average;
187  double min_zdr=0.8;
188  for(unsigned rg=0;rg<rho.beam_size;rg++)
189  {
190  if(zdr(az,rg)>0.2)average.feed(zdr(az,rg));
191  }
192  if(average.N) min_zdr=average.compute_mean();
193  else min_zdr=0.8;
194 
195  for(unsigned rg=0;rg<rho.beam_size;rg++)
196  {
197  //if(el==5)cout<<rg<<" "<<rho.beam_size<<"\t"<<az<<" "<<rho.beam_count<<endl;
198  if(rho(az,rg)>=0.9 && rho(az,rg)<=0.95 && HCA[el][az][rg].meteo_echo() && z.height(rg)>MIN_ML_H && z.height(rg)<MAX_ML_H) //TODO diminuisco la soglia minima di rho da 0.9 a 0.85 e la massima da 0.97 a 0.95
199  {
200  curr_rg=rg;
201  while(curr_rg<z.beam_size && z.diff_height(rg,curr_rg)<0.5 && !confirmed)
202  {
203  //if(el==5&&az==165&&rg==448)cout<<curr_rg<<endl;
204  //if(el==4)if(az==85)if(rg>200)if(rg<250)cout<<rg<<" "<<rho(az,rg)<<" "<<z(az,rg)<<" "<<zdr(az,rg)<<endl;
205  if(z(az,curr_rg)>20. && z(az,curr_rg)<47 && zdr(az,curr_rg)>min_zdr && //TODO cambio la soglia minima da 0.8 a metodo wise
206  zdr(az,curr_rg)<2.5 && z.height(rg)>melting_points.Hmin) // TODO diminuisco soglia minima z da 30 a 20
207  {
208  confirmed=true;
209  }
210  curr_rg++;
211  }
212  if(confirmed)
213  {
214  increment(melting_points,z,az,rg);
215  confirmed=false;
216  ML<<el<<"\t"<<az<<"\t"<<z.height(rg)<<endl;
217  }
218  }
219  }
220  }
221  }
222  }
223  OUT<<ML.str();
224  OUT.close();
225  cout<<endl<<"I punti ML trovati sono "<<melting_points.count<<endl;
226  cout<<"Cerco altri file di ML"<<endl;
227  seek4mlfile(vol_z.load_info->acq_date, melting_points);
228  cout<<"Ora i punti ML sono "<<melting_points.count<<endl;
229  melting_points.box_top_bottom(20.,0.2,0.8,bot,top);
230  fill_empty_azimuths();
231 
232 // cout<<"Altezza ML"<<endl;
233 // for(unsigned i=0;i<bot.size();i++)cout<<bot[i]<<"\t"<<top[i]<<endl;
234 }
Generic Class to perform statistical analysis Statistic object could be used as accumulator of data a...
Definition: statistics.h:21
PolarScan - structure to describe a polarScan containing a matrix of data and conversion factors...
Definition: volume.h:112
unsigned count
height max [km]
Definition: classifier.h:154
unsigned beam_count
Count of beams in this scan.
Definition: volume.h:30
Classes to compute hydrometeor classification.
T compute_mean(unsigned minimum=1)
Compute mean of the distribution of x values.
Definition: statistics.h:177
const unsigned beam_count
Number of beam_count used ast each elevations.
Definition: volume.h:419
std::shared_ptr< LoadInfo > load_info
Polar volume information.
Definition: volume.h:270
MLpoints Melting Layer Points matrix AzH.
Definition: classifier.h:149
unsigned beam_size
Number of samples in each beam.
Definition: volume.h:33
double diff_height(unsigned rg_start, unsigned rg_end)
Height difference in kilometers (legacy) between two range gates.
Definition: volume.cpp:32
double height(unsigned rg, double beam_half_width=0.0)
Height in kilometers (legacy) at range gate for beam elevation + beam_half_width. ...
Definition: volume.cpp:27
double elevation
Nominal elevation of this PolarScan, which may be different from the effective elevation of each sing...
Definition: volume.h:42