Elaboradar  0.1
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  float max_z,max_zdr;
163 
164  MLpoints melting_points(MIN_ML_H,MAX_ML_H,vol_z.beam_count,100);
165  top.resize(vol_z.beam_count);
166  bot.resize(vol_z.beam_count);
167  unsigned curr_rg=0,rg_maxz=0,rg_maxzdr=0;
168  bool confirmed=false;
169 
170  ostringstream ML;
171  ostringstream DATE;
172  DATE<<vol_z.load_info->acq_date<<".ml";
173  ofstream OUT;
174  OUT.open(DATE.str());
175 
176  for(unsigned el=0;el<vol_rhohv_1km.size();el++)
177  {
178  cout<<"\t\t ML el ";
179  PolarScan<double>& rho=vol_rhohv_1km.scan(el);
180  if(rho.elevation>3.&&rho.elevation<10.) //TODO abbasso la soglia minima di un'elevazione
181  {
182  cout<<el<<endl;
183  PolarScan<double>& z=vol_z_0_5km.scan(el);
184  PolarScan<double>& zdr=vol_zdr_1km.scan(el);
185  for(unsigned az=0;az<rho.beam_count;az++)
186  {
187  Statistic<double> average;
188  double min_zdr=0.8;
189  //calcolo la media sul raggio di zdr
190  for(unsigned rg=0;rg<rho.beam_size;rg++)
191  {
192  if(zdr(az,rg)>0.2) average.feed(zdr(az,rg));
193  }
194  if(average.N) min_zdr=average.compute_mean();
195  else min_zdr=0.8;
196  for(unsigned rg=0;rg<rho.beam_size;rg++)
197  {
198  //if(el==5)cout<<rg<<" "<<rho.beam_size<<"\t"<<az<<" "<<rho.beam_count<<endl;
199  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
200  {
201  curr_rg=rg;
202  max_z=-19.;
203  max_zdr=-5;
204  //scorro il raggio entro 0.5 km dal bin alla ricerca di un bin con valori accettabili di z e zdr
205  while (curr_rg<z.beam_size && z.diff_height(rg,curr_rg<0.5)){
206  if (z(az,curr_rg) > max_z)
207  {
208 
209  max_z=z(az,curr_rg);
210  rg_maxz=rg;
211 
212  }
213  if (zdr(az,curr_rg)>max_zdr)
214  {
215  max_zdr=zdr(az,curr_rg);
216  rg_maxzdr=rg;
217  }
218 
219  curr_rg++; }
220  curr_rg=rg;
221  while(curr_rg<z.beam_size && z.diff_height(rg,curr_rg)<0.5 && !confirmed)
222  {
223  //if(el==5&&az==165&&rg==448)cout<<curr_rg<<endl;
224  //if(el==4)if(az==85)if(rg>200)if(rg<250)cout<<rg<<" "<<rho(az,rg)<<" "<<z(az,rg)<<" "<<zdr(az,rg)<<endl;
225  if(z(az,rg_maxz)>20. && max_z <47 && max_zdr >min_zdr && //TODO cambio la soglia minima da 0.8 a metodo wise
226  max_zdr < 2.5 && z.height(rg)>melting_points.Hmin) // TODO diminuisco soglia minima z da 30 a 20
227 
228  {
229  confirmed=true;
230  }
231  curr_rg++;
232  }
233  if(confirmed)
234  {
235  increment(melting_points,z,az,rg);
236  confirmed=false;
237  ML<<el<<"\t"<<az<<"\t"<<z.height(rg)<<endl;
238  }
239  }
240  }
241  }
242  }
243  }
244  OUT<<ML.str();
245  OUT.close();
246  cout<<endl<<"I punti ML trovati sono "<<melting_points.count<<endl;
247  cout<<"Cerco altri file di ML"<<endl;
248  seek4mlfile(vol_z.load_info->acq_date, melting_points);
249  cout<<"Ora i punti ML sono "<<melting_points.count<<endl;
250  melting_points.box_top_bottom(20.,0.2,0.8,bot,top);
251  fill_empty_azimuths();
252 
253 // cout<<"Altezza ML"<<endl;
254 // for(unsigned i=0;i<bot.size();i++)cout<<bot[i]<<"\t"<<top[i]<<endl;
255 }
Classes to compute hydrometeor classification.
PolarScan - structure to describe a polarScan containing a matrix of data and conversion factors.
Definition: volume.h:113
T compute_mean(unsigned minimum=1)
Compute mean of the distribution of x values.
Definition: statistics.h:177
Generic Class to perform statistical analysis Statistic object could be used as accumulator of data a...
Definition: statistics.h:22
const unsigned beam_count
Number of beam_count used ast each elevations.
Definition: volume.h:434
unsigned count
height max [km]
Definition: classifier.h:154
MLpoints Melting Layer Points matrix AzH.
Definition: classifier.h:150
std::shared_ptr< LoadInfo > load_info
Polar volume information.
Definition: volume.h:272
String functions.
Definition: cart.cpp:4
double diff_height(unsigned rg_start, unsigned rg_end)
Height difference in kilometers (legacy) between two range gates.
Definition: volume.cpp:32
unsigned beam_size
Number of samples in each beam.
Definition: volume.h:33
double elevation
Nominal elevation of this PolarScan, which may be different from the effective elevation of each sing...
Definition: volume.h:42
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
unsigned beam_count
Count of beams in this scan.
Definition: volume.h:30