Daily_QC
estimate_missing_stations.c
Go to the documentation of this file.
1 
2 #include "prototypes_new.h"
3 
4 void estimate_missing_stations(int j)
5 
6 {
7 
8 extern int isom;
9 extern int isohyets_used;
10 extern struct pdata pdata[10];
11 extern struct station station[3000];
12 extern int max_stations;
13 extern int method;
14 int m,i,h,l,ii;
15 double lat1,lon1,fdist,fdata,fvalue[4],lat,lon,testdist,isoh,isoh1,padj,distlon;
16 float conv=.0174;
17 double fvalue24,fvalue06,fmult,stotal;
18 int num_missing;
19 
20 for(m=0;m<max_stations;m++) {
21 
22  /* only estimate missing data */
23 
24  if(pdata[j].stn[m].frain[4].data >= 0 &&
25  pdata[j].stn[m].frain[4].qual != 5 &&
26  pdata[j].stn[m].frain[4].qual != 1)
27  continue;
28 
29  for(h=4;h>=0;h--) {
30 
31  if(pdata[j].used[h]==0) {
32 
33  fvalue24=-1;
34  continue;
35 
36  }
37 
38 
39  lat1=station[m].lat;
40  lon1=station[m].lon;
41 
42  /*get isohyet */
43 
44  if(isohyets_used!=0)
45  isoh1=station[m].isoh[isom];
46 
47  fdist=0.0;
48  fdata=0.0;
49  l=0;
50 
51  for(ii=0;ii<30;ii++){
52 
53  i=station[m].index[ii];
54 
55  /* dont estimate unless good or forced good */
56 
57  if(pdata[j].stn[i].frain[h].qual!=0 &&
58  pdata[j].stn[i].frain[h].qual!=8 &&
59  pdata[j].stn[i].frain[h].qual!=6 &&
60  pdata[j].stn[i].frain[h].qual!=3 &&
61  pdata[j].stn[i].frain[h].qual!=2)
62  continue;
63 
64  /*dont use missing stations */
65 
66  if(pdata[j].stn[i].frain[h].data < 0)
67  continue;
68 
69  lat=station[i].lat;
70  lon=station[i].lon;
71 
72  if(isohyets_used != 0)
73  isoh=station[i].isoh[isom];
74 
75  distlon=(lon1-lon)*cos((lat1+lat)/2*conv);
76 
77  testdist= pow((double)(lat1-lat),2) +
78  pow((double)distlon,2);
79 
80  if(testdist==0.0)
81  testdist=.0001;
82 
83  if(method==2 && isoh > 0 && isoh1 > 0)
84  padj=pdata[j].stn[i].frain[h].data*isoh1/isoh;
85 
86  else
87  padj=pdata[j].stn[i].frain[h].data;
88 
89 
90  fdist=1/testdist+fdist;
91  fdata=padj/testdist + fdata;
92 
93  l++;
94 
95  }
96 
97 
98  if(l < 5) {
99 
100  fdist=0.0;
101  fdata=0.0;
102  l=0;
103 
104  for(i=0;i<max_stations;i++){
105 
106  if(i==m)
107  continue;
108 
109  /* dont estimate unless good or forced good */
110 
111  if(pdata[j].stn[i].frain[h].qual!=0 &&
112  pdata[j].stn[i].frain[h].qual!=8 &&
113  pdata[j].stn[i].frain[h].qual!=6 &&
114  pdata[j].stn[i].frain[h].qual!=3 &&
115  pdata[j].stn[i].frain[h].qual!=2)
116  continue;
117 
118  /*dont use missing stations */
119 
120  if(pdata[j].stn[i].frain[h].data < 0)
121  continue;
122 
123  lat=station[i].lat;
124  lon=station[i].lon;
125 
126  if(isohyets_used != 0)
127  isoh=station[i].isoh[isom];
128 
129  distlon=(lon1-lon)*cos((lat1+lat)/2*conv);
130 
131  testdist= pow((double)(lat1-lat),2) +
132  pow((double)distlon,2);
133 
134  if(testdist==0.0)
135  testdist=.0001;
136 
137  if(method==2 && isoh > 0 && isoh1 > 0)
138  padj=pdata[j].stn[i].frain[h].data*isoh1/isoh;
139 
140  else
141  padj=pdata[j].stn[i].frain[h].data;
142 
143 
144  fdist=1/testdist+fdist;
145  fdata=padj/testdist + fdata;
146 
147  l++;
148 
149  }
150 
151  }
152 
153  /* 24 hourly has stations */
154 
155  if(h==4 && l > 0)
156  fvalue24=fdata/fdist;
157 
158  /* 24 hourly has no stations */
159 
160  else if(h==4)
161  fvalue24=-1;
162 
163  /* 6 hourly has stations */
164 
165  else if(l > 0)
166  fvalue[h]=fdata/fdist;
167 
168  /* 6 hourly no stations */
169 
170  else
171  fvalue[h]=-1;
172 
173  }
174 
175 
176 /* have all 24 hour data */
177 
178 /* 24 hourly data available */
179 
180 if(fvalue24 >= 0) {
181 
182  fvalue06=0.0;
183  stotal=0.0;
184  num_missing=0;
185 
186  for(h=0;h<4;h++) {
187 
188  /* use good, forced good and questionable data to estimate */
189 
190  /* need to caculate partial total to ensure that 24 hourly
191  and 6 hourly data match */
192 
193  if((pdata[j].stn[m].frain[h].qual == 0 ||
194  pdata[j].stn[m].frain[h].qual == 8 ||
195  pdata[j].stn[m].frain[h].qual == 3 ||
196  pdata[j].stn[m].frain[h].qual == 2) &&
197  pdata[j].stn[m].frain[h].data >=0)
198  stotal=stotal+pdata[j].stn[m].frain[h].data;
199 
200  else {
201 
202  num_missing++;
203  fvalue06=fvalue06+fvalue[h];
204 
205  }
206 
207 
208  }
209 
210  /* stotal will now be difference between 24 hour estimate and 6 hourly
211  partial total */
212 
213  stotal=fvalue24-stotal;
214 
215  if(stotal <= 0)
216  stotal=0;
217 
218  if(fvalue06==0)
219  fmult=0;
220 
221  /* fmult is that ratio of actual summed (missing) 6 hourly rain to estimated
222  6 hourly rain */
223 
224  else
225  fmult=stotal/fvalue06;
226 
227  /* now rescale the estimates so they equal the 24 hour estimate */
228 
229  for(h=0;h<4;h++) {
230 
231  if((pdata[j].stn[m].frain[h].qual != 0 &&
232  pdata[j].stn[m].frain[h].qual != 8 &&
233  pdata[j].stn[m].frain[h].qual != 3 &&
234  pdata[j].stn[m].frain[h].qual != 2) ||
235  pdata[j].stn[m].frain[h].data < 0) {
236 
237  if(fvalue06 != 0)
238  pdata[j].stn[m].frain[h].data=fvalue[h]*fmult;
239 
240  else
241  pdata[j].stn[m].frain[h].data=stotal/(float)num_missing;
242 
243  pdata[j].stn[m].frain[h].qual=5;
244 
245  }
246 
247 
248  }
249 
250  pdata[j].stn[m].frain[4].qual=5;
251 
252  pdata[j].stn[m].frain[4].data=fvalue24;
253 
254 
255 }
256 
257 /* no 24 hour data - estimate done on 6 hour period only*/
258 
259 else {
260 
261  for(h=0;h<4;h++) {
262 
263  if(pdata[j].stn[m].frain[h].qual != 0 &&
264  pdata[j].stn[m].frain[h].qual != 8 &&
265  pdata[j].stn[m].frain[h].qual != 3 &&
266  pdata[j].stn[m].frain[h].qual != 2) {
267 
268  pdata[j].stn[m].frain[h].data=fvalue[h];
269  pdata[j].stn[m].frain[h].qual=5;
270 
271 
272  }
273 
274  }
275 
276 }
277 
278 }
279 
280 }
281 
max_stations
int max_stations
Definition: daily_qc.c:199
isohyets_used
int isohyets_used
Definition: daily_qc.c:192
isoh
struct isoh * isoh
Definition: daily_qc.c:287
isom
int isom
Definition: daily_qc.c:141
method
int method
Definition: daily_qc.c:254
qual
char qual[10]
Definition: display_data.c:29
estimate_missing_stations
void estimate_missing_stations(int j)
Definition: estimate_missing_stations.c:4
prototypes_new.h
fdata
Definition: misc_new.h:84
isoh
Definition: misc.h:465
pdata
Definition: misc.h:216
pdata::stn
struct stn stn[1500]
Definition: misc.h:222
rain::data
float data
Definition: misc.h:193
rain::qual
short int qual
Definition: misc.h:194
station
Definition: misc.h:232
station::lat
float lat
Definition: misc.h:239
station::isoh
float isoh
Definition: misc.h:234
station::lon
float lon
Definition: misc.h:240
station::index
short int index[30]
Definition: misc.h:238
stn
Definition: misc.h:200
stn::frain
struct rain frain[5]
Definition: misc.h:203
Generated by doxygen 1.9.1