1 /* Copyright (C) 1988-2005 by the Institute
2 of Global Environment and Society (IGES).
3
4 See file COPYRIGHT for more information. */
5
6 /* Authored by Jennifer Adams */
7
8 #ifdef HAVE_CONFIG_H
9 #include <config.h>
10
11 /* If autoconfed, only include malloc.h when it's presen */
12 #ifdef HAVE_MALLOC_H
13 #include <malloc.h>
14 #endif
15
16 #else /* undef HAVE_CONFIG_H */
17
18 #include <malloc.h>
19
20 #endif /* HAVE_CONFIG_H */
21
22 #include <stdio.h>
23 #include <string.h>
24 #include <math.h>
25 #include "grads.h"
26 #include "gx.h"
27
28 /* This routine parses a BUFR file, then collects reports that match a data request */
29
getbufr(struct gastn * stn)30 int getbufr (struct gastn *stn) {
31 char ch1[16],ch2[16],*ch,pout[256];
32 char datestr1[25],datestr2[25],datestr3[25];
33 int i,j,k,rc,mintim,maxtim,tim,oflg,keepreport,msgcnt;
34 int t,f,minf,maxf,minft,maxft,verb,toffneg;
35 float minlon,maxlon,minlat,maxlat,minlev,maxlev;
36 float hlon,hlat,hlev,htim;
37 double rfnew,rfold;
38 struct garpt *rpt;
39 struct bufrhdr subsethdr,rfhdr;
40 gabufr_msg *msg;
41 gabufr_val *val, *bval, *rfblockbeg;
42
43 /* For informational messages, set verb=1. For debugging, set verb=2. */
44 verb=1;
45
46 /* Get dimension limits */
47 mintim = stn->tmin; maxtim = stn->tmax;
48 minlon = stn->dmin[0]; maxlon = stn->dmax[0];
49 minlat = stn->dmin[1]; maxlat = stn->dmax[1];
50 minlev = stn->dmin[2]; maxlev = stn->dmax[2];
51 if (minlev > maxlev) {
52 minlev = stn->dmax[2];
53 maxlev = stn->dmin[2];
54 }
55 if (stn->rflag) {
56 minlon = minlon - stn->radius;
57 maxlon = maxlon + stn->radius;
58 minlat = minlat - stn->radius;
59 maxlat = maxlat + stn->radius;
60 }
61 stn->rnum = 0;
62
63 /* Loop through times -- How many data files are we going to have to open? */
64 maxf=maxft=0;
65 minf=minft=stn->pfi->dnum[3];
66 for (tim=mintim; tim<=maxtim; tim++) {
67 if (tim < 1) continue;
68 if (tim > stn->pfi->dnum[3]) break;
69 f = *(stn->pfi->fnums+tim-1); /* f is filenumber for this time */
70 if (f < minf) {
71 minf = f;
72 minft = tim;
73 }
74 if (f > maxf) {
75 maxf = f;
76 maxft = tim;
77 }
78 }
79
80 /* Loop through files */
81 for (f=minf; f<=maxf; f++) {
82
83 /* Find a time axis index that will open file f */
84 for (t=minft; t<=maxft; t++) {
85 if (*(stn->pfi->fnums+t-1) == f) {
86 tim=t;
87 break;
88 }
89 }
90
91 /* Call gaopfn to set the file name */
92 rc = gaopfn(tim,&oflg,stn->pfi);
93 if (rc==-99999) {
94 gaprnt(0,"getbufr error: gaopfn returned -99999\n");
95 goto err;
96 }
97 if (rc==-88888) continue;
98
99 /* Parse the BUFR data if it hasn't already been done */
100 if (!stn->pfi->bufrdset) {
101 gabufr_set_tbl_base_path(gxgnam("tables"));
102 if (stn->pfi->tmplat) {
103 if (verb) {
104 sprintf(pout,"Parsing BUFR file %s\n",stn->pfi->tempname);
105 gaprnt(2,pout);
106 }
107 stn->pfi->bufrdset = gabufr_open(stn->pfi->tempname);
108 }
109 else {
110 if (verb) {
111 sprintf(pout,"Parsing BUFR file %s\n",stn->pfi->name);
112 gaprnt(2,pout);
113 }
114 stn->pfi->bufrdset = gabufr_open(stn->pfi->name);
115 }
116 if (!stn->pfi->bufrdset) {
117 gaprnt(0,"Error from getbufr: gabufr_open failed\n");
118 goto err;
119
120 } else {
121 if (verb) gaprnt(2,"Finished parsing BUFR file\n");
122 }
123 }
124
125 msgcnt = -1;
126 /* Loop through bufr messages looking for valid reports */
127 for (msg = stn->pfi->bufrdset->msgs; msg != NULL; msg = msg->next) {
128 msgcnt++;
129 if (msg->is_new_tbl) continue;
130
131 /* loop through msg subsets */
132 for (i = 0; i < msg->subcnt; i++) {
133
134 /* Copy time vals from the msg header, intialize others */
135 subsethdr.tvals.yr = msg->year;
136 subsethdr.tvals.mo = msg->month;
137 subsethdr.tvals.dy = msg->day;
138 subsethdr.tvals.hr = msg->hour;
139 subsethdr.tvals.mn = msg->min;
140 toffneg = 0;
141 subsethdr.sec = 0;
142 subsethdr.toffvals.yr = 0; /* offset times are initially zero */
143 subsethdr.toffvals.mo = 0;
144 subsethdr.toffvals.dy = 0;
145 subsethdr.toffvals.hr = 0;
146 subsethdr.toffvals.mn = 0;
147 subsethdr.offsec = 0;
148 subsethdr.lon = -999;
149 subsethdr.lat = -999;
150 subsethdr.lev = -999;
151 for (k=0;k<8;k++) *(subsethdr.stid+k)='?';
152
153 /* Look for coordinate values in subset (First loop through subset vals) */
154 getbufrhdr(msg->subs[i], NULL, stn->pfi->bufrinfo, &subsethdr, 0, &toffneg);
155
156 /* Sort gabufr_vals into blocks according to their repetition factors (rf) */
157 val = msg->subs[i]; /* first val in subset */
158 rfblockbeg = val; /* first val in initial rfblock */
159 rfold = val->z; /* initial rf */
160
161 while (1) {
162
163 /* (!val) occurs when all gabufr_vals in the subset have the same rf */
164 /* (val->z != rfold) marks the end of a block of gabufr_vals that have the same rf */
165 if ( (!val) || (val->z != rfold) ) {
166
167 /* rfblock is a set of gabufr_vals that have the same repetition factor */
168 /* Copy bufrhdr values from the subset, look in the rfblock for more */
169 rfhdr.tvals = subsethdr.tvals;
170 rfhdr.toffvals = subsethdr.toffvals;
171 rfhdr = subsethdr;
172 getbufrhdr(rfblockbeg, val, stn->pfi->bufrinfo, &rfhdr, 1, &toffneg);
173
174 /* Determine if we want this report */
175 keepreport=1;
176 /* Convert seconds to minutes and add fields */
177 rfhdr.tvals.mn += rfhdr.sec/60;
178 rfhdr.toffvals.mn += rfhdr.offsec/60;
179
180 /* Merge the time values and time offset values to get report time */
181 if (toffneg) {
182 timsub(&rfhdr.tvals,&rfhdr.toffvals);
183 } else {
184 timadd(&rfhdr.tvals,&rfhdr.toffvals);
185 }
186
187 /* Get the report time in grid coordinates */
188 htim = t2gr(stn->tvals,&rfhdr.toffvals);
189
190 /* Check if time is within range*/
191 if (stn->ftmin==stn->ftmax) {
192 if (fabs(htim-stn->ftmin)>0.5) {
193 keepreport=0;
194 if (verb==2) {
195 printf("report time (%4.1f) is outside range; dmin/dmax=%4.1f tim=%d\n",
196 htim,stn->ftmin,mintim);
197 }
198 }
199 } else {
200 if (htim<stn->ftmin || htim>stn->ftmax) {
201 keepreport=0;
202 if (verb==2) {
203 printf("report time (%4.1f) is outside range; dmin=%4.1f dmax=%4.1f tmin=%d tmax=%d\n",
204 htim,stn->ftmin,stn->ftmax,mintim,maxtim);
205 }
206 }
207 }
208
209 if (keepreport) {
210 if (stn->sflag) {
211 /* check if stids match */
212 for (k=0; k<8; k++) *(ch1+k) = tolower(*(rfhdr.stid+k));
213 for (k=0; k<8; k++) *(ch2+k) = *(stn->stid+k);
214 if (!cmpwrd(ch1,ch2)) {
215 keepreport=0;
216 if (verb==2) printf("report stid doesn't match\n");
217 }
218
219 } else {
220 /* check if stid is still the initialized value */
221 for (k=0;k<8;k++) if (*(rfhdr.stid+k) == '?') {
222 keepreport=0;
223 if (verb==2) printf("report has no stid\n");
224 }
225 /* check if lat and lon are within range */
226 hlon = (float)rfhdr.lon;
227 hlat = (float)rfhdr.lat;
228 if (hlon<minlon) hlon+=360.0;
229 if (hlon>maxlon) hlon-=360.0;
230 if (hlon<minlon || hlon>maxlon || hlat<minlat || hlat>maxlat) {
231 keepreport=0;
232 if (verb==2) printf("report not in lat/lon domain\n");
233 }
234 if (keepreport && stn->rflag &&
235 hypot(hlon-minlon,hlat-minlat)>stn->radius) {
236 keepreport=0;
237 if (verb==2) printf("report not within radius of lat/lon location\n");
238 }
239 }
240 }
241
242 /* loop through rfblock to get a data value */
243 if (keepreport) {
244 for (bval=rfblockbeg; bval != val; bval=bval->next) {
245 if (bval->undef) continue;
246
247 /* Non-replicated surface report */
248 if ((stn->pvar->levels==0) && (bval->z == -1)) {
249
250 /* If variable x,y matches, chain report off stn block */
251 if ((bval->x == stn->pvar->units[0]) && (bval->y == stn->pvar->units[1])) {
252 rpt = gaarpt(stn);
253 if (rpt==NULL) {
254 gaprnt(0,"getbufr error: gaarpt returned NULL\n");
255 goto err;
256 }
257 rpt->lat = hlat;
258 rpt->lon = hlon;
259 rpt->lev = stn->pfi->undef;
260 rpt->tim = htim;
261 rpt->val = bval->val;
262 for (k=0; k<8; k++) *(rpt->stid+k) = *(rfhdr.stid+k);
263 stn->rnum++;
264 break; /* quit loop now that we've got non-replicated report */
265 }
266
267 }
268 /* Replicated surface report */
269 else if ((stn->pvar->levels==2) && (bval->z != -1)) {
270
271 /* If variable x,y matches, chain report off stn block */
272 if ((bval->x == stn->pvar->units[0]) && (bval->y == stn->pvar->units[1])) {
273 rpt = gaarpt(stn);
274 if (rpt==NULL) {
275 gaprnt(0,"getbufr error: gaarpt returned NULL\n");
276 goto err;
277 }
278 rpt->lat = hlat;
279 rpt->lon = hlon;
280 rpt->lev = stn->pfi->undef;
281 rpt->tim = htim;
282 rpt->val = bval->val;
283 for (k=0; k<8; k++) *(rpt->stid+k) = *(rfhdr.stid+k);
284 stn->rnum++;
285 }
286
287 }
288 /* Replicated upper air report */
289 else if ((stn->pvar->levels==1) && (bval->z != -1)) {
290
291 /* check if level is within range */
292 hlev = rfhdr.lev;
293 if (minlev==maxlev) {
294 if (fabs(hlev-minlev)>0.01) {
295 keepreport=0;
296 if (verb==2) printf("report level doesn't match\n");
297 }
298 } else {
299 if (hlev<minlev || hlev>maxlev) {
300 keepreport=0;
301 if (verb==2) printf("report level is out of range\n");
302 }
303 }
304 if (keepreport) {
305 /* If variable x,y matches, chain report off stn block */
306 if (bval->x == stn->pvar->units[0] && bval->y == stn->pvar->units[1]) {
307 rpt = gaarpt (stn);
308 if (rpt==NULL) {
309 gaprnt(0,"getbufr error: gaarpt returned NULL\n");
310 goto err;
311 }
312 rpt->lat = hlat;
313 rpt->lon = hlon;
314 rpt->lev = hlev;
315 rpt->tim = htim;
316 rpt->val = bval->val;
317 for (k=0; k<8; k++) *(rpt->stid+k) = *(rfhdr.stid+k);
318 stn->rnum++;
319 }
320 }
321 } /* Matches if (stn->pvar->levels==0) { ... } else { */
322 } /* Matches for (bval=rfblockbeg; bval != val; bval=bval->next) { */
323 } /* Matches if (keepreport) { */
324
325 /* If we've gotten here then we've reached the end of the subset */
326 if (!val) break;
327
328 /* reset markers */
329 rfblockbeg = val;
330 rfold = val->z;
331 }
332 val = val->next;
333 } /* end of while loop */
334 } /* end of loop through message subsets */
335 } /* end of loop through messages */
336 } /* Matches for (f=minf; f<=maxf; f++) { */
337
338 stn->rpt = sortrpt(stn->rpt);
339 return(0);
340
341 err:
342 for (i=0; i<BLKNUM; i++) {
343 if (stn->blks[i] != NULL) free (stn->blks[i]);
344 }
345 return (1);
346 }
347
getbufrhdr(gabufr_val * first,gabufr_val * last,struct bufrinfo * info,struct bufrhdr * hdr,int flag,int * toffneg)348 void getbufrhdr (gabufr_val *first, gabufr_val *last, struct bufrinfo *info, struct bufrhdr *hdr, int flag, int *toffneg) {
349 int k,toffhr;
350 char bigstr[100];
351 double pval;
352 float tofffrac;
353
354 gabufr_val *val;
355
356 for (val = first; val != last; val = val->next) {
357 if (!val) break;
358 if (val->undef) continue;
359 /* flag should be 0 for subsets, 1 for rfblocks */
360 if ((!flag) && (val->z != -1)) continue;
361
362 /* YEAR */
363 if (val->x == info->base.yrxy[0] && val->y == info->base.yrxy[1]) {
364 if (val->sval == NULL) hdr->tvals.yr = (int)val->val;
365 }
366 if (val->x == info->offset.yrxy[0] && val->y == info->offset.yrxy[1]) {
367 if (val->sval == NULL) hdr->toffvals.yr = (int)val->val;
368 }
369 /* MONTH */
370 if (val->x == info->base.moxy[0] && val->y == info->base.moxy[1]) {
371 if (val->sval == NULL) hdr->tvals.mo = (int)val->val;
372 }
373 if (val->x == info->offset.moxy[0] && val->y == info->offset.moxy[1]) {
374 if (val->sval == NULL) hdr->toffvals.mo = (int)val->val;
375 }
376 /* DAY */
377 if (val->x == info->base.dyxy[0] && val->y == info->base.dyxy[1]) {
378 if (val->sval == NULL) hdr->tvals.dy = (int)val->val;
379 }
380 if (val->x == info->offset.dyxy[0] && val->y == info->offset.dyxy[1]) {
381 if (val->sval == NULL) hdr->toffvals.dy = (int)val->val;
382 }
383 /* HOUR */
384 if (val->x == info->base.hrxy[0] && val->y == info->base.hrxy[1]) {
385 if (val->sval == NULL) hdr->tvals.hr = (int)val->val;
386 }
387 if (val->x == info->offset.hrxy[0] && val->y == info->offset.hrxy[1]) {
388 if (val->sval == NULL) {
389 /* If offset is negative, trip flag and then use absolute value */
390 if (val->val < 0) *toffneg = 1;
391 pval = fabs(val->val);
392 /* If offset contains fractional hours, update minutes too */
393 toffhr = (int)pval;
394 tofffrac = pval - (float)toffhr;
395 hdr->toffvals.hr = toffhr;
396 hdr->toffvals.mn = (int)(0.5+(tofffrac*60.0));
397 }
398 }
399 /* MINUTE */
400 if (val->x == info->base.mnxy[0] && val->y == info->base.mnxy[1]) {
401 if (val->sval == NULL) hdr->tvals.mn = (int)val->val;
402 }
403 if (val->x == info->offset.mnxy[0] && val->y == info->offset.mnxy[1]) {
404 if (val->sval == NULL) hdr->toffvals.mn = (int)val->val;
405 }
406 /* SECONDS */
407 if (val->x == info->base.scxy[0] && val->y == info->base.scxy[1]) {
408 if (val->sval == NULL) hdr->sec = (int)val->val;
409 }
410 if (val->x == info->offset.scxy[0] && val->y == info->offset.scxy[1]) {
411 if (val->sval == NULL) hdr->offsec = (int)val->val;
412 }
413 /* STATION ID */
414 if (val->x == info->stidxy[0] && val->y == info->stidxy[1]) {
415 if (val->sval != NULL) {
416 /* copy string */
417 for (k=0; k<8; k++) {
418 if (*(val->sval+k) == '\0') break;
419 *(hdr->stid+k) = *(val->sval+k);
420 }
421 /* pad with spaces */
422 while (k<8) {
423 *(hdr->stid+k) = ' ';
424 k++;
425 }
426 } else {
427 sprintf(bigstr,"%-10d",(int)val->val);
428 for (k=0; k<8; k++) *(hdr->stid+k) = *(bigstr+k);
429 }
430 }
431 /* LATITUDE */
432 if (val->x == info->latxy[0] && val->y == info->latxy[1]) {
433 if (val->sval == NULL) hdr->lat = val->val;
434 }
435 /* LONGITUDE */
436 if (val->x == info->lonxy[0] && val->y == info->lonxy[1]) {
437 if (val->sval == NULL) hdr->lon = val->val;
438 }
439 /* LEVEL */
440 if (val->x == info->levxy[0] && val->y == info->levxy[1]) {
441 if (val->sval == NULL) hdr->lev = val->val;
442 }
443 }
444 }
445
446 /*
447 * Code for sorting a linked list of station reports so they are
448 * in increasing time order. The algorithm used is Mergesort.
449 * The sort function returns the new head of the list.
450 *
451 * This code is copyright 2001 Simon Tatham.
452 *
453 * Permission is hereby granted, free of charge, to any person
454 * obtaining a copy of this software and associated documentation
455 * files (the "Software"), to deal in the Software without
456 * restriction, including without limitation the rights to use,
457 * copy, modify, merge, publish, distribute, sublicense, and/or
458 * sell copies of the Software, and to permit persons to whom the
459 * Software is furnished to do so, subject to the following
460 * conditions:
461 *
462 * The above copyright notice and this permission notice shall be
463 * included in all copies or substantial portions of the Software.
464 *
465 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
466 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
467 * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
468 * NONINFRINGEMENT. IN NO EVENT SHALL SIMON TATHAM BE LIABLE FOR
469 * ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
470 * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
471 * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
472 * SOFTWARE.
473 *
474 */
sortrpt(struct garpt * list)475 struct garpt * sortrpt(struct garpt *list) {
476 struct garpt *p, *q, *e, *tail;
477 int insize, nmerges, psize, qsize, i;
478
479 if (!list) return NULL;
480 insize = 1;
481 while (1) {
482 p = list;
483 list = NULL;
484 tail = NULL;
485 nmerges = 0; /* count number of merges we do in this pass */
486 while (p) {
487 nmerges++; /* there exists a merge to be done */
488 /* step `insize' places along from p */
489 q = p;
490 psize = 0;
491 for (i = 0; i < insize; i++) {
492 psize++;
493 q = q->rpt;
494 if (!q) break;
495 }
496
497 /* if q hasn't fallen off end, we have two lists to merge */
498 qsize = insize;
499
500 /* now we have two lists; merge them */
501 while (psize > 0 || (qsize > 0 && q)) {
502
503 /* decide whether next rpt of merge comes from p or q */
504 if (psize == 0) {
505 /* p is empty; e must come from q. */
506 e = q;
507 q = q->rpt;
508 qsize--;
509 } else if (qsize == 0 || !q) {
510 /* q is empty; e must come from p. */
511 e = p;
512 p = p->rpt;
513 psize--;
514 } else if ((p->tim - q->tim) <= 0) {
515 /* First rpt of p is lower (or same); e must come from p. */
516 e = p;
517 p = p->rpt;
518 psize--;
519 } else {
520 /* First garpt of q is lower; e must come from q. */
521 e = q;
522 q = q->rpt;
523 qsize--;
524 }
525
526 /* add the next rpt to the merged list */
527 if (tail) {
528 tail->rpt = e;
529 } else {
530 list = e;
531 }
532 tail = e;
533 }
534
535 /* now p has stepped `insize' places along, and q has too */
536 p = q;
537 }
538 tail->rpt = NULL;
539
540 /* If we have done only one merge, we're finished. */
541 if (nmerges <= 1) /* allow for nmerges==0, the empty list case */
542 return list;
543
544 /* Otherwise repeat, merging lists twice the size */
545 insize *= 2;
546 }
547 }
548