Daily_QC
estimate_partial_stations.c
Go to the documentation of this file.
1 
2 #include "prototypes_new.h"
3 
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,k,i,l,ii;
15 double lat1,lon1,fdist,fdata,fvalue[4],lat,lon,testdist,fmult,ftotal,isoh,isoh1,padj,stotal;
16 int num_missing;
17 
18 if(pdata[j].data_time==0)
19  return;
20 
21 /* this routine will estimate 6 hourly periods when 24 hour rain exists */
22 
23 for(m=0;m<max_stations;m++) {
24 
25  /* search for a good 6 hourly period */
26 
27  for(k=0;k<4;k++) {
28 
29  if(pdata[j].stn[m].rain[k].data >= 0 &&
30  (pdata[j].stn[m].frain[k].qual==0 ||
31  pdata[j].stn[m].frain[k].qual==8 ||
32  pdata[j].stn[m].frain[k].qual==3 ||
33  pdata[j].stn[m].frain[k].qual==2))
34  break;
35 
36  }
37 
38  if(k==4) {
39 
40  /* all 6 hourly periods missing or set bad */
41  /* need to re-estimate 24 hour data */
42 
43  if(pdata[j].stn[m].frain[4].data >= 0 &&
44  pdata[j].stn[m].frain[4].qual == 4) {
45 
46  pdata[j].stn[m].frain[4].qual =0;
47  pdata[j].stn[m].frain[4].data =-1;
48  }
49 
50  continue;
51 
52 
53  }
54 
55  /* dont estimate if 24 hour station available */
56 
57  if(pdata[j].stn[m].frain[4].data >= 0 &&
58  pdata[j].stn[m].frain[4].qual != 4 &&
59  pdata[j].stn[m].frain[4].qual != 5)
60  continue;
61 
62  /* at least one missing 6 hourly period found */
63 
64  lat1=station[m].lat;
65  lon1=station[m].lon;
66 
67  /*get isohyet */
68 
69  if(isohyets_used!=0)
70  isoh1=station[m].isoh[isom];
71 
72  /* first */
73 
74  for(k=0;k<4;k++) {
75 
76  fdist=0.0;
77  fdata=0.0;
78 
79  l=0;
80 
81  for(ii=0;ii<30;ii++){
82 
83  i=station[m].index[ii];
84 
85  /* dont estimate unless good or forced good */
86 
87  if(pdata[j].stn[i].frain[k].qual!=0 &&
88  pdata[j].stn[i].frain[k].qual!=8 &&
89  pdata[j].stn[i].frain[k].qual!=3 &&
90  pdata[j].stn[i].frain[k].qual!=6 &&
91  pdata[j].stn[i].frain[k].qual!=2)
92  continue;
93 
94  /*dont use missing stations */
95 
96  if(pdata[j].stn[i].frain[k].data < 0)
97  continue;
98 
99  lat=station[i].lat;
100  lon=station[i].lon;
101 
102  if(isohyets_used != 0)
103  isoh=station[i].isoh[isom];
104 
105  testdist= pow((double)(lat1-lat),2) +
106  pow((double)(lon1-lon),2);
107 
108  if(testdist==0.0)
109  testdist=.000001;
110 
111  if(method==2 && isoh > 0 && isoh1 > 0)
112  padj=pdata[j].stn[i].frain[k].data*isoh1/isoh;
113 
114  else
115  padj=pdata[j].stn[i].frain[k].data;
116 
117  fdist=1/testdist+fdist;
118  fdata=padj/testdist + fdata;
119  l++;
120 
121  }
122 
123  if(l < 5) {
124 
125  fdist=0.0;
126  fdata=0.0;
127 
128  l=0;
129 
130  for(i=0;i<max_stations;i++){
131 
132  if(i==m)
133  continue;
134 
135  if(pdata[j].stn[i].frain[k].qual!=0 &&
136  pdata[j].stn[i].frain[k].qual!=8 &&
137  pdata[j].stn[i].frain[k].qual!=3 &&
138  pdata[j].stn[i].frain[k].qual!=6 &&
139  pdata[j].stn[i].frain[k].qual!=2)
140  continue;
141 
142  if(pdata[j].stn[i].frain[k].data < 0)
143  continue;
144 
145  lat=station[i].lat;
146  lon=station[i].lon;
147 
148  if(isohyets_used != 0)
149  isoh=station[i].isoh[isom];
150 
151  testdist= pow((double)(lat1-lat),2) +
152  pow((double)(lon1-lon),2);
153 
154  if(testdist==0.0)
155  testdist=.0001;
156 
157  if(method==2 && isoh > 0 && isoh1 > 0)
158  padj=pdata[j].stn[i].frain[k].data*isoh1/isoh;
159 
160  else
161  padj=pdata[j].stn[i].frain[k].data;
162 
163  fdist=1/testdist+fdist;
164  fdata=padj/testdist + fdata;
165  l++;
166  }
167 
168  }
169 
170  if(l != 0)
171  fvalue[k]=fdata/fdist;
172 
173  else
174  fvalue[k]=-9999;
175 
176  }
177 
178 
179  /* have good and estimated periods */
180 
181 
182  if(fvalue[0]==-9999 ||
183  fvalue[1]==-9999 ||
184  fvalue[2]==-9999 ||
185  fvalue[3]==-9999)
186  continue;
187 
188  ftotal=0;
189 
190  for(k=0;k<4;k++) {
191 
192  if(pdata[j].stn[m].rain[k].data >= 0 &&
193  (pdata[j].stn[m].frain[k].qual==0 ||
194  pdata[j].stn[m].frain[k].qual==8 ||
195  pdata[j].stn[m].frain[k].qual==3 ||
196  pdata[j].stn[m].frain[k].qual==2))
197  ftotal=ftotal+pdata[j].stn[m].rain[k].data;
198 
199  else {
200 
201  ftotal=ftotal+fvalue[k];
202 
203  /* if(pdata[j].stn[m].frain[k].qual != 1) {*/
204 
205  pdata[j].stn[m].frain[k].data=fvalue[k];
206  pdata[j].stn[m].frain[k].qual=5;
207 
208  /* } */
209 
210  }
211 
212 
213  }
214 
215 
216  pdata[j].stn[m].frain[k].data=ftotal;
217  pdata[j].stn[m].frain[k].qual=4;
218 
219 
220 }
221 
222 
223 }
224 
int max_stations
Definition: daily_qc.c:199
int isohyets_used
Definition: daily_qc.c:192
struct isoh * isoh
Definition: daily_qc.c:287
int isom
Definition: daily_qc.c:141
int method
Definition: daily_qc.c:254
char qual[10]
Definition: display_data.c:29
void estimate_partial_stations(int j)
Definition: misc_new.h:84
Definition: misc.h:465
Definition: misc.h:216
struct stn stn[1500]
Definition: misc.h:222
Definition: misc.h:191
float data
Definition: misc.h:193
short int qual
Definition: misc.h:194
Definition: misc.h:232
float lat
Definition: misc.h:239
float isoh
Definition: misc.h:234
float lon
Definition: misc.h:240
short int index[30]
Definition: misc.h:238
Definition: misc.h:200
struct rain rain[5]
Definition: misc.h:202
struct rain frain[5]
Definition: misc.h:203