Elaboradar 0.1
Caricamento in corso...
Ricerca in corso...
Nessun risultato
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
22using namespace radarelab;
23using namespace volume;
24using namespace std;
25
26void 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
34void 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
81void 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
113void 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
140MeltingLayer::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
Homogeneous volume with a common beam count for all PolarScans.
Definition: volume.h:431
unsigned count
height max [km]
Definition: classifier.h:154
double Hmax
height min [km]
Definition: classifier.h:153
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