1 /* Copyright (C) 1988-2005 by Brian Doty and the Institute
2 of Global Environment and Society (IGES).
3
4 See file COPYRIGHT for more information. */
5
6 /* Authored by B. Doty and 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 "grads.h"
23 #include <math.h>
24 #include <stdio.h>
25 #if USESDF == 1
26 #include <netcdf.h>
27 #if USEHDF == 1
28 #include <mfhdf.h>
29 #endif
30 #endif
31
32 /* Do garead function prototype here */
33 int garead (off_t, int, float *);
34
35 /* Do gafcorlf function prototype here */
36 off_t gafcorlf (int, int, int, int);
37
38 /* global struct for warning level setting */
39 extern struct gamfcmn mfcmn;
40
41 int cvhdr (unsigned char *, struct rpthdr *);
42 int cvflt (unsigned char *, int);
43
44 static char pout[256]; /* Build error msgs here */
45
46 /* Global pointers for this file */
47
48 static struct gafile *pfi;
49 static struct gagrid *pgr;
50 static struct gavar *pvr;
51 static struct gaindx *pindx;
52 static int timerr;
53 static int msgflg=1;
54
55 /* GRIB I/O caching. GRIB data is chached, as well as the bit
56 maps, if present. Sometimes the expanded bit map is cached. */
57
58 static char *cache; /* I/O cache for GRIB */
59 static char *bcache; /* Bit map cache */
60 static int cflag=0; /* cache flag */
61 static int bcflag=0; /* Bit cache flag */
62 static int bpsav = -999; /* Bit cache pointer */
63 static int bssav = -999; /* Bit cache size */
64 static int *bpcach; /* expanded bit map cache */
65
66 /* Station data I/O caching. We will cache fairly small I/O
67 requests that fit within the specified size buffer. If the buffer
68 gets overfilled, we just forget the caching. */
69
70 static int scflg = 0; /* Anything cached? */
71 static int scuca = 0; /* Can use cache for this request */
72 static int scerr = 0; /* Buffer full? */
73 static int scok; /* Ok to fill buffer */
74 static int scpnt; /* Current cache offset */
75 static int scseq; /* File sequence of last request */
76 static struct gastn scstn; /* Previous request */
77 static char *scbuf=NULL; /* Cache */
78 /* Size of cache */
79 #define SCNUM 50000
80
81
82 /* Routine resets flag to allow warning in regards to interpolation */
83
gaiomg()84 void gaiomg () {
85 msgflg = 1;
86 }
87
88 /* Routine to obtain a grid. The addresses of the gagrid
89 structure is passed to this routine. The storage for the
90 grid is obtained and the grid is filled with data. */
91
gaggrd(struct gagrid * pgrid)92 int gaggrd (struct gagrid *pgrid) {
93 int fipnt;
94 float *gr;
95 int x,i,id,jd,d[4],dx[4];
96 float ulo, uhi, undef;
97 int incr,rc,dflag;
98 int size,ssz;
99 #if USESDF == 1
100 /* prototype for SDF reader */
101 int gagsdf(struct gagrid *gridptr, float grid[]) ;
102 #endif
103
104 if (cflag) free(cache);
105 cache = NULL;
106 cflag = 0;
107 if (bcflag) {
108 free(bcache);
109 free(bpcach);
110 }
111 bcache = NULL;
112 bpcach = NULL;
113 bcflag = 0;
114 bssav = -999;
115 bpsav = -999;
116
117
118 pgr = pgrid;
119 pvr = pgr->pvar;
120 pfi = pgr->pfile;
121 timerr = 0;
122 if (pfi->idxflg) pindx = pfi->pindx;
123
124 if (pfi->ppflag && msgflg) {
125 gaprnt (3,"Notice: Automatic Grid Interpolation Taking Place\n");
126 msgflg = 0;
127 }
128
129 if (pfi->type==4) {
130 rc = gagdef();
131 return (rc);
132 }
133
134 /* Check dimensions we were given */
135
136 if (pgr->idim < -1 || pgr->idim > 3 ||
137 pgr->jdim < -1 || pgr->jdim > 3 ||
138 ( pgr->idim == -1 && pgr->jdim!=-1 ) ) {
139 sprintf (pout,"Internal logic check 16: %i %i \n", pgr->idim,
140 pgr->jdim);
141 gaprnt (0,pout);
142 return (16);
143 }
144
145 /* Calc sizes and get storage for the grid */
146
147 id = pgr->idim;
148 jd = pgr->jdim;
149 if (id > -1) pgr->isiz = pgr->dimmax[id] - pgr->dimmin[id] + 1;
150 else pgr->isiz = 1;
151 if (jd > -1) pgr->jsiz = pgr->dimmax[jd] - pgr->dimmin[jd] + 1;
152 else pgr->jsiz = 1;
153
154 size = pgr->isiz*pgr->jsiz;
155 if (size>1) {
156 ssz = size * sizeof(float);
157 gr = (float *)malloc(ssz);
158 if (gr==NULL) {
159 gaprnt (0,"Memory Allocation Error: grid storage \n");
160 return (1);
161 }
162 pgr->grid = gr;
163 } else {
164 pgr->grid = &(pgr->rmin);
165 gr = pgr->grid;
166 }
167
168 /* Handle predefined variable */
169
170 if (pvr->levels<-900) {
171 rc = gagpre();
172 return (rc);
173 }
174
175 for (i=0; i<4; i++) {d[i] = pgr->dimmin[i]; dx[i] = pfi->dnum[i];}
176 dx[2] = pvr->levels;
177 if (dx[2]==0) {
178 if (id==2 || jd==2) goto nozdat;
179 dx[2] = 1;
180 d[2] = 1;
181 }
182
183 incr = pgr->isiz;
184
185 /* If X does not vary, make sure the X coordinate is normalized. */
186
187 if (id!=0 && pfi->wrap) {
188 x=pgr->dimmin[0];
189 while (x<1) x=x+dx[0];
190 while (x>dx[0]) x=x-dx[0];
191 pgr->dimmin[0]=x;
192 pgr->dimmax[0]=x;
193 d[0] = x;
194 }
195
196 /* If any of the non-varying dimensions are out of bounds of the
197 file dimension limits, then we have a grid of missing data.
198 Check for this. */
199
200 for (i=0; i<4; i++) {
201 if (id!=i && jd!=i &&
202 (d[i]<1 || d[i]>dx[i]) ) goto nodat;
203 }
204
205 #if USESDF == 1
206 /* If we don't have a regular FILE*, but we do have an IO_STD*, use SDF read */
207 /* if (!(pfi->infile) && (pfi->sdf_ptr)) { */
208 if (pfi->is_a_SDF) {
209 rc = gagsdf(pgr, gr);
210
211 if (rc == 0) {
212 if (SETMISS) {
213 undef = pgr->undef;
214 ulo = fabs(undef/EPSILON);
215 uhi = undef + ulo;
216 ulo = undef - ulo;
217 for (i = 0; i < size; i++) {
218 if (*gr > ulo && *gr < uhi) {
219 *gr = undef;
220 }
221 gr++;
222 }
223 }
224 }
225 return rc;
226 }
227 #endif
228
229 /* Handle case where X varies. */
230
231 dflag = 0;
232 if ( id == 0 ) {
233 if (jd<0) jd = 1;
234 for (d[jd]=pgr->dimmin[jd]; d[jd]<=pgr->dimmax[jd]; d[jd]++) {
235 if (d[jd]<1 || d[jd]>dx[jd]) {
236 for (i=0; i<incr; i++) *(gr+i) = pgr->undef;
237 } else {
238 rc = gagrow(gr, d);
239 if (rc > 0 ) return (1);
240 if (rc==0) dflag=1;
241 }
242 gr += incr;
243 }
244 if (!dflag) goto nodatmsg;
245 return (0);
246 }
247
248 /* Handle cases where X does not vary. We will have to read
249 each point in the grid seperately. */
250
251 if (jd<0) {
252 if (id<0) { id=0; jd=1; }
253 else jd=0;
254 }
255
256 for (d[jd]=pgr->dimmin[jd]; d[jd]<=pgr->dimmax[jd]; d[jd]++) {
257 if (d[jd]<1 || d[jd]>dx[jd]) {
258 for (i=0; i<incr; i++,gr++) *gr = pgr->undef;
259 } else {
260 for (d[id]=pgr->dimmin[id]; d[id]<=pgr->dimmax[id]; d[id]++) {
261 if (d[id]<1 || d[id]>dx[id]) *gr = pgr->undef;
262 else {
263 rc = garrow (d[0], d[1], d[2], d[3], 1, gr);
264 if (rc != 0 ) return (1);
265 dflag=1;
266 }
267 gr++;
268 }
269 }
270 }
271 if (!dflag) goto nodatmsg;
272 return (0);
273
274 nozdat:
275 if(mfcmn.warnflg>0) {
276 gaprnt (1,"Data Request Warning: Varying Z dimension environment...\n");
277 gaprnt (1," but the requested variable has no Z dimension\n");
278 gaprnt (2," Entire grid contents are set to missing data \n");
279 }
280 for (i=0; i<size; i++,gr++) *gr = pgr->undef; /* yeah, i duped it */
281 return (-1);
282
283 nodat:
284 for (i=0; i<size; i++,gr++) *gr = pgr->undef;
285
286 nodatmsg:
287 if(mfcmn.warnflg>0) {
288 gaprnt (1,"Data Request Warning: Request beyond file limits\n");
289 gaprnt (2," Entire grid contents are set to missing data \n");
290 sprintf (pout," Dimension ranges are: X = %i %i Y = %i %i ",
291 d[0],dx[0],d[1],dx[1]);
292 gaprnt (2,pout);
293 sprintf (pout," Z = %i %i T = %i %i \n",d[2],dx[2],d[3],dx[3]);
294 gaprnt (2,pout);
295 }
296 return (-1);
297
298 }
299
300
301 /* gagrow gets a row of data from the file. The row of data can
302 be 'wrapped' if the x direction of the grid spans the globe. */
303
304
gagrow(float * gr,int * d)305 int gagrow ( float *gr, int *d ) {
306 int fpnt;
307 int rc,i,x,j;
308 int y,z,t;
309
310
311 y = *(d+1);
312 z = *(d+2);
313 t = *(d+3);
314
315 /* If the needed data is within the bounds of the file dimensions
316 then read the data directly. */
317
318 if (pgr->dimmin[0] >= 1 && pgr->dimmax[0] <= pfi->dnum[0]) {
319 rc = garrow (pgr->dimmin[0], y, z, t,
320 (pgr->dimmax[0]-pgr->dimmin[0]+1), gr);
321 if (rc != 0 ) return (1);
322 return (0);
323 }
324
325 /* If the file does not wrap, then read the data directly, if
326 possible. If the requested data lies outside the file's bounds,
327 fill in with missing data where appropriate. */
328
329 if (!pfi->wrap) {
330 if ( pgr->dimmin[0]>=1 && pgr->dimmax[0]<=pfi->dnum[0] ) {
331 rc = garrow (pgr->dimmin[0], y, z, t,
332 (pgr->dimmax[0]-pgr->dimmin[0]+1), gr);
333 if (rc != 0 ) return (1);
334 return (0);
335 }
336
337 for (i=0; i<pgr->isiz; i++) *(gr+i) = pgr->undef;
338 if (pgr->dimmin[0]<1 && pgr->dimmax[0]<1 ) return (-1);
339 if (pgr->dimmin[0]>pfi->dnum[0] &&
340 pgr->dimmax[0]>pfi->dnum[0] ) return (-1);
341 i = 1 - pgr->dimmin[0];
342 if (i>0) gr+=i;
343 i = 1;
344 if (pgr->dimmin[0]>1) i = pgr->dimmin[0];
345 j = pgr->dimmax[0];
346 if (j > pfi->dnum[0]) j = pfi->dnum[0];
347 j = 1 + (j - i);
348 rc = garrow (i, y, z, t, j, gr);
349 if (rc != 0 ) return (1);
350 return (0);
351 }
352
353 /* When the file wraps, we read the entire row into the row buffer, and
354 copy the values as needed into locations in the requested row. */
355 rc = garrow (1, y, z, t, pfi->dnum[0], pfi->rbuf);
356 if (rc != 0 ) return (1);
357
358 for (x=pgr->dimmin[0];x<=pgr->dimmax[0];x++) {
359 i=x;
360 while (i<1) i = i + pfi->dnum[0];
361 while (i>pfi->dnum[0]) i = i-(pfi->dnum[0]); /* Best way??? */
362 *gr = *((pfi->rbuf)+i-1);
363 gr++;
364 }
365 return (0);
366 }
367
gafcor(int x,int y,int z,int t)368 long gafcor ( int x, int y, int z, int t) {
369 off_t pos;
370 long ltmpz,ltmpy,ltmpt;
371 int yy,zz;
372 int levs;
373
374 levs=pvr->levels;
375 if(levs == 0) levs=1;
376 if (pfi->tlpflg) {
377 t = t + pfi->tlpst;
378 if (t > pfi->dnum[3]) t = t - pfi->dnum[3];
379 }
380
381 if (pfi->yrflg) yy = pfi->dnum[1] - y;
382 else yy = y-1;
383
384 if (pfi->zrflg) {
385 if (pvr->levels==0) zz=0;
386 else zz = pvr->levels-z;
387 } else zz = z-1;
388
389 if(pvr->var_t) {
390 pos = (t-1)*(pfi->gsiz)*levs +
391 pvr->offset +
392 zz*(pfi->gsiz) +
393 yy*(pfi->dnum[0]) +
394 (x-1);
395 } else {
396 ltmpt=(long)(t-1)*((long)pfi->tsiz);
397 ltmpz=zz*(pfi->gsiz)*((long)pvr->var_z);
398 ltmpy=yy*(pfi->dnum[0]);
399 pos = ltmpt+
400 (long)pvr->offset +
401 ltmpz +
402 ltmpy +
403 ((long)x-1);
404 }
405
406 if (pfi->xyhdr) pos = pos + pfi->xyhdr;
407 if (pfi->thdr) pos = pos + pfi->thdr;
408 return (pos);
409 }
410
411
412
gafcorlf(int x,int y,int z,int t)413 off_t gafcorlf (int x, int y, int z, int t) {
414 off_t pos;
415 off_t ltmpz,ltmpy,ltmpt;
416 off_t yy,zz;
417 off_t levs;
418 off_t xl, yl, zl, tl;
419
420 xl = x;
421 yl = y;
422 zl = z;
423 tl = t;
424 levs=(off_t)pvr->levels;
425 if(levs == 0) levs=1;
426 if (pfi->tlpflg) {
427 tl = tl + (off_t)pfi->tlpst;
428 if (tl > (off_t)pfi->dnum[3]) tl = tl - (off_t)pfi->dnum[3];
429 }
430
431 if (pfi->yrflg) {
432 yy = (off_t)pfi->dnum[1] - yl;
433 }
434 else {
435 yy = yl - 1;
436 }
437
438 if (pfi->zrflg) {
439 if (levs==0) {
440 zz = 0;
441 }
442 else {
443 zz = levs - zl;
444 }
445 }
446 else {
447 zz = zl - 1;
448 }
449 if (pvr->var_t) {
450 pos = (tl-1)*(off_t)(pfi->gsiz)*levs +
451 (off_t)pvr->offset +
452 zz*(off_t)(pfi->gsiz) +
453 yy*(off_t)(pfi->dnum[0]) +
454 (xl-1);
455 } else {
456 ltmpt=(tl-1)*((off_t)pfi->tsiz);
457 ltmpz=zz*(off_t)(pfi->gsiz)*((off_t)pvr->var_z);
458 ltmpy=yy*((off_t)pfi->dnum[0]);
459 pos = ltmpt +
460 (off_t)pvr->offset +
461 ltmpz +
462 ltmpy +
463 (xl-1);
464 }
465
466 if (pfi->xyhdr) pos = pos + (off_t)pfi->xyhdr;
467 if (pfi->thdr) pos = pos + (off_t)pfi->thdr;
468 return (pos);
469 }
470
gafcyx(int x,int y,int z,int t)471 long gafcyx ( int x, int y, int z, int t) {
472 off_t pos;
473 int yy,zz;
474
475 if (pfi->tlpflg) {
476 t = t + pfi->tlpst;
477 if (t > pfi->dnum[3]) t = t - pfi->dnum[3];
478 }
479
480 if (pfi->yrflg) yy = pfi->dnum[0] - y;
481 else yy = y-1;
482
483 if (pfi->zrflg) {
484 if (pvr->levels==0) zz=0;
485 else zz = pvr->levels-z;
486 } else zz = z-1;
487
488 if(pvr->var_t) {
489 pos = (t-1)*(pfi->gsiz)*(pvr->levels) +
490 pvr->offset +
491 zz*(pfi->gsiz) +
492 (x-1)*(pfi->dnum[1]) +
493 yy;
494 } else {
495 pos = (t-1)*(pfi->tsiz) +
496 pvr->offset +
497 zz*(pfi->gsiz)*(pvr->var_z) +
498 (x-1)*(pfi->dnum[1]) +
499 yy;
500 }
501
502 if (pfi->xyhdr) pos = pos + pfi->xyhdr;
503 if (pfi->thdr) pos = pos + pfi->thdr;
504 return (pos);
505 }
506
507
gardyx(off_t fpos,int len,float * gr)508 int gardyx (off_t fpos, int len, float *gr) {
509 char *ch1,*ch2,*ch3,*ch4,cc1,cc2;
510 int rc,i,j;
511 off_t pos;
512 pos = fpos;
513
514 for (i=1;i<=len;i++) {
515
516 rc = fseeko(pfi->infile, pos*sizeof(float)+pfi->fhdr, 0);
517 if (rc!=0) {
518 gaprnt (0,"Low Level I/O Error: Seek error on data file \n");
519 sprintf (pout,"%d rc=%ld pos=%ld pfi->fhdr =%ld \n",__LINE__,rc,pos,pfi->fhdr);
520 gaprnt (0,pout);
521 sprintf (pout," Data file name = %s \n",pfi->name);
522 gaprnt (0,pout);
523 sprintf (0," Error occurred when seeking to byte %ld \n",fpos);
524 gaprnt (0,pout);
525 return (1);
526 }
527 rc = fread (gr, sizeof(float), 1, pfi->infile);
528 if (rc<1) {
529 gaprnt (0,"Low Level I/O Error: Read error on data file \n");
530 sprintf (pout," Data file name = %s \n",pfi->name);
531 gaprnt (0,pout);
532 sprintf (pout," Error reading %i bytes at location %ld \n",
533 len, fpos);
534 gaprnt (0,pout);
535 return (1);
536 }
537
538 /* Do byte swapping if needed */
539
540 if (pfi->bswap) {
541 ch1 = (char *)gr;
542 ch2 = ch1+1;
543 ch3 = ch2+1;
544 ch4 = ch3+1;
545 for (i=0; i<len; i++) {
546 cc1 = *ch1;
547 cc2 = *ch2;
548 *ch1 = *ch4;
549 *ch2 = *ch3;
550 *ch3 = cc2;
551 *ch4 = cc1;
552 ch1+=4; ch2+=4; ch3+=4; ch4+=4;
553 }
554 }
555
556 /* Set missing data values to exact value if specified */
557 if (SETMISS) {
558 for (j=0;j<len;j++) {
559 if (*gr >= pfi->ulow && *gr <= pfi->uhi) *gr = pgr->undef;
560 gr++;
561 }
562 }
563 pos=pos+pfi->dnum[1];
564
565 }
566
567 return (0);
568
569 }
570
571
572 /* Basic read of a row of data elements -- a row is always
573 in the X direction, which for grads binary is the fastest
574 varying dimension */
575
garrow(int x,int y,int z,int t,int len,float * gr)576 int garrow (int x, int y, int z, int t, int len, float *gr) {
577 int rc,i=0,tt,oflg;
578 off_t fposlf;
579
580 tt = t;
581 if (pfi->tmplat) {
582 tt = gaopfn(t,&oflg,pfi);
583 if (tt==-99999) return(1);
584 if (tt==-88888) {
585 for (i=0; i<len; i++) *(gr+i) = pfi->undef;
586 return (0);
587 }
588 if (oflg) bpsav = -999; /* Force new bit map cache if new file opened */
589 }
590
591 if (pfi->ppflag) { /* Preprojected (pdef) */
592 rc = gaprow (x, y, z, t, tt, len, gr);
593 return (rc);
594 }
595
596 #if USESDF == 1
597 if (pfi->ncflg==1) { /* netcdf */
598 rc = ganrow(x,y,z,tt,len,gr);
599 return(rc);
600 }
601 #endif
602
603 # if USEHDF == 1
604 if (pfi->ncflg==2) { /* HDF-SDS */
605 rc = gahrow(x,y,z,tt,len,gr);
606 return(rc);
607 }
608 #endif
609
610 if (pfi->idxflg) { /* Indexed (grib) */
611 rc = gairow (x, y, z, t, i, len, gr);
612 return (rc);
613 }
614
615 if(pvr->y_x) {
616 fposlf = gafcyx (x,y,z,tt);
617 rc = gardyx (fposlf, len, gr);
618 } else { /* if none of the above... binary */
619 fposlf = gafcorlf (x,y,z,tt);
620 rc = garead (fposlf, len, gr);
621 }
622 return (rc);
623 }
624
garead(off_t fpos,int len,float * gr)625 int garead (off_t fpos, int len, float *gr) {
626 char *ch1,*ch2,*ch3,*ch4,cc1,cc2;
627 unsigned char *uch1,*uch2,ucc1,ucc2;
628 int rc,i;
629 int j,cnt,ival,*ig;
630 float *fg;
631
632 unsigned char *igr;
633 unsigned char *cgr;
634 signed char chsign;
635 unsigned char chunsign;
636 off_t ffpos;
637
638 if(pvr->dfrm == 1) {
639 ffpos = fpos*(off_t)sizeof(char)+(off_t)pfi->fhdr;
640 } else if(pvr->dfrm == 2 || pvr->dfrm == -2 ) {
641 ffpos = fpos*2ll + (off_t)pfi->fhdr;
642 } else if(pvr->dfrm == 4) {
643 ffpos = fpos*(off_t)sizeof(int)+(off_t)pfi->fhdr;
644 } else if(pfi->cray_ieee) {
645 ffpos = fpos*4ll+(off_t)pfi->fhdr;
646 } else {
647 ffpos = fpos*(off_t)sizeof(float)+(off_t)pfi->fhdr;
648 }
649
650 rc = fseeko(pfi->infile, ffpos, 0);
651
652 if (rc!=0) {
653 gaprnt (0,"Low Level I/O Error: Seek error on data file \n");
654 sprintf (pout," Data file name = %s \n",pfi->name);
655 gaprnt (0,pout);
656 sprintf (pout,"%d rc=%d pos=%ld pfi->fhdr =%d\n",__LINE__,rc,fpos,pfi->fhdr);
657 gaprnt (0,pout);
658 sprintf (0," Error occurred when seeking to byte %ld \n",fpos);
659 gaprnt (0,pout);
660 return (1);
661 }
662
663 if(pvr->dfrm == 1) {
664
665 j = len*sizeof(char);
666 igr = (unsigned char *)malloc(len*sizeof(char));
667 if (igr==NULL) {
668 gaprnt (0,"Memory Allocation Error: char grid storage \n");
669 return (1);
670 }
671
672 rc = fread (igr, 1, len, pfi->infile);
673 for(i=0;i<len;i++) {
674 *(gr+i) = (float)(*(igr+i));
675 }
676 free(igr);
677
678 /*mf 961127 handle integer*2 for NCEP CPC does byteswapping too mf*/
679
680 } else if(pvr->dfrm == 2 || pvr->dfrm == -2 ) {
681
682 j = len*2;
683 cgr = (unsigned char *)malloc(len*2);
684 if (cgr==NULL) {
685 gaprnt (0,"Memory Allocation Error: integer*2 storage \n");
686 return (1);
687 }
688
689 rc = fread (cgr, 2, len, pfi->infile);
690 cnt=0;
691
692 /*mf - 961127 - byteswapping mf*/
693
694 if(pfi->bswap) {
695 uch1 = cgr;
696 uch2 = uch1+1;
697 for (i=0; i<len; i++) {
698 ucc1 = *uch1;
699 ucc2 = *uch2;
700 *uch1 = ucc2;
701 *uch2 = ucc1;
702 uch1+=2; uch2+=2;
703 }
704 }
705
706 /*mf - 961127 - signed integer*2 mf*/
707
708 if(pvr->dfrm == -2) {
709
710 for(i=0;i<len;i++) {
711 ival=(int)(*(cgr+cnt)*256) + (int)((*(cgr+cnt+1))) - 65536*(*(cgr+cnt)>>7);
712 *(gr+i) = (float)ival;
713 cnt+=2;
714 }
715
716 /*mf - 961127 - unsigned integer*2 mf*/
717
718 } else {
719
720 for(i=0;i<len;i++) {
721 ival=(int)(*(cgr+cnt)*256) + (int)((*(cgr+cnt+1)));
722 *(gr+i) = (float)ival;
723 cnt+=2;
724 }
725
726 }
727
728 free(cgr);
729
730 } else if(pvr->dfrm == 4) {
731
732 j = len*sizeof(int);
733 ig = (int *)malloc(len*sizeof(int));
734 if (ig==NULL) {
735 gaprnt (0,"Memory Allocation Error: integer*4 storage \n");
736 return (1);
737 }
738
739 rc = fread (ig, sizeof(int), len, pfi->infile);
740
741 for(i=0;i<len;i++) {
742 if (pfi->bswap) {
743 ch1 = (char *)(ig+i);
744 ch2 = ch1+1;
745 ch3 = ch2+1;
746 ch4 = ch3+1;
747 cc1 = *ch1;
748 cc2 = *ch2;
749 *ch1 = *ch4;
750 *ch2 = *ch3;
751 *ch3 = cc2;
752 *ch4 = cc1;
753 }
754 *(gr+i) = (float)(*(ig+i));
755 }
756 free(ig);
757
758 } else if(pfi->cray_ieee) {
759
760 /*mf 970112 -- use wesley ebisuzaki's routines ---mf*/
761
762 cgr = (unsigned char *)malloc(len*4);
763 if (cgr==NULL) {
764 gaprnt (0,"Memory Allocation Error: 32-bit grid storage on cray \n");
765 sprintf (pout,"len*4 = %d\n",len*4);
766 gaprnt (0,pout);
767 return (1);
768 }
769
770 rc = fread (cgr, 4, len, pfi->infile);
771 if (rc<len) {
772 gaprnt (0,"Low Level I/O Error: Read error on 32-bit grid data file on cray \n");
773 sprintf (pout," Data file name = %s \n",pfi->name);
774 gaprnt (0,pout);
775 sprintf (pout," Error reading %i bytes at location %li \n",
776 len, fpos);
777 gaprnt (0,pout);
778 free(cgr);
779 return (1);
780 }
781
782 for(i=0;i<len;i++) {
783
784 if (pfi->bswap) {
785 ch1 = (char *)cgr;
786 ch2 = ch1+1;
787 ch3 = ch2+1;
788 ch4 = ch3+1;
789 cc1 = *ch1;
790 cc2 = *ch2;
791 *ch1 = *ch4;
792 *ch2 = *ch3;
793 *ch3 = cc2;
794 *ch4 = cc1;
795 }
796
797 *(gr+i) = ieee2flt(cgr);
798 cgr+=4;
799
800 }
801
802 free(cgr);
803
804 } else {
805 /* standard direct access read */
806 rc = fread (gr, sizeof(float), len, pfi->infile);
807 }
808
809 if (rc<len) {
810 gaprnt (0,"Low Level I/O Error: Read error on data file \n");
811 sprintf (pout," Data file name = %s \n",pfi->name);
812 gaprnt (0,pout);
813 sprintf (pout," Error reading %i bytes at location %li \n",
814 len, fpos);
815 gaprnt (0,pout);
816 return (1);
817 }
818
819 /* Do byte swapping if needed, only if !CRAY and !cray_ieee*/
820
821 if (pfi->bswap && !GRADS_CRAY && ( pvr->dfrm != 4 ) && ( abs(pvr->dfrm) != 2 ) && ( pvr->dfrm != 1 ) && !pfi->cray_ieee ) {
822 ch1 = (char *)gr;
823 ch2 = ch1+1;
824 ch3 = ch2+1;
825 ch4 = ch3+1;
826 for (i=0; i<len; i++) {
827 cc1 = *ch1;
828 cc2 = *ch2;
829 *ch1 = *ch4;
830 *ch2 = *ch3;
831 *ch3 = cc2;
832 *ch4 = cc1;
833 ch1+=4; ch2+=4; ch3+=4; ch4+=4;
834 }
835 }
836
837 /* Set missing data values to exact value if specified */
838 if (SETMISS) {
839 for (i=0;i<len;i++) {
840 if (*gr >= pfi->ulow && *gr <= pfi->uhi) *gr = pgr->undef;
841 gr++;
842 }
843 }
844
845 return (0);
846 }
847
848 /* Handle a station data request */
849
gagstn(struct gastn * stn)850 int gagstn (struct gastn *stn) {
851 struct garpt *rpt;
852 struct stninf sts;
853 struct rpthdr ghdr, *hdr;
854 int i,j,k,rc,ii,flag,tim,fnum,rtot,rdw,nsiz;
855 int dpos,selflg,oflg;
856 off_t fpos;
857 int sizhdrf,sizhdrd,sizf,sizd,idum;
858 float lnmin,lnmax,ltmin,ltmax,hlon;
859 char ch1[16],ch2[16],*ch;
860 char rec[256];
861
862
863 stn->rpt = NULL;
864 for (i=0; i<BLKNUM; i++) {
865 stn->blks[i] = NULL;
866 }
867
868 if (stn->pfi->bufrflg) { /* bufrstn */
869 rc = getbufr(stn); /* bufrstn */
870 return (rc); /* bufrstn */
871 }
872
873 #if USEGADODS
874 if (stn->pfi->dhandle > -999) { /* dodstn */
875 rc = dodget(stn); /* dodstn */
876 return (rc); /* dodstn */
877 }
878 #endif
879
880 /* Determine cache situation */
881 if (scbuf==NULL && !scerr) {
882 scbuf = (char *)malloc(SCNUM);
883 if (scbuf==NULL) {
884 gaprnt (0,"Memory allocation error: Stn data cache buffer\n");
885 gaprnt (0," Station data cache disabled\n");
886 scerr = 1;
887 }
888 scflg = 0;
889 }
890
891 scuca = 0;
892 scpnt = 0;
893 if (!scerr) scok = 1;
894 if (scflg && stn->pfi->fseq != -999 && scseq == stn->pfi->fseq &&
895 scstn.pfi==stn->pfi && scstn.idim==stn->idim &&
896 scstn.jdim==stn->jdim && scstn.tmin==stn->tmin &&
897 scstn.tmax==stn->tmax && scstn.rflag==stn->rflag &&
898 scstn.sflag==stn->sflag) {
899 rc = 1;
900 if (stn->rflag && scstn.radius!=stn->radius) rc = 0;
901 if (stn->sflag) {
902 for (i=0; i<8; i++) if (scstn.stid[i]!=stn->stid[i]) rc=0;
903 } else {
904 if (scstn.dmin[0]!=stn->dmin[0]) rc = 0;
905 if (scstn.dmin[1]!=stn->dmin[1]) rc = 0;
906 if (scstn.dmax[0]!=stn->dmax[0]) rc = 0;
907 if (scstn.dmax[1]!=stn->dmax[1]) rc = 0;
908 }
909 if (rc) {
910 scuca = 1;
911 scok = 0;
912 }
913 }
914
915 pvr = stn->pvar;
916 pfi = stn->pfi;
917 hdr = &ghdr;
918
919 lnmin = stn->dmin[0]; lnmax = stn->dmax[0];
920 ltmin = stn->dmin[1]; ltmax = stn->dmax[1];
921 if (stn->rflag) {
922 lnmin = lnmin - stn->radius;
923 lnmax = lnmax + stn->radius;
924 ltmin = ltmin - stn->radius;
925 ltmax = ltmax + stn->radius;
926 }
927 stn->rnum = 0;
928
929 /* set size of the file and data hdr */
930 sizhdrf = sizeof(struct rpthdr);
931 sizhdrd = sizeof(struct rpthdr);
932 if(pfi->cray_ieee) {
933 sizhdrf=8+((sizeof(struct rpthdr)-8)/2);
934 }
935
936 /* Loop through times looking for appropriate stations */
937 for (tim=stn->tmin; tim<=stn->tmax; tim++) {
938 if (tim<1) continue;
939 if (tim > pfi->dnum[3]) break;
940
941 if (!scuca && pfi->tmplat) {
942 rc = gaopfn(tim,&oflg,pfi);
943 if (rc==-99999) goto err;
944 if (rc==-88888) {
945 if (scok) {
946 hdr->nlev = 0;
947 gacstn((char *)hdr, NULL, 0, sizhdrd);
948 }
949 continue;
950 }
951 }
952
953 /* Loop through stations for this time looking for valid reports */
954 if (!scuca) {
955 fpos = *(pfi->tstrt+tim-1);
956 rc = gasstn(fpos);
957 if (rc) goto err;
958 }
959
960 while (1) {
961 /* get the header */
962 if (scuca) { /* from the cache */
963 gagcst (sizhdrd, (char *)hdr);
964 } else { /* from the file */
965 if (pfi->seqflg) {
966 rc = garstn(4,(char *)(&rdw),fpos); /*mf changed 4 to sizeof(int) */
967 if (rc) goto err;
968 if (pfi->bswap) gabswp((float *)(&rdw),1);
969 if(pfi->cray_ieee) {
970 idum=be_int2int((unsigned char*)(&rdw));
971 rdw=idum;
972 }
973
974 }
975
976 rc = garstn (sizhdrf, (char *)hdr, fpos);
977 if (rc) goto err;
978 if (pfi->bswap) gahswp(hdr);
979 if(pfi->cray_ieee) {
980 memcpy(rec,hdr,sizhdrf);
981 rc=cvhdr((unsigned char *)rec,hdr);
982 }
983 }
984
985 if (hdr->nlev==0) break; /* END OF DATA CHECK */
986
987 /* Left justify the station id if there are leading blanks */
988
989 j = 0;
990 while (j<7 && hdr->id[0]==' ') {
991 for (i=0; i<7; i++) hdr->id[i] = hdr->id[i+1];
992 hdr->id[7] = ' ';
993 j++;
994 }
995
996 /* Determine if we want to read the data portion of this report */
997 selflg = 1;
998 if (stn->sflag) {
999 getwrd (ch1,hdr->id,8);
1000 lowcas(ch1);
1001 getwrd(ch2,stn->stid,8);
1002 if (!cmpwrd(ch1,ch2)) selflg = 0;
1003 } else {
1004 hlon = hdr->lon;
1005 if (hlon<lnmin) hlon+=360.0;
1006 if (hlon>lnmax) hlon-=360.0;
1007 if (hlon<lnmin || hlon>lnmax ||
1008 hdr->lat<ltmin || hdr->lat>ltmax ) selflg=0;
1009 if (selflg && stn->rflag &&
1010 hypot(hlon-stn->dmin[0],hdr->lat-stn->dmin[1])>stn->radius) {
1011 selflg = 0;
1012 }
1013 }
1014
1015 /* Determine size of the data portion of this report */
1016 if (hdr->flag) {
1017 fnum = (hdr->nlev-1) * (pfi->lvnum+1) + pfi->ivnum;
1018 } else {
1019 fnum = hdr->nlev * (pfi->lvnum+1);
1020 }
1021
1022 /* calc size of floating point data section in the FILE not the machine */
1023 sizd=fnum*sizeof(float);
1024 if(pfi->cray_ieee) {
1025 sizf= fnum*4;
1026 } else {
1027 sizf=fnum*sizeof(float);
1028 }
1029
1030 /* Read the data portion of this report, byteswap it if needed,
1031 and set exact missing data values if specified.*/
1032 if (selflg) {
1033 if (scuca) { /* from the cache */
1034 gagcst (sizd, (char *)pfi->rbuf);
1035 } else { /* from the file */
1036 if (pfi->seqflg) {
1037 ch = (char *)(pfi->rbuf);
1038 nsiz = rdw - sizhdrf;
1039 if (nsiz>0) {
1040 rc = garstn(nsiz,ch,fpos);
1041 if (rc) goto err;
1042 ch += nsiz;
1043 }
1044 rtot = rdw;
1045 nsiz = sizf + sizhdrf;
1046 while (rtot<=nsiz) {
1047 fpos = fpos + rdw + 8;
1048 rc = gasstn(fpos);
1049 if (rc) goto err;
1050 if (rtot==nsiz) break;
1051 rc = garstn(4,(char *)(&rdw),fpos); /*mf changed 4 to sizeof(int) */
1052 if (rc) goto err;
1053 if (pfi->bswap) gabswp((float *)(&rdw),1);
1054 if(pfi->cray_ieee) {
1055 idum=be_int2int((unsigned char*)(&rdw));
1056 rdw=idum;
1057 }
1058 rtot +=rdw;
1059 if (rtot>nsiz) break;
1060 rc = garstn(rdw,ch,fpos);
1061
1062 if (rc) goto err;
1063 idum=rdw*2;
1064 ch += idum;
1065 }
1066 if (rtot>nsiz) {
1067 gaprnt (0,"Low Level I/O Error: Sequential read error\n");
1068 gaprnt (0," Record size exceeds report size\n");
1069 sprintf (pout," Data file name = %s \n",pfi->name);
1070 gaprnt (0,pout);
1071 goto err;
1072 }
1073
1074 /* normal read -- NON sequential */
1075 } else {
1076 rc = garstn (sizf, (char *)pfi->rbuf, fpos+sizhdrf);
1077 fpos = fpos + sizf + sizhdrf;
1078 if (rc) goto err;
1079 }
1080 if (pfi->bswap) gabswp(pfi->rbuf,fnum);
1081 /* convert to float on cray */
1082 if(pfi->cray_ieee) {
1083 rc=cvflt((unsigned char*)pfi->rbuf,fnum);
1084 }
1085 if (SETMISS) {
1086 for (i=0; i<fnum; i++) {
1087 if ((*(pfi->rbuf+i) >= pfi->ulow) && (*(pfi->rbuf+i) <= pfi->uhi))
1088 *(pfi->rbuf+i) = pfi->undef;
1089 }
1090 }
1091 }
1092
1093 /* Check the data portion for any matches. */
1094 rc = gaglvs (tim,hdr,stn);
1095 if (rc) goto err;
1096
1097 /* Cache this report if appropriate */
1098 if (scok) gacstn((char *)hdr,(char *)pfi->rbuf,sizd,sizhdrd);
1099
1100 /* Skip the data portion of this report.*/
1101 } else {
1102 if (scuca) {
1103 gaprnt (0,"Logic Error 8 in gaio\n");
1104 goto err;
1105 }
1106 if (pfi->seqflg) {
1107 rtot = rdw;
1108 sizf += sizhdrf;
1109 while (rtot<=sizf) {
1110 fpos = fpos + rdw + 8;
1111 rc = gasstn(fpos);
1112 if (rc) goto err;
1113 if (rtot==sizf) break;
1114 rc = garstn(4,(char *)(&rdw),fpos); /*mf changed 4 to sizeof(int) */
1115 if (rc) goto err;
1116 if (pfi->bswap) gabswp((float *)(&rdw),1);
1117 if(pfi->cray_ieee) {
1118 idum=be_int2int((unsigned char*)(&rdw));
1119 rdw=idum;
1120 }
1121 rtot +=rdw;
1122 }
1123 if (rtot>sizf) {
1124 gaprnt (0,"Low Level I/O Error: Sequential read error\n");
1125 gaprnt (0," Record size exceeds report size\n");
1126 sprintf (pout," Data file name = %s \n",pfi->name);
1127 gaprnt (0,pout);
1128 goto err;
1129 }
1130 } else {
1131 fpos = fpos + sizf + sizhdrf;
1132 rc = gasstn(fpos);
1133 if (rc) goto err;
1134 }
1135 } /* END OF if (scuca) -- use the cache or not */
1136 } /* END OF while (1) */
1137
1138 if (scok) {
1139 hdr->nlev = 0;
1140 gacstn((char *)hdr, NULL, 0,sizhdrd);
1141 }
1142 }
1143 if (scok) {
1144 scflg = 1;
1145 scstn = *stn;
1146 scseq = stn->pfi->fseq;
1147 } else scflg = 0;
1148 if (scuca) scflg = 1;
1149 return (0);
1150
1151 err:
1152 for (i=0; i<BLKNUM; i++) {
1153 if (stn->blks[i] != NULL) free (stn->blks[i]);
1154 }
1155 return (1);
1156 }
1157
1158 /* Select appropriate variable and levels from report, and chain
1159 them off the stn block. */
1160
gaglvs(int tim,struct rpthdr * hdr,struct gastn * stn)1161 int gaglvs (int tim, struct rpthdr *hdr, struct gastn *stn) {
1162 struct garpt *rpt;
1163 float *vals,*pval;
1164 int i,k,voff,mlev;
1165
1166 vals = pfi->rbuf;
1167 voff = pvr->offset;
1168 if (pvr->levels==0) {
1169 if (hdr->flag) {
1170 pval = vals+voff;
1171 rpt = gaarpt (stn);
1172 if (rpt==NULL) return(1);
1173 rpt->lat = hdr->lat;
1174 rpt->lon = hdr->lon;
1175 rpt->lev = -9.99e33;
1176 rpt->tim = tim + hdr->t;
1177 rpt->val = *pval;
1178 for (k=0; k<8; k++) *(rpt->stid+k) = *(hdr->id+k);
1179 stn->rnum++;
1180 }
1181 } else {
1182 if (hdr->flag) vals = vals + pfi->ivnum;
1183 mlev = hdr->nlev;
1184 if (hdr->flag) mlev--;
1185 for (i=0; i<mlev; i++) {
1186 pval = vals+(i*(pfi->lvnum+1));
1187 if (stn->dmax[2]==stn->dmin[2]) {
1188 if (fabs(*pval-stn->dmin[2])>0.01) continue;
1189 } else {
1190 if (*pval<stn->dmax[2] || *pval>stn->dmin[2]) continue;
1191 }
1192 rpt = gaarpt (stn);
1193 if (rpt==NULL) return(1);
1194 rpt->lat = hdr->lat;
1195 rpt->lon = hdr->lon;
1196 rpt->lev = *pval;
1197 rpt->tim = tim + hdr->t;
1198 rpt->val = *(pval+voff+1);
1199 for (k=0; k<8; k++) *(rpt->stid+k) = *(hdr->id+k);
1200 stn->rnum++;
1201 }
1202 }
1203 return (0);
1204 }
1205
1206 /* Allocate a rpt structure, return pointer to allocated buffer.
1207 On the first request, stn->rpt should be set to NULL. */
1208
gaarpt(struct gastn * stn)1209 struct garpt *gaarpt (struct gastn *stn) {
1210 struct garpt *rpt;
1211 int i;
1212
1213 /* First time through, define the static variables. */
1214
1215 if (stn->rpt == NULL) {
1216 stn->prev = &(stn->rpt);
1217 for (i=0; i<BLKNUM; i++) {
1218 stn->blks[i] = NULL;
1219 }
1220 stn->rptcnt = RPTNUM; /* Force new block allocation */
1221 stn->blkcnt = -1;
1222 }
1223
1224 stn->rptcnt++;
1225 rpt = stn->crpt;
1226
1227 if (stn->rptcnt>=RPTNUM) {
1228 stn->blkcnt++;
1229 if (stn->blkcnt==BLKNUM) {
1230 printf ("Out of memory blocks to allocate \n");
1231 return(NULL);
1232 }
1233 rpt = (struct garpt *)malloc(sizeof(struct garpt)*(RPTNUM+2));
1234 if (rpt==NULL) {
1235 printf ("Couldn't allocate memory for stn block \n");
1236 return(NULL);
1237 }
1238 stn->blks[stn->blkcnt] = rpt;
1239 stn->rptcnt = 0;
1240 } else rpt++;
1241
1242 *(stn->prev) = rpt;
1243 stn->prev = &(rpt->rpt);
1244 rpt->rpt = NULL;
1245 stn->crpt = rpt;
1246 return(rpt);
1247 }
1248
gacstn(char * hdr,char * rdat,int siz,int sizhdr)1249 void gacstn (char *hdr, char *rdat, int siz, int sizhdr) {
1250 int i;
1251 if (scpnt+sizhdr*2+siz+10 > SCNUM) {
1252 scok = 0;
1253 } else {
1254 for (i=0; i<sizhdr; i++) *(scbuf+scpnt+i) = *(hdr+i);
1255 scpnt += sizhdr;
1256 if (siz>0) {
1257 for (i=0; i<siz; i++) *(scbuf+scpnt+i) = *(rdat+i);
1258 scpnt += siz;
1259 }
1260 }
1261 }
1262
1263 /* Return info from the station data cache */
1264
gagcst(int siz,char * ch)1265 void gagcst (int siz, char *ch) {
1266 int i;
1267 for (i=0; i<siz; i++) *(ch+i) = *(scbuf+scpnt+i);
1268 scpnt += siz;
1269 }
1270
1271 /* Seek to specified location in a station data file */
1272
gasstn(off_t fpos)1273 int gasstn (off_t fpos) {
1274 int rc;
1275
1276 rc = fseeko(pfi->infile, fpos, 0);
1277 if (rc!=0) {
1278 gaprnt (0,"Low Level I/O Error: Seek error on data file \n");
1279 sprintf (pout," Data file name = %s \n",pfi->name);
1280 gaprnt (0,pout);
1281 sprintf (pout,"%d rc=%ld pos=%ld pfi->fhdr =%ld \n",__LINE__,rc,fpos,pfi->fhdr);
1282 gaprnt (0,pout);
1283 sprintf (pout," Error occurred when seeking to byte %li \n",fpos);
1284 gaprnt (0,pout);
1285 return (1);
1286 }
1287 return (0);
1288 }
1289
1290 /* Read specified amount of data from a station data file */
1291
garstn(int siz,char * val,off_t fpos)1292 int garstn (int siz, char *val, off_t fpos) {
1293
1294 int rc;
1295
1296 rc = fread (val, siz, 1, pfi->infile);
1297 if (rc<1) {
1298 gaprnt (0,"Low Level I/O Error: Read error on data file \n");
1299 sprintf (pout," Data file name = %s \n",pfi->name);
1300 gaprnt (0,pout);
1301 sprintf (pout," Error reading %i bytes at location %li \n",
1302 siz, fpos);
1303 gaprnt (0,pout);
1304 return (1);
1305 }
1306 return (0);
1307 }
1308
1309 /* Obtain user requested grid from defined variable */
1310
gagdef(void)1311 int gagdef (void) {
1312 int id, jd, i, flag;
1313 int ys,zs,ts,siz,pos;
1314 int d[4],d1min,d1max,d2min,d2max,xt,yt;
1315 float *v;
1316
1317 /* If a dimension is a fixed dimension in the defined
1318 variable, it must be a fixed dimension in the output
1319 grid. */
1320
1321 id = pgr->idim;
1322 jd = pgr->jdim;
1323 if (jd>-1) {
1324 if (pfi->dnum[jd]==1) {
1325 jd = -1;
1326 pgr->jdim = -1;
1327 pgr->jsiz = -1;
1328 }
1329 }
1330 if (id>-1) {
1331 if (pfi->dnum[id]==1) {
1332 id = jd;
1333 pgr->idim = pgr->jdim;
1334 pgr->isiz = pgr->jsiz;
1335 pgr->igrab = pgr->jgrab;
1336 pgr->iabgr = pgr->jabgr;
1337 pgr->ivals = pgr->jvals;
1338 pgr->iavals = pgr->javals;
1339 pgr->ilinr = pgr->jlinr;
1340 jd = -1;
1341 pgr->jdim = -1;
1342 pgr->jsiz = 1;
1343 }
1344 }
1345
1346 /* Set up constants for array subscripting */
1347
1348 ys = pfi->dnum[0];
1349 zs = ys * pfi->dnum[1];
1350 ts = zs * pfi->dnum[2];
1351
1352 /* Set up dimension ranges */
1353
1354 for (i=0; i<4; i++) d[i] = pgr->dimmin[i] - pfi->dimoff[i] - 1;
1355 for (i=0; i<4; i++) if (pfi->dnum[i]==1) d[i] = 0;
1356 if (id>-1) {
1357 d1min = d[id];
1358 d1max = pgr->dimmax[id] - pfi->dimoff[id] - 1;
1359 }
1360 if (jd>-1) {
1361 d2min = d[jd];
1362 d2max = pgr->dimmax[jd] - pfi->dimoff[jd] - 1;
1363 }
1364
1365 /* Get storage for output grid */
1366
1367 pgr->isiz = 1;
1368 pgr->jsiz = 1;
1369 if (id>-1) pgr->isiz = 1 + d1max - d1min;
1370 if (jd>-1) pgr->jsiz = 1 + d2max - d2min;
1371 siz = pgr->isiz * pgr->jsiz;
1372 if (siz>1) {
1373 pgr->grid = (float *)malloc(sizeof(float)*siz);
1374 if (pgr->grid==NULL) {
1375 gaprnt (0,"Memory Allocation Error: Grid Request\n");
1376 return (2);
1377 }
1378 } else {
1379 pgr->grid = &(pgr->rmin);
1380 }
1381
1382 /* Normalize time coordinate if not varying */
1383 /* This does not handle leap years properly!!!! Gotta fix this
1384 someday */
1385
1386 if (pfi->climo && id!=3 && jd!=3) clicyc(d+3);
1387
1388 /* Check for entirely undefined grid */
1389
1390 flag = 0;
1391 for (i=0; i<4; i++) {
1392 if (i!=id && i!=jd && (d[i]<0 || d[i]>=pfi->dnum[i])) flag = 1;
1393 }
1394 if (flag) {
1395 for (i=0; i<siz; i++) *(pgr->grid+i) = pfi->undef;
1396 return (0);
1397 }
1398
1399 /* Move appropriate grid values */
1400
1401 if (id==-1 && jd==-1) {
1402 pos = d[0] + d[1]*ys + d[2]*zs + d[3]*ts;
1403 pgr->rmin = *(pfi->rbuf+pos);
1404 return (0);
1405 }
1406
1407 v = pgr->grid;
1408
1409 if (jd==-1) {
1410 for (xt=d1min; xt<=d1max; xt++) {
1411 d[id] = xt;
1412 if (id==3 && pfi->climo) clicyc(d+3);
1413 if (d[id]<0 || d[id]>=pfi->dnum[id]) *v = pfi->undef;
1414 else {
1415 pos = d[0] + d[1]*ys + d[2]*zs + d[3]*ts;
1416 *v = *(pfi->rbuf+pos);
1417 }
1418 v++;
1419 }
1420 return (0);
1421 }
1422
1423 for (yt=d2min; yt<=d2max; yt++) {
1424 d[jd] = yt;
1425 if (jd==3 && pfi->climo) clicyc(d+3);
1426 for (d[id]=d1min; d[id]<=d1max; d[id]++) {
1427 if (d[jd]<0 || d[jd]>=pfi->dnum[jd] ||
1428 d[id]<0 || d[id]>=pfi->dnum[id]) *v = pfi->undef;
1429 else {
1430 pos = d[0] + d[1]*ys + d[2]*zs + d[3]*ts;
1431 *v = *(pfi->rbuf+pos);
1432 }
1433 v++;
1434 }
1435 }
1436 return(0);
1437 }
1438
clicyc(int * ti)1439 void clicyc (int *ti) {
1440 struct dt dtim,otim;
1441 float tt;
1442
1443 if (pfi->climo>0) {
1444 while (*ti>pfi->cysiz-1) *ti = *ti - pfi->cysiz;
1445 while (*ti<0) *ti = *ti + pfi->cysiz;
1446 }
1447 }
1448
1449 /* Fill in grid for predefined variable */
1450
gagpre(void)1451 int gagpre (void) {
1452 float (*conv)(float *, float);
1453 int d[4],id,jd,i,j,dim;
1454 float *gr,*vals,t;
1455
1456 id = pgr->idim;
1457 jd = pgr->jdim;
1458 for (i=0; i<4; i++) d[i] = pgr->dimmin[i];
1459
1460 dim = pvr->offset;
1461 conv = pfi->gr2ab[dim];
1462 vals = pfi->grvals[dim];
1463
1464 gr = pgr->grid;
1465
1466 if (id>-1 && jd>-1) {
1467 for (d[jd]=pgr->dimmin[jd]; d[jd]<=pgr->dimmax[jd]; d[jd]++) {
1468 for (d[id]=pgr->dimmin[id]; d[id]<=pgr->dimmax[id]; d[id]++) {
1469 t = (float)(d[dim]);
1470 *gr = conv(vals, t);
1471 gr++;
1472 }
1473 }
1474 } else if (id>-1) {
1475 for (d[id]=pgr->dimmin[id]; d[id]<=pgr->dimmax[id]; d[id]++) {
1476 t = (float)(d[dim]);
1477 *gr = conv(vals, t);
1478 gr++;
1479 }
1480 } else {
1481 t = (float)(d[dim]);
1482 *gr = conv(vals, t);
1483 }
1484 return (0);
1485 }
1486
1487
1488 /* Read index data, in this case GRIB type data.
1489 Currently assumes no pole point, and only one record
1490 per grid. */
1491
gairow(int x,int y,int z,int t,int offset,int len,float * gr)1492 int gairow (int x, int y, int z, int t, int offset, int len, float *gr) {
1493 int irec,ioff,bstrt,bend,blen,cstrt,cend,clen,rc;
1494 int ival,i,yy,bpos,boff,siz,gtyp,xsiz,ysiz;
1495 off_t fpos;
1496 float dsf,bsf,ref;
1497
1498 /* Figure out position and length of the I/O */
1499
1500 gtyp = *(pindx->hipnt+3);
1501 irec = (t-1)*pfi->trecs + pvr->recoff + z - 1;
1502 if (pfi->ppflag) {xsiz = pfi->ppisiz; ysiz = pfi->ppjsiz;}
1503 else {xsiz = pfi->dnum[0]; ysiz = pfi->dnum[1];}
1504 if (gtyp==29) {
1505 xsiz = 145;
1506 irec = irec*6;
1507 if (y<37) y--;
1508 else {irec+=3; y-=37; }
1509 yy = y;
1510 } else {
1511 irec = irec*3;
1512 if (pfi->yrflg) yy = ysiz - y;
1513 else yy = y-1;
1514 }
1515 if (pfi->ppflag) ioff = offset;
1516 else ioff = yy*xsiz + x - 1;
1517 boff = ioff;
1518 blen = *(pindx->intpnt + irec + 2);
1519 if (blen<0) {
1520 for (i=0; i<len; i++) *(gr+i) = pfi->undef;
1521 return (0);
1522 }
1523 bpos = *(pindx->intpnt + irec + 1);
1524 dsf = *(pindx->fltpnt+irec);
1525 bsf = *(pindx->fltpnt+irec+1);
1526 ref = *(pindx->fltpnt+irec+2);
1527 if (bpos>-900 && bpos!=bpsav) {
1528 bpsav = bpos;
1529 siz = 2+(xsiz*ysiz)/8;
1530 if (siz>bssav) {
1531 if (bcflag) {
1532 free(bcache);
1533 free(bpcach);
1534 }
1535 bcache = (char *)malloc(siz);
1536 bpcach = (int *)malloc(sizeof(int)*(xsiz*ysiz+1));
1537 if (bcache==NULL||bpcach==NULL) {
1538 gaprnt(0,"Memory Allocation Error During GRIB I/O\n");
1539 return (1);
1540 }
1541 bssav = siz;
1542 bcflag = 1;
1543 }
1544 rc = fseeko(pfi->infile, bpos, 0);
1545 rc = fread(bcache,1,siz,pfi->infile);
1546 if (rc!=siz) {
1547 gaprnt(0,"GRIB I/O Error: Bit Map I/O\n");
1548 return(1);
1549 }
1550 boff=1;
1551 for (i=0; i<xsiz*ysiz; i++) {
1552 if (gagbb(bcache,i,1)) {
1553 *(bpcach+i) = boff;
1554 boff++;
1555 } else {
1556 *(bpcach+i) = -1*boff;
1557 }
1558 }
1559 *(bpcach+xsiz*ysiz) = boff; /* Provide an ending offset */
1560 }
1561 if (bpos>-900) {
1562 boff = *(bpcach+ioff);
1563 if (boff<0) boff = -1*boff;
1564 boff--;
1565 bstrt = blen * boff;
1566 boff = *(bpcach+ioff+len);
1567 if (boff<0) boff = -1*boff;
1568 boff--;
1569 bend = blen * boff - 1;
1570 } else {
1571 bstrt = blen * boff;
1572 bend = bstrt + blen*len;
1573 }
1574 cstrt = bstrt/8;
1575 cend = bend/8;
1576 clen = cend-cstrt+2;
1577 fpos = *(pindx->intpnt+irec);
1578 rc = gaird(fpos,cstrt,clen,xsiz,ysiz,blen);
1579 if (rc) return(rc);
1580 bstrt = bstrt - cstrt*8;
1581 for (i=0; i<len; i++) {
1582 if (bpos>-900 && *(bpcach+ioff+i)<0) *(gr+i) = pfi->undef;
1583 else {
1584 ival = gagbb(pfi->pbuf,bstrt,blen);
1585 *(gr+i) = ( ref + (float)ival * bsf )/dsf;
1586 bstrt += blen;
1587 }
1588 }
1589 return (0);
1590 }
1591
gaird(off_t fpos,int cstrt,int clen,int xsiz,int ysiz,int blen)1592 int gaird (off_t fpos, int cstrt, int clen, int xsiz, int ysiz, int blen) {
1593 int rc,siz,i;
1594 if (pfi->ppflag && pgr->idim==0 && pgr->jdim==1) {
1595 if (!cflag) {
1596 cflag = 1;
1597 siz = 5 + xsiz*ysiz*blen/8; /* qqq Warning: siz calc does not */
1598 /* qqq take into account bms!!! */
1599 cache = (char *)malloc(siz);
1600 if (cache==NULL) {
1601 gaprnt(0,"GRIB Memory Allocation Error\n");
1602 return (1);
1603 }
1604 rc = fseeko(pfi->infile, fpos, 0);
1605 rc = fread(cache,sizeof(char),siz,pfi->infile);
1606 if (rc==0) {
1607 sprintf (pout,"GRIB I/O Error reading %i bytes at %ld\n",siz,fpos); /* xxx */
1608 gaprnt (0,pout);
1609 gaprnt (0," File name is: ");
1610 if (pfi->tempname) gaprnt(0,pfi->tempname);
1611 else gaprnt(0,pfi->name);
1612 gaprnt (0,"\n");
1613 return (1);
1614 }
1615 }
1616 if (cache==NULL) return(1);
1617 for (i=0; i<clen; i++) {
1618 *(pfi->pbuf+i) = *(cache+cstrt+i);
1619 }
1620 } else {
1621 rc = fseeko(pfi->infile, fpos+cstrt, 0);
1622 rc = fread (pfi->pbuf, sizeof(char), clen, pfi->infile);
1623 if (rc==0) {
1624 sprintf (pout,"GRIB I/O Error reading %i bytes at %ld\n",clen,fpos+cstrt);
1625 gaprnt (0,pout);
1626 gaprnt (0," File name is: ");
1627 if (pfi->tempname) gaprnt(0,pfi->tempname);
1628 else gaprnt(0,pfi->name);
1629 gaprnt (0,"\n");
1630 return(1);
1631 }
1632 }
1633 return(0);
1634 }
1635
1636 /* Routine to open appropriate file when using file templates */
1637 /* Warning -- changes time value to time with respect to this file */
1638
gaopfn(int t,int * oflg,struct gafile * pfi)1639 int gaopfn(int t, int *oflg, struct gafile *pfi) {
1640 int i,rc;
1641 struct dt dtim, dtimi;
1642 char *fn;
1643
1644 *oflg = 0;
1645 if (t<1 || t>pfi->dnum[3]) return(-99999);
1646 i = *(pfi->fnums+t-1);
1647
1648 /* The current file is not the one we need; close it */
1649 if (i != pfi->fnumc) {
1650 /* close SDF file */
1651 if (pfi->ncflg) {
1652 #if USESDF == 1
1653 if (pfi->ncflg==1) {
1654 if (pfi->ncid != -999) ncclose(pfi->ncid); /* no error checking */
1655 }
1656 #if USEHDF == 1
1657 else if (pfi->ncflg==2) {
1658 if (pfi->sdid != -999) SDend(pfi->sdid); /* no error checking */
1659 }
1660 #endif
1661 #endif
1662 }
1663 /* close BUFR file*/
1664 else if (pfi->bufrflg) {
1665 if (pfi->bufrdset) {
1666 gabufr_close(pfi->bufrdset); /* free memory */
1667 pfi->bufrdset=NULL; /* reset the pointer */
1668 }
1669 }
1670 /* close non-SDF, non-BUFR file */
1671 else {
1672 if (pfi->infile!=NULL) fclose(pfi->infile);
1673 }
1674
1675 /* find the filename that goes with current time */
1676 if (pfi->tempname!=NULL) free(pfi->tempname);
1677 gr2t(pfi->grvals[3], (float)t, &dtim);
1678 gr2t(pfi->grvals[3], 1.0, &dtimi);
1679 fn = gafndt(pfi->name, &dtim, &dtimi, pfi->abvals[3], pfi->pchsub1, t);
1680 if (fn==NULL) return (-99999);
1681
1682 /* Open the data file */
1683 rc = 0;
1684 pfi->tempname = fn;
1685 pfi->fnumc = i;
1686 /* open netcdf */
1687 if (pfi->ncflg==1) {
1688 rc = gaopnc (pfi,1,0);
1689 if (rc) pfi->ncid = -999;
1690 }
1691 /* open hdfsds */
1692 else if (pfi->ncflg==2) {
1693 rc = gaophdf (pfi,1,0);
1694 if (rc) pfi->sdid = -999;
1695 }
1696 /* open all others except BUFR */
1697 else if (!pfi->bufrflg) {
1698 pfi->infile = fopen (fn, "rb");
1699 if (pfi->infile == NULL) rc = 1;
1700 }
1701
1702 /* Error checking on file open */
1703 if (rc) {
1704 if (pfi->errflg && timerr!=t) {
1705 gaprnt(1,"Warning: Open error, fn = ");
1706 gaprnt(1,fn);
1707 gaprnt(1,"\n");
1708 timerr = t;
1709 }
1710 pfi->fnumc = 0;
1711 return(-88888);
1712 }
1713 *oflg = 1;
1714 }
1715 t = 1 + t - pfi->fnumc;
1716 return (t);
1717 }
1718
1719 /* Read in a row of data from a pre-projected grid data set.
1720 This involves doing interpolation to the lat-lon
1721 grid */
1722
gaprow(int x,int y,int z,int t,int tt,int len,float * gr)1723 int gaprow (int x, int y, int z, int t, int tt, int len, float *gr) {
1724 float p[4],dx,dy,g1,g2;
1725 float vals[9],wts[9],sum,wt;
1726 int ioffs[9],cnt,ig0,goflg;
1727 int rc,i,j,ig,ioff,ix,iy,ncig,ncjg;
1728 off_t pos,pos0;
1729
1730 /* Handle generalized arbitrary points + weights */
1731
1732 if (pfi->ppflag==8) {
1733 cnt = (int)(pfi->ppvals[0]+0.1);
1734 pos0 = (tt-1)*(pfi->tsiz) + pvr->offset + (z-1)*(pfi->gsiz);
1735 ig0 = (y-1) * pfi->dnum[0] + x - 1;
1736 for (i=0; i<len; i++) {
1737 ig = ig0 + i;
1738 goflg = 0;
1739 for (j=0; j<cnt; j++) {
1740 ioffs[j] = *(pfi->ppi[j]+ig);
1741 if (ioffs[j] >= 0 && ioffs[j] <= pfi->gsiz) goflg = 1;
1742 wts[j] = *(pfi->ppf[j]+ig);
1743 }
1744 if (!goflg) *gr = pgr->undef;
1745 else {
1746 goflg = 1;
1747 j = 0;
1748 sum = 0.0; wt = 0.0;
1749 while (j<cnt) {
1750 if (ioffs[j] > 0) {
1751 pos = pos0 + ioffs[j] - 1;
1752 if (pfi->idxflg) {
1753 rc = gairow(x,y,z,t,ioffs[j],1,vals+j);
1754 }
1755 else if (pfi->ncflg==1) {
1756 /* The y-argument to ganrow/gahrow is hardcoded to 1 for PDEF file option */
1757 rc = ganrow(ioffs[j],1,z,tt,1,vals+j); /* PDEF FILE for netcdf */
1758 }
1759 else if (pfi->ncflg==2) {
1760 rc = gahrow(ioffs[j],1,z,tt,1,vals+j); /* PDEF FILE for hdfsds */
1761 }
1762 else
1763 rc = garead(pos,1,vals+j);
1764 if (rc) return(rc);
1765 if ( *(vals+j)==pgr->undef) {
1766 goflg = 0;
1767 break;
1768 }
1769 sum = sum + *(vals+j) * *(wts+j);
1770 wt = wt + *(wts+j);
1771 }
1772 j++;
1773 }
1774 if (goflg && wt!=0.0) *gr = sum/wt;
1775 else *gr = pgr->undef;
1776 }
1777 gr++;
1778 }
1779 } else {
1780 for (i=0; i<len; i++) {
1781 ig = (y-1) * pfi->dnum[0] + x + i - 1;
1782 ioff = *(pfi->ppi[0]+ig);
1783 if (ioff<0) *gr = pgr->undef;
1784 else {
1785 dx = *(pfi->ppf[0]+ig);
1786 dy = *(pfi->ppf[1]+ig);
1787 pos = (tt-1)*(pfi->tsiz) + pvr->offset + (z-1)*(pfi->gsiz) + ioff;
1788
1789 /* Get the first two pre-projected grid values */
1790 if (pfi->idxflg) {
1791 rc = gairow(x,y,z,t,ioff,2,p);
1792 }
1793 else if (pfi->ncflg==1) {
1794 ncig = (int)(1 + ioff%pfi->ppisiz);
1795 ncjg = (int)(1 + ioff/pfi->ppisiz);
1796 rc = ganrow(ncig,ncjg,z,tt,2,p); /* PDEF BILIN for netcdf */
1797 }
1798 else if (pfi->ncflg==2) {
1799 gaprnt(0,"PDEF not fully implemented for dtype hdfsds\n");
1800 rc = 1;
1801 }
1802 else {
1803 rc = garead(pos,2,p);
1804 }
1805 if (rc) return(rc);
1806
1807 /* Get the second two pre-projected grid values */
1808 if (pfi->idxflg) {
1809 rc = gairow(x,y,z,t,ioff+pfi->ppisiz,2,p+2);
1810 }
1811 else if (pfi->ncflg==1) {
1812 ncjg++;
1813 rc = ganrow(ncig,ncjg,z,tt,2,p+2); /* PDEF BILIN for netcdf */
1814 }
1815 else if (pfi->ncflg==2) {
1816 gaprnt(0,"PDEF not fully implemented for dtype hdfsds\n");
1817 rc = 1;
1818 }
1819 else {
1820 rc = garead(pos+pfi->ppisiz,2,p+2);
1821 }
1822 if (rc) return(rc);
1823
1824 /* Do the bilinear interpolation, as long as we have no undefs */
1825 if (p[0]==pgr->undef || p[1]==pgr->undef ||
1826 p[2]==pgr->undef || p[3]==pgr->undef) *gr = pgr->undef;
1827 else {
1828 g1 = p[0] + (p[1]-p[0])*dx;
1829 g2 = p[2] + (p[3]-p[2])*dx;
1830 *gr = g1 + (g2-g1)*dy;
1831 }
1832 }
1833 gr++;
1834 }
1835 }
1836 return(0);
1837 }
1838
1839 /* Convert header */
1840
cvhdr(unsigned char * rec,struct rpthdr * hdr)1841 int cvhdr (unsigned char *rec, struct rpthdr *hdr) {
1842 int rc;
1843 memcpy(hdr->id,&rec[0],8);
1844 hdr->lat=ieee2flt(&rec[8]);
1845 hdr->lon=ieee2flt(&rec[12]);
1846 hdr->t=ieee2flt(&rec[16]);
1847 hdr->nlev=be_int2int(&rec[20]);
1848 hdr->flag=be_int2int(&rec[24]);
1849 return (0);
1850 }
1851
1852
cvflt(unsigned char * rec,int fnum)1853 int cvflt(unsigned char *rec, int fnum) {
1854 float *val;
1855 int i;
1856 val=(float *)malloc(sizeof(float)*fnum);
1857 for(i=0;i<fnum;i++) {
1858 *(val+i)=ieee2flt(&rec[i*4]);
1859 }
1860 memcpy(rec,val,fnum*sizeof(float));
1861 free(val);
1862 return (0);
1863 }
1864
1865
1866 /* Open a netcdf file */
gaopnc(struct gafile * pfil,int tflag,int eflag)1867 int gaopnc (struct gafile *pfil, int tflag, int eflag) {
1868 #if USESDF == 1
1869 struct gavar *pvarl;
1870 int i,rc;
1871
1872 #if USEHDF == 1
1873 if (tflag) {
1874 i = ncopen (pfil->tempname, NC_NOWRITE);
1875 }
1876 else {
1877 i = ncopen (pfil->name, NC_NOWRITE);
1878 }
1879 if (i == -1) {
1880 if (eflag) {
1881 gaprnt(0,"NetCDF Error (ncopen): Unable to open file\n");
1882 }
1883 return(1);
1884 }
1885 #else
1886 if (tflag) {
1887 rc = nc_open(pfil->tempname, NC_NOWRITE, &i);
1888 }
1889 else {
1890 rc = nc_open(pfil->name, NC_NOWRITE, &i);
1891 }
1892 if (rc != NC_NOERR) {
1893 if (eflag) {
1894 gaprnt(0,"NetCDF Error (nc_open): Unable to open file\n");
1895 }
1896 return (1);
1897 }
1898 #endif
1899
1900 pfil->ncid = i;
1901 pvarl = pfil->pvar1;
1902 for (i=0; i<pfil->vnum; i++) {
1903 if (pvarl->ncvid != -888) pvarl->ncvid = -999;
1904 pvarl++;
1905 }
1906 return (0);
1907 #else
1908 return(1);
1909 #endif
1910 }
1911
1912
1913 /* Read a row varying in the X direction from a netcdf grid */
ganrow(int x,int y,int z,int t,int len,float * gr)1914 int ganrow (int x, int y, int z, int t, int len, float *gr) {
1915 #if USESDF == 1
1916 int vid,error,rc,dim,i,j,yy,zz;
1917 size_t start[16],count[16];
1918 nc_type data_dtype,attr_dtype;
1919 short *sval,sundef,sscale,sadd;
1920 long *lval,lundef,lscale,ladd;
1921 double dscale,dadd;
1922 float val,ulow,uhi;
1923 int oldncopts ; /* to save and restore setting for automatic error handling */
1924
1925 /* Turn off automatic error handling. */
1926 ncopts = NC_VERBOSE ;
1927 oldncopts = ncopts ;
1928 ncopts = 0;
1929
1930 /* Get the varid if we haven't already done that for this file */
1931 if (pvr->ncvid == -888) {
1932 ncopts = oldncopts ;
1933 return(1); /* already tried and failed */
1934 }
1935 if (pvr->ncvid == -999) {
1936 error=0;
1937 #if USEHDF == 1
1938 if (pvr->longnm) {
1939 vid = ncvarid(pfi->ncid, pvr->longnm);
1940 }
1941 else {
1942 vid = ncvarid(pfi->ncid, pvr->abbrv);
1943 }
1944 if (vid == -1) error=1;
1945 #else
1946 if (pvr->longnm) {
1947 rc = nc_inq_varid(pfi->ncid, pvr->longnm, &vid);
1948 }
1949 else {
1950 rc = nc_inq_varid(pfi->ncid, pvr->abbrv, &vid);
1951 }
1952 if (rc != NC_NOERR) error=1;
1953 #endif
1954 if (error) {
1955 pvr->ncvid = -888; /* set flag so we won't try this variable ever again */
1956 if (pvr->longnm) {
1957 sprintf(pout,"Error: Variable %s not in netcdf file\n",pvr->longnm);
1958 }
1959 else {
1960 sprintf(pout,"Error: Variable %s not in netcdf file\n",pvr->abbrv);
1961 }
1962 gaprnt (0,pout);
1963 ncopts = oldncopts ;
1964 return (1);
1965 }
1966 /* No errors, so we can set the varid in the gavar structure */
1967 pvr->ncvid = vid;
1968
1969 /* If undef attribute name is given, get the undef value */
1970 if (pfi->undefattrflg) {
1971 #if USEHDF == 1
1972 /* get the undef attribute type */
1973 if (ncattinq(pfi->ncid, pvr->ncvid, pfi->undefattr, &attr_dtype, NULL) == -1) {
1974 sprintf(pout,"Warning: Could not retrieve data type for \"%s\" -- using %g instead\n",
1975 pfi->undefattr,pfi->undef);
1976 gaprnt(1,pout);
1977 pvr->undef = pfi->undef;
1978 }
1979 /* get the undef attribute value */
1980 switch(attr_dtype)
1981 {
1982 case (NC_SHORT):
1983 if (ncattget(pfi->ncid, pvr->ncvid, pfi->undefattr, &sundef) == -1) {
1984 sprintf(pout,"Warning: Could not retrieve \"%s\" -- using %g instead\n",
1985 pfi->undefattr,pfi->undef);
1986 gaprnt(1,pout);
1987 pvr->undef = pfi->undef;
1988 }
1989 else {
1990 pvr->undef = (float)sundef; /* Cast short to float */
1991 }
1992 break;
1993 case (NC_LONG):
1994 if (ncattget(pfi->ncid, pvr->ncvid, pfi->undefattr, &lundef) == -1) {
1995 sprintf(pout,"Warning: Could not retrieve \"%s\" -- using %g instead\n",
1996 pfi->undefattr,pfi->undef);
1997 gaprnt(1,pout);
1998 pvr->undef = pfi->undef;
1999 }
2000 else {
2001 pvr->undef = (float)lundef; /* Cast long to float */
2002 }
2003 break;
2004 case (NC_FLOAT):
2005 if (ncattget(pfi->ncid, pvr->ncvid, pfi->undefattr, &val) == -1) {
2006 sprintf(pout,"Warning: Could not retrieve \"%s\" -- using %g instead\n",
2007 pfi->undefattr,pfi->undef);
2008 gaprnt(1,pout);
2009 pvr->undef = pfi->undef;
2010 }
2011 else {
2012 pvr->undef = val; /* Value is already a float */
2013 }
2014 break;
2015 default:
2016 sprintf(pout,"Warning: data type for \"%s\" not handled -- using %g instead\n",
2017 pfi->undefattr,pfi->undef);
2018 gaprnt(1,pout);
2019 pvr->undef = pfi->undef;
2020 };
2021 #else
2022 if (nc_get_att_float(pfi->ncid, pvr->ncvid, pfi->undefattr, &val) != NC_NOERR) {
2023 sprintf(pout,"Warning: Could not retrieve \"%s\" -- using %g instead\n",
2024 pfi->undefattr,pfi->undef);
2025 gaprnt(1,pout);
2026 pvr->undef = pfi->undef;
2027 }
2028 else {
2029 pvr->undef = val;
2030 }
2031 #endif
2032 }
2033 /* If no undef attribute name is given, copy the file-wide undef */
2034 else {
2035 pvr->undef = pfi->undef;
2036 }
2037
2038 /* If data are packed, get the scale factor and offset attribute values */
2039 if (pfi->packflg) {
2040 /* initialize values */
2041 pvr->scale=1.0;
2042 pvr->add=0.0;
2043
2044 /* get scale factor */
2045 #if USEHDF == 1
2046 /* get the scale factor attribute type */
2047 if (ncattinq(pfi->ncid, pvr->ncvid, pfi->scattr, &attr_dtype, NULL) == -1) {
2048 gaprnt(1,"Warning: Could not retrieve scale factor data type -- setting scale factor to 1.0\n");
2049 pvr->scale = 1.0;
2050 }
2051 else {
2052 /* get the scale factor attribute value */
2053 switch(attr_dtype)
2054 {
2055 case (NC_SHORT):
2056 if (ncattget(pfi->ncid, pvr->ncvid, pfi->scattr, &sscale) == -1) {
2057 gaprnt(1,"Warning: Could not retrieve scale factor -- setting to 1.0\n");
2058 pvr->scale = 1.0;
2059 }
2060 else {
2061 pvr->scale = (float)sscale; /* Cast short to float */
2062 }
2063 break;
2064 case (NC_LONG):
2065 if (ncattget(pfi->ncid, pvr->ncvid, pfi->scattr, &lscale) == -1) {
2066 gaprnt(1,"Warning: Could not retrieve scale factor -- setting to 1.0\n");
2067 pvr->scale = 1.0;
2068 }
2069 else {
2070 pvr->scale = (float)lscale; /* Cast long to float */
2071 }
2072 break;
2073 case (NC_FLOAT):
2074 if (ncattget(pfi->ncid, pvr->ncvid, pfi->scattr, &val) == -1) {
2075 gaprnt(1,"Warning: Could not retrieve scale factor -- setting to 1.0\n");
2076 pvr->scale = 1.0;
2077 }
2078 else {
2079 pvr->scale = val; /* Value is already a float */
2080 }
2081 break;
2082 case (NC_DOUBLE):
2083 if (ncattget(pfi->ncid, pvr->ncvid, pfi->scattr, &dscale) == -1) {
2084 gaprnt(1,"Warning: Could not retrieve scale factor -- setting to 1.0\n");
2085 pvr->scale = 1.0;
2086 }
2087 else {
2088 pvr->scale = (float)dscale; /* Cast double to float */
2089 }
2090 break;
2091 default:
2092 gaprnt(1,"Warning: scale factor data type not handled -- setting scale factor to 1.0\n");
2093 pvr->scale = 1.0;
2094 };
2095 }
2096 #else
2097 /* get the scale factor attribute value */
2098 if (nc_get_att_float(pfi->ncid, pvr->ncvid, pfi->scattr, &val) != NC_NOERR) {
2099 gaprnt(1,"Warning: Could not retrieve scale factor -- setting to 1.0\n");
2100 pvr->scale = 1.0;
2101 }
2102 else {
2103 pvr->scale = val;
2104 }
2105 #endif
2106
2107 /* get add offset if required */
2108 if (pfi->packflg == 2) {
2109 #if USEHDF == 1
2110 /* get the add offset attribute type */
2111 if (ncattinq(pfi->ncid, pvr->ncvid, pfi->ofattr, &attr_dtype, NULL) == -1) {
2112 gaprnt(1,"Warning: Could not retrieve add offset data type -- setting add offset to 0.0\n");
2113 pvr->add = 0.0;
2114 }
2115 else {
2116 /* get the add offset attribute value */
2117 switch(attr_dtype)
2118 {
2119 case (NC_SHORT):
2120 if (ncattget(pfi->ncid, pvr->ncvid, pfi->ofattr, &sadd) == -1) {
2121 gaprnt(1,"Warning: Could not retrieve add offset -- setting to 0.0\n");
2122 pvr->add = 0.0;
2123 }
2124 else {
2125 pvr->add = (float)sadd; /* Cast short to float */
2126 }
2127 break;
2128 case (NC_LONG):
2129 if (ncattget(pfi->ncid, pvr->ncvid, pfi->ofattr, &ladd) == -1) {
2130 gaprnt(1,"Warning: Could not retrieve add offset -- setting to 0.0\n");
2131 pvr->add = 0.0;
2132 }
2133 else {
2134 pvr->add = (float)ladd; /* Cast long to float */
2135 }
2136 break;
2137 case (NC_FLOAT):
2138 if (ncattget(pfi->ncid, pvr->ncvid, pfi->ofattr, &val) == -1) {
2139 gaprnt(1,"Warning: Could not retrieve add offset -- setting to 0.0\n");
2140 pvr->add = 0.0;
2141 }
2142 else {
2143 pvr->add = val; /* Value is already a float */
2144 }
2145 break;
2146 case (NC_DOUBLE):
2147 if (ncattget(pfi->ncid, pvr->ncvid, pfi->ofattr, &dadd) == -1) {
2148 gaprnt(1,"Warning: Could not retrieve add offset -- setting to 0.0\n");
2149 pvr->add = 0.0;
2150 }
2151 else {
2152 pvr->add = (float)dadd; /* Cast double to float */
2153 }
2154 break;
2155 default:
2156 gaprnt(1,"Warning: add offset data type not handled -- setting add offset to 0.0\n");
2157 pvr->add = 0.0;
2158 };
2159 }
2160 #else
2161 /* get the add offset attribute value */
2162 if (nc_get_att_float(pfi->ncid, pvr->ncvid, pfi->ofattr, &val) != NC_NOERR) {
2163 gaprnt(1,"Warning: Could not retrieve add offset -- setting to 0.0\n");
2164 pvr->add = 0.0;
2165 }
2166 else {
2167 pvr->add = val;
2168 }
2169 #endif
2170 } /* matches if (pfi->packflg == 2) { */
2171 } /* matches if (pfi->packflg) { */
2172 }
2173
2174 /* Change the Y indexes if yrev flag is set */
2175 if (pfi->yrflg) {
2176 yy = pfi->dnum[1] - y;
2177 }
2178 else {
2179 yy = y-1;
2180 }
2181 /* Change the Z indexes if zrev flag is set */
2182 if (pfi->zrflg) {
2183 if (pvr->levels==0) {
2184 zz=0;
2185 }
2186 else {
2187 zz = pvr->levels-z;
2188 }
2189 }
2190 else {
2191 zz = z-1;
2192 }
2193
2194 /* Set up the start and count array. The units values
2195 provided for each variable indicate the mapping of the
2196 netcdf variable shape into the grads dimensions */
2197 for (i=0; i<16; i++) {
2198 start[i] = -999;
2199 count[i] = -999;
2200 if (pvr->units[i] == -100) { start[i] = x-1; count[i] = len; }
2201 if (pvr->units[i] == -101) { start[i] = yy; count[i] = 1; }
2202 if (pvr->units[i] == -102) { start[i] = zz; count[i] = 1; }
2203 if (pvr->units[i] == -103) { start[i] = t-1; count[i] = 1; }
2204 if (pvr->units[i] >=0) { start[i] = pvr->units[i]; count[i] = 1; }
2205 }
2206
2207 /* Now we are ready to do the I/O */
2208 #if USEHDF == 1
2209
2210 /* Get the data type */
2211 if (ncvarinq (pfi->ncid, pvr->ncvid, NULL, &data_dtype, NULL, NULL, NULL) == -1) {
2212 gaprnt(0,"NetCDF Error (ncvarinq): Unable to retrieve variable data type\n");
2213 ncopts = oldncopts ;
2214 return (1);
2215 }
2216 switch(data_dtype)
2217 {
2218 case (NC_SHORT):
2219 sval = (short *)malloc(len * sizeof(NC_SHORT));
2220 if (sval==NULL) {
2221 gaprnt(0,"Unable to allocate memory for data type NC_SHORT\n");
2222 ncopts = oldncopts ;
2223 return(1);
2224 }
2225 if (ncvarget (pfi->ncid, pvr->ncvid, (long *)start, (long *)count, (void*)sval) == -1) {
2226 gaprnt(0,"NetCDF Read Error (ncvarget) for data type NC_SHORT\n");
2227 ncopts = oldncopts ;
2228 return (1);
2229 }
2230 for (i=0; i<len; i++) gr[i] = (float)sval[i]; /* Cast short to float */
2231 free(sval);
2232 break;
2233
2234 case (NC_LONG):
2235 lval = (long *)malloc(len * sizeof(NC_LONG));
2236 if (lval==NULL) {
2237 gaprnt(0,"Unable to allocate memory for data type NC_LONG\n");
2238 ncopts = oldncopts ;
2239 return(1);
2240 }
2241 if (ncvarget (pfi->ncid, pvr->ncvid, (long *)start, (long *)count, (void*)lval) == -1) {
2242 gaprnt(0,"NetCDF Read Error (ncvarget) for data type NC_LONG\n");
2243 ncopts = oldncopts ;
2244 return (1);
2245 }
2246 for (i=0; i<len; i++) gr[i] = (float)lval[i]; /* Cast long to float */
2247 free(lval);
2248 break;
2249
2250 case (NC_FLOAT):
2251 if (ncvarget (pfi->ncid, pvr->ncvid, (long *)start, (long *)count, gr) == -1) {
2252 gaprnt(0,"NetCDF Read Error (ncvarget) for data type NC_FLOAT\n");
2253 ncopts = oldncopts ;
2254 return (1);
2255 }
2256 break;
2257
2258 default:
2259 sprintf(pout,"NetCDF Error: Data type %d not handled\n", data_dtype);
2260 gaprnt(0,pout);
2261 ncopts = oldncopts ;
2262 return(1);
2263 };
2264 #else
2265 rc = nc_get_vara_float(pfi->ncid, pvr->ncvid, start, count, gr);
2266 if (rc != NC_NOERR) {
2267 sprintf(pout,"NetCDF Error (nc_get_vara_float): %s\n",nc_strerror(rc));
2268 gaprnt(0,pout);
2269 ncopts = oldncopts ;
2270 return (1);
2271 }
2272 #endif
2273
2274 /* Set missing data values to gafile undef and then unpack if necessary */
2275 if (SETMISS) {
2276 /* use the gavar undef to set the fuzzy test limits */
2277 /* If gavar undef equals zero, change it to 1/EPSILON */
2278 if (pvr->undef == 0) {
2279 ulow = 1e-5;
2280 }
2281 else {
2282 ulow = fabs(pvr->undef/EPSILON);
2283 }
2284 uhi = pvr->undef + ulow;
2285 ulow = pvr->undef - ulow;
2286 /* set the gagrid undef equal to the gafile undef */
2287 pgr->undef = pfi->undef;
2288
2289 /* Do the fuzzy test for undef values before unpacking */
2290 for (i=0;i<len;i++) {
2291
2292 /* set the missing data to the gafile undef */
2293 if (*(gr+i) >= ulow && *(gr+i) <= uhi) {
2294 *(gr+i) = pfi->undef;
2295 }
2296 else {
2297 /* Data is not missing, so unpack with scale and offset if necessary */
2298 if (pfi->packflg) {
2299 *(gr+i) = *(gr+i)*pvr->scale + pvr->add;
2300 }
2301 }
2302 }
2303 }
2304
2305
2306 ncopts = oldncopts ;
2307 return(0);
2308 #else
2309 gaprnt(0,"Reading NetCDF files is not supported in this build\n");
2310 return(1);
2311 #endif
2312 }
2313
2314 /* Open an HDF-SDS file */
gaophdf(struct gafile * pfil,int tflag,int eflag)2315 int gaophdf (struct gafile *pfil, int tflag, int eflag) {
2316 #if USESDF == 1
2317 #if USEHDF == 1
2318 struct gavar *pvarl;
2319 int i,rc;
2320 int32 sd_id;
2321
2322 sd_id = -999;
2323 if (tflag) {
2324 sd_id = SDstart(pfil->tempname, DFACC_READ);
2325 }
2326 else {
2327 sd_id = SDstart(pfil->name, DFACC_READ);
2328 }
2329 if (sd_id==FAIL) {
2330 if (eflag) {
2331 gaprnt(0,"HDF-SDS Error: Unable to open file\n");
2332 }
2333 return (1);
2334 }
2335 pfil->sdid = sd_id;
2336 pvarl = pfil->pvar1;
2337 for (i=0; i<pfil->vnum; i++) {
2338 if (pvarl->sdvid != -888) pvarl->sdvid = -999;
2339 pvarl++;
2340 }
2341 return (0);
2342 #endif
2343 #endif
2344 gaprnt(0,"Opening HDF-SDS files is not supported in this build\n");
2345 return(1);
2346 }
2347
2348
2349 /* Read a row varying in the X direction from an HDF-SDS grid */
gahrow(int x,int y,int z,int t,int len,float * gr)2350 int gahrow (int x, int y, int z, int t, int len, float *gr) {
2351 #if USESDF == 1
2352 #if USEHDF == 1
2353 int rc,i,j,dim,size,yy,zz;
2354 int32 start[16],count[16];
2355 int32 sd_id, v_id, sds_id, attr_index;
2356 int32 data_dtype, n_atts, rank, dim_sizes[MAX_VAR_DIMS];
2357 int32 attr_dtype, attr_count;
2358 float val,ulow,uhi;
2359
2360 int8 *bval;
2361 uint8 *ubval;
2362 int16 *sval;
2363 uint16 *usval;
2364 int32 *ival;
2365 uint32 *uival;
2366 char *cattr_val;
2367 float32 *fattr_val;
2368 float64 *dattr_val;
2369
2370 /* Get the vid if we haven't already done that for this file */
2371 if (pvr->sdvid == -888) return(1); /* already tried and failed */
2372
2373 sd_id = pfi->sdid; /* this might be risky, casting int to int32 */
2374 if (pvr->sdvid == -999) {
2375
2376 /* Get the variable index number from the variable name */
2377 if (pvr->longnm) {
2378 v_id = SDnametoindex(sd_id, pvr->longnm);
2379 }
2380 else {
2381 v_id = SDnametoindex(sd_id, pvr->abbrv);
2382 }
2383 if (v_id==FAIL) {
2384 pvr->sdvid = -888;
2385 if (pvr->longnm) {
2386 sprintf(pout,"Error: Variable %s not in HDF-SDS file\n",pvr->longnm);
2387 }
2388 else {
2389 sprintf(pout,"Error: Variable %s not in HDF-SDS file\n",pvr->abbrv);
2390 }
2391 gaprnt(0,pout);
2392 return (1);
2393 }
2394 pvr->sdvid = v_id;
2395
2396 /* If undef attribute name is used, get the undef value */
2397 if (pfi->undefattrflg) {
2398 /* Select the variable (get sds_id) */
2399 v_id = pvr->sdvid;
2400 sds_id = SDselect(sd_id,v_id);
2401 if (sds_id==FAIL) {
2402 if (pvr->longnm) {
2403 sprintf(pout,"Error: SDselect failed for %s \n",pvr->longnm);
2404 }
2405 else {
2406 sprintf(pout,"Error: SDselect failed for %s \n",pvr->abbrv);
2407 }
2408 gaprnt(0,pout);
2409 return (1);
2410 }
2411 /* Retrieve the HDF undef attribute value */
2412 if (hdfattr(sds_id, pfi->undefattr, &val) != 0) {
2413 sprintf(pout,"Warning: Could not retrieve undef attribute \"%s\" -- using %g instead\n",
2414 pfi->undefattr,pfi->undef);
2415 gaprnt(1,pout);
2416 pvr->undef = pfi->undef;
2417 }
2418 else {
2419 pvr->undef = val;
2420 }
2421 }
2422 /* If undef attribute name is not given, copy the file-wide undef */
2423 else {
2424 pvr->undef = pfi->undef;
2425 }
2426
2427
2428 /* If data are packed, get the scale factor and offset attribute values */
2429 if (pfi->packflg) {
2430 /* initialize values */
2431 pvr->scale=1.0;
2432 pvr->add=0.0;
2433
2434 /* Select the variable (get sds_id) */
2435 v_id = pvr->sdvid;
2436 sds_id = SDselect(sd_id,v_id);
2437 if (sds_id==FAIL) {
2438 if (pvr->longnm) {
2439 sprintf(pout,"Error: SDselect failed for %s \n",pvr->longnm);
2440 }
2441 else {
2442 sprintf(pout,"Error: SDselect failed for %s \n",pvr->abbrv);
2443 }
2444 gaprnt(0,pout);
2445 return (1);
2446 }
2447 /* Retrieve the scale factor attribute value */
2448 if (hdfattr(sds_id, pfi->scattr, &pvr->scale) != 0) {
2449 sprintf(pout,"Warning: Could not retrieve \"%s\" -- setting to 1.0\n",pfi->scattr);
2450 gaprnt(1,pout);
2451 pvr->scale = 1.0;
2452 }
2453 /* Retrieve the add offset attribute value if required */
2454 if (pfi->packflg == 2) {
2455 if (hdfattr(sds_id, pfi->ofattr, &pvr->add) != 0) {
2456 sprintf(pout,"Warning: Could not retrieve \"%s\" -- setting to 0.0\n",pfi->ofattr);
2457 gaprnt(1,pout);
2458 pvr->add = 0.0;
2459 }
2460 }
2461 }
2462 }
2463
2464 /* Select the variable (get sds_id) */
2465 v_id = pvr->sdvid;
2466 sds_id = SDselect(sd_id,v_id);
2467
2468 if (sds_id==FAIL) {
2469 if (pvr->longnm) {
2470 sprintf (pout,"Error: SDselect failed for %s \n",pvr->longnm);
2471 }
2472 else {
2473 sprintf (pout,"Error: SDselect failed for %s \n",pvr->abbrv);
2474 }
2475 gaprnt(0,pout);
2476 return (1);
2477 }
2478
2479 /* Change the Y indexes if yrev flag is set */
2480 if (pfi->yrflg) {
2481 yy = pfi->dnum[1] - y;
2482 }
2483 else {
2484 yy = y-1;
2485 }
2486
2487 /* Change the Z indexes if zrev flag is set */
2488 if (pfi->zrflg) {
2489 if (pvr->levels==0) {
2490 zz=0;
2491 }
2492 else {
2493 zz = pvr->levels-z;
2494 }
2495 }
2496 else {
2497 zz = z-1;
2498 }
2499
2500 /* Set up the start and count array. The units records in the
2501 descriptor file for each variable indicate the mapping of the
2502 hdf-sds variable shape into the grads dimensions */
2503 for (i=0; i<16; i++) {
2504 start[i] = -999;
2505 count[i] = -999;
2506 if (pvr->units[i] == -100) { start[i] = x-1; count[i] = len;}
2507 if (pvr->units[i] == -101) { start[i] = yy; count[i] = 1; }
2508 if (pvr->units[i] == -102) { start[i] = zz; count[i] = 1; }
2509 if (pvr->units[i] == -103) { start[i] = t-1; count[i] = 1; }
2510 if (pvr->units[i] >= 0) { start[i] = pvr->units[i]; count[i] = 1; }
2511 }
2512
2513 /* Get the data type */
2514 if (pvr->longnm) {
2515 rc = SDgetinfo(sds_id, pvr->longnm, &rank, dim_sizes, &data_dtype, &n_atts);
2516 }
2517 else {
2518 rc = SDgetinfo(sds_id, pvr->abbrv, &rank, dim_sizes, &data_dtype, &n_atts);
2519 }
2520
2521 /* Data types that are handled are 8-bit unsigned ints (uint8), shorts (int16),
2522 ints (int32) and float. shorts and ints are converted to float.
2523 Unpacking done after I/O is finished */
2524 switch (data_dtype)
2525 {
2526 case (DFNT_INT8): /* definition value 20 */
2527 bval = (int8 *)malloc(len * sizeof (DFNT_INT8));
2528 if (bval==NULL) {
2529 gaprnt(0,"HDF-SDS Error: unable to allocate memory for dtype INT8\n");
2530 return(1);
2531 }
2532 if (SDreaddata(sds_id, start, NULL, count, (VOIDP *)bval) != 0) {
2533 gaprnt(0,"HDF-SDS Read Error for dtype INT8\n");
2534 return(1);
2535 }
2536 else {
2537 for (i=0; i<len; i++) gr[i] = (float)bval[i]; /* Cast int8 to float */
2538 }
2539 free(bval);
2540 break;
2541
2542 case (DFNT_UINT8): /* definition value 21 */
2543 ubval = (uint8 *)malloc(len * sizeof (DFNT_UINT8));
2544 if (ubval==NULL) {
2545 gaprnt(0,"HDF-SDS Error: unable to allocate memory for dtype UINT8\n");
2546 return(1);
2547 }
2548 if (SDreaddata(sds_id, start, NULL, count, (VOIDP *)ubval) != 0) {
2549 gaprnt(0,"HDF-SDS Read Error for dtype UINT8\n");
2550 return(1);
2551 }
2552 else {
2553 for (i=0; i<len; i++) gr[i] = (float)ubval[i]; /* Cast uint8 to float */
2554 }
2555 free(ubval);
2556 break;
2557
2558 case (DFNT_INT16): /* definition value 22 */
2559 sval = (int16 *)malloc(len * sizeof(DFNT_INT16));
2560 if (sval==NULL) {
2561 gaprnt(0,"HDF-SDS Error: unable to allocate memory for dtype INT16\n");
2562 return(1);
2563 }
2564 if (SDreaddata(sds_id, start, NULL, count, (VOIDP *)sval) != 0) {
2565 gaprnt(0,"HDF-SDS Read Error for dtype INT16\n");
2566 return(1);
2567 }
2568 else {
2569 for (i=0; i<len; i++) gr[i] = (float)sval[i]; /* Cast int16 to float */
2570 }
2571 free(sval);
2572 break;
2573
2574 case (DFNT_UINT16): /* definition value 23 */
2575 usval = (uint16 *)malloc(len * sizeof (DFNT_UINT16));
2576 if (usval==NULL) {
2577 gaprnt(0,"HDF-SDS Error: unable to allocate memory for dtype UINT16\n");
2578 return(1);
2579 }
2580 if (SDreaddata(sds_id, start, NULL, count, (VOIDP *)usval) != 0) {
2581 gaprnt(0,"HDF-SDS Read Error for dtype UINT16\n");
2582 return(1);
2583 }
2584 else {
2585 for (i=0; i<len; i++) gr[i] = (float)usval[i]; /* Cast uint16 to float */
2586 }
2587 free(usval);
2588 break;
2589
2590 case (DFNT_INT32): /* definition value 24 */
2591 ival = (int32 *)malloc(len * sizeof (DFNT_INT32));
2592 if (ival==NULL) {
2593 gaprnt(0,"HDF-SDS Error: unable to allocate memory for dtype INT32\n");
2594 return(1);
2595 }
2596 if (SDreaddata(sds_id, start, NULL, count, (VOIDP *)ival) != 0) {
2597 gaprnt(0,"HDF-SDS Read Error for dtype INT32\n");
2598 return(1);
2599 }
2600 else {
2601 for (i=0; i<len; i++) gr[i] = (float)ival[i]; /* Cast int32 to float */
2602 }
2603 free(ival);
2604 break;
2605
2606 case (DFNT_UINT32): /* definition value 25 */
2607 uival = (uint32 *)malloc(len * sizeof (DFNT_UINT32));
2608 if (uival==NULL) {
2609 gaprnt(0,"HDF-SDS Error: unable to allocate memory for dtype UINT32\n");
2610 return(1);
2611 }
2612 if (SDreaddata(sds_id, start, NULL, count, (VOIDP *)uival) != 0) {
2613 gaprnt(0,"HDF-SDS Read Error for dtype UINT32\n");
2614 return(1);
2615 }
2616 else {
2617 for (i=0; i<len; i++) gr[i] = (float)uival[i]; /* Cast uint32 to float */
2618 }
2619 free(uival);
2620 break;
2621
2622 case (DFNT_FLOAT32): /* definition value 5 */
2623 if (SDreaddata(sds_id, start, NULL, count, (VOIDP *)gr) != 0) {
2624 gaprnt(0,"HDF-SDS Read Error for dtype FLOAT32\n");
2625 return(1);
2626 }
2627 break;
2628
2629 default:
2630 sprintf(pout,"HDF-SDS Error: Data type %d not handled\n", data_dtype);
2631 gaprnt(0,pout);
2632 return(1);
2633 };
2634
2635 /* Set missing data values to exact value if specified */
2636 if (SETMISS) {
2637 /* Use the gavar undef to set the fuzzy test limits */
2638 /* If gavar undef equals zero, change it to 1/EPSILON */
2639 if (pvr->undef == 0) {
2640 ulow = 1e-5;
2641 }
2642 else {
2643 ulow = fabs(pvr->undef/EPSILON);
2644 }
2645 uhi = pvr->undef + ulow;
2646 ulow = pvr->undef - ulow;
2647 /* set the gagrid undef equal to the gafile undef */
2648 pgr->undef = pfi->undef;
2649
2650 /* Do the fuzzy test for undef values before unpacking */
2651 for (i=0;i<len;i++) {
2652 /* set the missing data to the gafile undef */
2653 if (*(gr+i) >= ulow && *(gr+i) <= uhi) {
2654 *(gr+i) = pfi->undef;
2655 }
2656 else {
2657 /* Data is not missing, so unpack with scale and offset if necessary */
2658 if (pfi->packflg) {
2659 *(gr+i) = *(gr+i)*pvr->scale + pvr->add;
2660 }
2661 }
2662 }
2663 }
2664
2665 return (0);
2666
2667 #endif
2668 #endif
2669 gaprnt(0,"Reading HDF-SDS files is not supported in this build\n");
2670 return(1);
2671 }
2672
2673
2674 /* Retrieves a non-character HDF-SDS Attribute.
2675 Attribute values that are short, int, or double are converted into floats. */
hdfattr(int sds_id,char * attr_name,float * value)2676 int hdfattr(int sds_id, char *attr_name, float *value) {
2677 #if USESDF == 1
2678 #if USEHDF == 1
2679 int32 rc, attr_index, attr_dtype, attr_count;
2680 uint8 *battr_val;
2681 int16 *sattr_val;
2682 uint16 *usattr_val;
2683 int32 *iattr_val;
2684 uint32 *uiattr_val;
2685 float32 *fattr_val;
2686 float64 *dattr_val;
2687
2688 /* Get the attribute index number from its name */
2689 attr_index = SDfindattr(sds_id, attr_name);
2690 if (attr_index == -1) {
2691 sprintf(pout,"Warning: HDF attribute named \"%s\" does not exist\n",attr_name);
2692 gaprnt(1,pout);
2693 return(1);
2694 }
2695
2696 /* Get info about the attribute, make sure there's only one value */
2697 if (SDattrinfo(sds_id, attr_index, attr_name, &attr_dtype, &attr_count) == -1) {
2698 gaprnt(1,"Warning: SDattrinfo failed\n");
2699 return(1);
2700 }
2701 else {
2702 if (attr_count != 1) {
2703 sprintf(pout,"Warning: HDF attribute named \"%s\" has more than one value\n",attr_name);
2704 gaprnt(1,pout);
2705 return(1);
2706 }
2707
2708 /* Read the attribute value */
2709 switch (attr_dtype)
2710 {
2711 case (DFNT_UINT8): /* definition value 21 */
2712 battr_val = malloc (attr_count * sizeof (DFNT_UINT8));
2713 if (SDreadattr(sds_id, attr_index, battr_val) == -1) {
2714 gaprnt(1,"Warning: SDreadattr failed for attribute type UINT8\n");
2715 return(1);
2716 }
2717 else {
2718 *value = *battr_val;
2719 }
2720 free(battr_val);
2721 break;
2722
2723 case (DFNT_INT16): /* definition value 22 */
2724 sattr_val = malloc (attr_count * sizeof (DFNT_INT16));
2725 if (SDreadattr(sds_id, attr_index, sattr_val) == -1) {
2726 gaprnt(1,"Warning: SDreadattr failed for attribute type INT16\n");
2727 return(1);
2728 }
2729 else {
2730 *value = *sattr_val;
2731 }
2732 free(sattr_val);
2733 break;
2734
2735 case (DFNT_UINT16): /* definition value 23 */
2736 usattr_val = malloc (attr_count * sizeof (DFNT_UINT16));
2737 if (SDreadattr(sds_id, attr_index, usattr_val) == -1) {
2738 gaprnt(1,"Warning: SDreadattr failed for attribute type UINT16\n");
2739 return(1);
2740 }
2741 else {
2742 *value = *usattr_val;
2743 }
2744 free(usattr_val);
2745 break;
2746
2747 case (DFNT_INT32): /* definition value 24 */
2748 iattr_val = malloc (attr_count * sizeof (DFNT_INT32));
2749 if (SDreadattr(sds_id, attr_index, iattr_val) == -1) {
2750 gaprnt(1,"Warning: SDreadattr failed for attribute type INT32\n");
2751 return(1);
2752 }
2753 else {
2754 *value = *iattr_val;
2755 }
2756 free(iattr_val);
2757 break;
2758
2759 case (DFNT_UINT32): /* definition value 25 */
2760 uiattr_val = malloc (attr_count * sizeof (DFNT_UINT32));
2761 if (SDreadattr(sds_id, attr_index, uiattr_val) == -1) {
2762 gaprnt(1,"Warning: SDreadattr failed for attribute type UINT32\n");
2763 return(1);
2764 }
2765 else {
2766 *value = *uiattr_val;
2767 }
2768 free(uiattr_val);
2769 break;
2770
2771 case (DFNT_FLOAT32): /* definition value 5 */
2772 fattr_val = malloc (attr_count * sizeof (DFNT_FLOAT32));
2773 if (SDreadattr(sds_id, attr_index, fattr_val) == -1) {
2774 gaprnt(1,"Warning: SDreadattr failed for attribute type FLOAT32\n");
2775 return(1);
2776 }
2777 else {
2778 *value = *fattr_val;
2779 }
2780 free(fattr_val);
2781 break;
2782
2783 case (DFNT_FLOAT64): /* definition value 6 */
2784 dattr_val = malloc (attr_count * sizeof (DFNT_FLOAT64));
2785 if (SDreadattr(sds_id, attr_index, dattr_val) == -1) {
2786 gaprnt(1,"Warning: SDreadattr failed for attribute type FLOAT64\n");
2787 return(1);
2788 }
2789 else {
2790 *value = *dattr_val;
2791 }
2792 free(dattr_val);
2793 break;
2794
2795 default:
2796 sprintf(pout,"Warning: HDF Attribute \"%s\" is not a numeric data type (%d)\n",
2797 attr_name, attr_dtype);
2798 gaprnt(1,pout);
2799 return(1);
2800 };
2801 }
2802 return(0);
2803
2804 #endif
2805 #endif
2806 gaprnt(0,"Reading HDF-SDS files is not supported in this build\n");
2807 return(1);
2808 }
2809
2810
2811 /* Subroutine to print out NetCDF attributes */
ncpattrs(int ncid,char * varnam,char * abbrv,int hdrflg,int fnum,char * ftit)2812 int ncpattrs(int ncid, char *varnam, char *abbrv, int hdrflg, int fnum, char* ftit) {
2813 #if USESDF == 1
2814 int attr_index, attr_count, rc, i, truncflag, varid, n_atts;
2815 char attr_name[MAX_NC_NAME];
2816 nc_type attr_dtype;
2817 char *battr_val;
2818 char *cattr_val;
2819 short *sattr_val;
2820 long *iattr_val;
2821 float *fattr_val;
2822 double *dattr_val;
2823 int error=0;
2824
2825 /* Get the variable id and number of attributes */
2826 #if USEHDF == 1
2827 if (cmpwrd("global",abbrv)) {
2828 varid = NC_GLOBAL;
2829 rc = ncinquire(ncid, NULL, NULL, &n_atts, NULL);
2830 if (rc == -1) error=1;
2831 }
2832 else {
2833 varid = ncvarid(ncid, varnam);
2834 if (varid == -1) error=1;
2835 if (!error) {
2836 rc = ncvarinq(ncid, varid, NULL, NULL, NULL, NULL, &n_atts);
2837 if (rc == -1) error=1;
2838 }
2839 }
2840 #else
2841 if (cmpwrd("global",abbrv)) {
2842 varid = NC_GLOBAL;
2843 rc = nc_inq_natts(ncid,&n_atts);
2844 if (rc != NC_NOERR) error=1;
2845 }
2846 else {
2847 rc = nc_inq_varid(ncid, varnam, &varid);
2848 if (rc != NC_NOERR) error=1;
2849 if (!error) {
2850 rc = nc_inq_varnatts(ncid, varid, &n_atts);
2851 if (rc != NC_NOERR) error=1;
2852 }
2853 }
2854 #endif
2855
2856 /* Print out the header */
2857 if (!error) {
2858 if (hdrflg) {
2859 if (n_atts > 0) {
2860 sprintf(pout,"Native Attributes for File %i : %s \n",fnum,ftit);
2861 gaprnt(2,pout);
2862 }
2863 }
2864 }
2865 else {
2866 return(0); /* zero attributes printed */
2867 }
2868
2869 /* Loop through list of attributes, print the name of each one */
2870 for (attr_index = 0 ; attr_index < n_atts ; attr_index++) {
2871
2872 /* Get current attribute's name */
2873 #if USEHDF == 1
2874 if (ncattname(ncid, varid, attr_index, attr_name) == -1) {
2875 sprintf(pout,"ncattname failed for variable %s, attribute number %d\n", abbrv, attr_index);
2876 gaprnt(2,pout);
2877 }
2878 #else
2879 if (nc_inq_attname(ncid, varid, attr_index, attr_name) == -1) {
2880 sprintf(pout,"nc_inq_attname failed for variable %s, attribute number %d\n", abbrv, attr_index);
2881 gaprnt(2,pout);
2882 }
2883 #endif
2884 else {
2885 /* Get current attribute's data type and length */
2886 #if USEHDF == 1
2887 if (ncattinq(ncid, varid, attr_name, &attr_dtype, &attr_count) == -1){
2888 sprintf(pout,"ncattinq failed for variable %s, attribute number %d\n", abbrv, attr_index);
2889 gaprnt(2,pout);
2890 }
2891 #else
2892 if (nc_inq_att(ncid, varid, attr_name, &attr_dtype, (size_t*)&attr_count) == -1) {
2893 sprintf(pout,"nc_inq_att failed for variable %s, attribute number %d\n", abbrv, attr_index);
2894 gaprnt(2,pout);
2895 }
2896 #endif
2897 else {
2898 if (attr_count>0) {
2899 /* Retrieve and print out the attribute */
2900 switch (attr_dtype)
2901 {
2902 case (NC_BYTE):
2903 battr_val = malloc ((attr_count+1) * sizeof (NC_BYTE));
2904 #if USEHDF == 1
2905 if (ncattget(ncid, varid, attr_name, battr_val) == -1) {
2906 gaprnt(2,"ncattget failed for type NC_BYTE\n");
2907 }
2908 #else
2909 if (nc_get_att_schar(ncid, varid, attr_name, (signed char*)battr_val) == -1) {
2910 gaprnt(2,"nc_get_att_schar failed for type NC_BYTE\n");
2911 }
2912 #endif
2913 else {
2914 gaprnt(2,abbrv);
2915 gaprnt(2," Byte ");
2916 gaprnt(2,attr_name);
2917 gaprnt(2," ");
2918 for (i=0; i<attr_count; i++) {
2919 sprintf(pout,"%d ", (int)(battr_val[i]));
2920 gaprnt(2,pout);
2921 }
2922 gaprnt(2,"\n");
2923 }
2924 free(battr_val);
2925 break;
2926 case (NC_CHAR):
2927 cattr_val = malloc ((attr_count+1) * sizeof (NC_CHAR));
2928 #if USEHDF == 1
2929 if (ncattget(ncid, varid, attr_name, cattr_val) == -1) {
2930 gaprnt(2,"ncattget failed for type NC_CHAR\n");
2931 }
2932 #else
2933 if (nc_get_att_text(ncid, varid, attr_name, cattr_val) == -1) {
2934 gaprnt(2,"nc_get_att_text failed for type NC_CHAR\n");
2935 }
2936 #endif
2937 else {
2938 cattr_val[attr_count]='\0';
2939 gaprnt(2,abbrv);
2940 gaprnt(2," String ");
2941 gaprnt(2,attr_name);
2942 gaprnt(2," ");
2943 prntwrap(abbrv, attr_name, cattr_val);
2944 }
2945 free(cattr_val);
2946 break;
2947 case (NC_SHORT):
2948 sattr_val = malloc (attr_count * sizeof (NC_SHORT));
2949 #if USEHDF == 1
2950 if (ncattget(ncid, varid, attr_name, sattr_val) == -1) {
2951 gaprnt(2,"ncattget failed for type NC_SHORT\n");
2952 }
2953 #else
2954 if (nc_get_att_short(ncid, varid, attr_name, sattr_val) == -1) {
2955 gaprnt(2,"nc_get_att_short failed for type NC_SHORT\n");
2956 }
2957 #endif
2958 else {
2959 gaprnt(2,abbrv);
2960 gaprnt(2," Int16 ");
2961 gaprnt(2,attr_name);
2962 gaprnt(2," ");
2963 for (i=0; i<attr_count; i++) {
2964 sprintf(pout,"%d ", (int)(sattr_val[i]));
2965 gaprnt(2,pout);
2966 }
2967 gaprnt(2,"\n");
2968 }
2969 free(sattr_val);
2970 break;
2971 case (NC_LONG):
2972 iattr_val = malloc (attr_count * sizeof (NC_LONG));
2973 #if USEHDF == 1
2974 if (ncattget(ncid, varid, attr_name, iattr_val) == -1) {
2975 gaprnt(2,"ncattget failed for type NC_LONG\n");
2976 }
2977 #else
2978 if (nc_get_att_long(ncid, varid, attr_name, iattr_val) == -1) {
2979 gaprnt(2,"nc_get_att_long failed for type NC_LONG\n");
2980 }
2981 #endif
2982 else {
2983 gaprnt(2,abbrv);
2984 gaprnt(2," Int32 ");
2985 gaprnt(2,attr_name);
2986 gaprnt(2," ");
2987 for (i=0; i<attr_count; i++) {
2988 sprintf(pout,"%d ", iattr_val[i]);
2989 gaprnt(2,pout);
2990 }
2991 gaprnt(2,"\n");
2992 }
2993 free(iattr_val);
2994 break;
2995 case (NC_FLOAT):
2996 fattr_val = malloc (attr_count * sizeof (NC_FLOAT));
2997 #if USEHDF == 1
2998 if (ncattget(ncid, varid, attr_name, fattr_val) == -1) {
2999 gaprnt(2,"ncattget failed for type NC_FLOAT\n");
3000 }
3001 #else
3002 if (nc_get_att_float(ncid, varid, attr_name, fattr_val) == -1) {
3003 gaprnt(2,"nc_get_att_float failed for type NC_FLOAT\n");
3004 }
3005 #endif
3006 else {
3007 gaprnt(2,abbrv);
3008 gaprnt(2," Float32 ");
3009 gaprnt(2,attr_name);
3010 gaprnt(2," ");
3011 for (i=0; i<attr_count; i++) {
3012 sprintf(pout,"%g ", fattr_val[i]);
3013 gaprnt(2,pout);
3014 }
3015 gaprnt(2,"\n");
3016 }
3017 free(fattr_val);
3018 break;
3019 case (NC_DOUBLE):
3020 dattr_val = malloc (attr_count * sizeof (NC_DOUBLE));
3021 #if USEHDF == 1
3022 if (ncattget(ncid, varid, attr_name, dattr_val) == -1) {
3023 gaprnt(2,"ncattget failed for type NC_DOUBLE\n");
3024 }
3025 #else
3026 if (nc_get_att_double(ncid, varid, attr_name, dattr_val) == -1) {
3027 gaprnt(2,"nc_get_att_double failed for type NC_FLOAT\n");
3028 }
3029 #endif
3030 else {
3031 gaprnt(2,abbrv);
3032 gaprnt(2," Float64 ");
3033 gaprnt(2,attr_name);
3034 gaprnt(2," ");
3035 for (i=0; i<attr_count; i++) {
3036 sprintf(pout,"%g ", dattr_val[i]);
3037 gaprnt(2,pout);
3038 }
3039 gaprnt(2,"\n");
3040 }
3041 free(dattr_val);
3042 break;
3043 default:
3044 sprintf(pout,"Failed to retrieve attribute %d of type %d \n", attr_index, attr_dtype);
3045 gaprnt(2,pout);
3046 };
3047 } /* end of if statement for attr_count >0 */
3048 } /* end of if-else statement for ncattinq */
3049 } /* end of if-else statement for ncattname */
3050 } /* end of for loop on attr_index */
3051 return(n_atts);
3052 #endif
3053 gaprnt(0,"Reading NetCDF attributes is not supported in this build\n");
3054 return(1);
3055 }
3056
3057
3058 /* Subroutine to print out HDF attributes */
hdfpattrs(int sdid,char * varname,char * abbrv,int hdrflg,int fnum,char * ftit)3059 int hdfpattrs(int sdid, char *varname, char *abbrv, int hdrflg, int fnum, char* ftit) {
3060 #if USESDF == 1
3061 #if USEHDF == 1
3062 int attr_index, rc, i;
3063 char *pos, *line;
3064 char attr_name[MAX_NC_NAME];
3065 int32 attr_dtype, attr_count;
3066 char8 *cattr_val;
3067 uchar8 *ucattr_val;
3068 int8 *icattr_val;
3069 uint8 *uicattr_val;
3070 int16 *sattr_val;
3071 uint16 *usattr_val;
3072 int32 *iattr_val;
3073 uint32 *uiattr_val;
3074 float32 *fattr_val;
3075 float64 *dattr_val;
3076 int error=0;
3077 char name[MAX_NC_NAME];
3078 int32 sds_id, n_atts, n_dsets, varid, rank, type, dim_sizes[4];
3079
3080
3081 /* Get the dataset id and number of attributes */
3082 if (cmpwrd("global",abbrv)) {
3083 sds_id = sdid;
3084 rc = SDfileinfo(sdid, &n_dsets, &n_atts);
3085 if (rc == -1) error=1;
3086 }
3087 else {
3088 sds_id = SDnametoindex(sdid, varname);
3089 if (sds_id == -1) error=1;
3090 if (!error) {
3091 sds_id = SDselect(sdid,sds_id);
3092 rc = SDgetinfo(sds_id, name, &rank, dim_sizes, &type, &n_atts);
3093 if (rc == -1) error=1;
3094 }
3095 }
3096 /* Print out the header */
3097 if (!error) {
3098 if (hdrflg) {
3099 if (n_atts > 0) {
3100 sprintf(pout,"Native Attributes for File %i : %s \n",fnum,ftit);
3101 gaprnt(2,pout);
3102 }
3103 }
3104 }
3105 else {
3106 return(0); /* zero attributes printed */
3107 }
3108
3109 /* Loop through list of attributes, print the name of each one */
3110 for (attr_index = 0 ; attr_index < n_atts ; attr_index++) {
3111
3112 /* Get info about the current attribute and then print out Name, Type, and Value */
3113 if (SDattrinfo(sds_id, attr_index, attr_name, &attr_dtype, &attr_count) == -1) {
3114 sprintf(pout,"SDattrinfo failed for variable %s, attribute number %d\n", abbrv, attr_index);
3115 gaprnt(2,pout);
3116 }
3117 else {
3118 switch (attr_dtype)
3119 {
3120 case (DFNT_CHAR8): /* definition value 4 */
3121 cattr_val = malloc ((attr_count+1) * sizeof (char8));
3122 if (SDreadattr(sds_id, attr_index, cattr_val) == -1) {
3123 gaprnt(2,"SDreadattr failed for type CHAR8\n");
3124 }
3125 else {
3126 cattr_val[attr_count]='\0';
3127 gaprnt(2,abbrv);
3128 gaprnt(2," String ");
3129 gaprnt(2,attr_name);
3130 gaprnt(2," ");
3131 prntwrap(abbrv, attr_name, cattr_val);
3132 }
3133 free(cattr_val);
3134 break;
3135 case (DFNT_UCHAR8): /* definition value 3 */
3136 ucattr_val = malloc ((attr_count+1) * sizeof (uchar8));
3137 if (SDreadattr(sds_id, attr_index, ucattr_val) == -1) {
3138 gaprnt(2,"SDreadattr failed for type UCHAR8\n");
3139 }
3140 else {
3141 ucattr_val[attr_count]='\0';
3142 gaprnt(2,abbrv);
3143 gaprnt(2," String ");
3144 gaprnt(2,attr_name);
3145 gaprnt(2," ");
3146 prntwrap(abbrv, attr_name, (char*)ucattr_val);
3147 }
3148 free(ucattr_val);
3149 break;
3150 case (DFNT_INT8): /* definition value 20 */
3151 icattr_val = malloc ((attr_count+1) * sizeof (int8));
3152 if (SDreadattr(sds_id, attr_index, icattr_val) == -1) {
3153 gaprnt(2,"SDreadattr failed for type INT8\n");
3154 }
3155 else {
3156 gaprnt(2,abbrv);
3157 gaprnt(2," Byte ");
3158 gaprnt(2,attr_name);
3159 gaprnt(2," ");
3160 for (i=0; i<attr_count; i++) {
3161 sprintf(pout,"%d ", (int)(icattr_val[i]));
3162 gaprnt(2,pout);
3163 }
3164 gaprnt(2,"\n");
3165 }
3166 free(cattr_val);
3167 break;
3168 case (DFNT_UINT8): /* definition value 21 */
3169 uicattr_val = malloc ((attr_count) * sizeof (uint8));
3170 if (SDreadattr(sds_id, attr_index, cattr_val) == -1) {
3171 gaprnt(2,"SDreadattr failed for type UINT8\n");
3172 }
3173 else {
3174 gaprnt(2,abbrv);
3175 gaprnt(2," Byte ");
3176 gaprnt(2,attr_name);
3177 gaprnt(2," ");
3178 for (i=0; i<attr_count; i++) {
3179 sprintf(pout,"%u ", (unsigned int)(uicattr_val[i]));
3180 gaprnt(2,pout);
3181 }
3182 gaprnt(2,"\n");
3183 }
3184 free(uicattr_val);
3185 break;
3186 case (DFNT_INT16): /* definition value 22 */
3187 sattr_val = malloc (attr_count * sizeof (int16));
3188 if (SDreadattr(sds_id, attr_index, sattr_val) == -1) {
3189 gaprnt(2,"SDreadattr failed for type INT16\n");
3190 }
3191 else {
3192 gaprnt(2,abbrv);
3193 gaprnt(2," Int16 ");
3194 gaprnt(2,attr_name);
3195 gaprnt(2," ");
3196 for (i=0; i<attr_count; i++) {
3197 sprintf(pout,"%d ", (int)(sattr_val[i]));
3198 gaprnt(2,pout);
3199 }
3200 gaprnt(2,"\n");
3201 }
3202 free(sattr_val);
3203 break;
3204 case (DFNT_UINT16): /* definition value 23 */
3205 usattr_val = malloc (attr_count * sizeof (uint16));
3206 if (SDreadattr(sds_id, attr_index, usattr_val) == -1) {
3207 gaprnt(2,"SDreadattr failed for type UINT16\n");
3208 }
3209 else {
3210 gaprnt(2,abbrv);
3211 gaprnt(2," UInt16 ");
3212 gaprnt(2,attr_name);
3213 gaprnt(2," ");
3214 for (i=0; i<attr_count; i++) {
3215 sprintf(pout,"%u ", (unsigned int)(usattr_val[i]));
3216 gaprnt(2,pout);
3217 }
3218 gaprnt(2,"\n");
3219 }
3220 free(usattr_val);
3221 break;
3222 case (DFNT_INT32): /* definition value 24 */
3223 iattr_val = malloc (attr_count * sizeof (int32));
3224 if (SDreadattr(sds_id, attr_index, iattr_val) == -1) {
3225 gaprnt(2,"SDreadattr failed for type INT32\n");
3226 }
3227 else {
3228 gaprnt(2,abbrv);
3229 gaprnt(2," Int32 ");
3230 gaprnt(2,attr_name);
3231 gaprnt(2," ");
3232 for (i=0; i<attr_count; i++) {
3233 sprintf(pout,"%d ", iattr_val[i]);
3234 gaprnt(2,pout);
3235 }
3236 gaprnt(2,"\n");
3237 }
3238 free(iattr_val);
3239 break;
3240 case (DFNT_UINT32): /* definition value 25 */
3241 uiattr_val = malloc (attr_count * sizeof (uint32));
3242 if (SDreadattr(sds_id, attr_index, uiattr_val) == -1) {
3243 gaprnt(2,"SDreadattr failed for type UINT32\n");
3244 }
3245 else {
3246 gaprnt(2,abbrv);
3247 gaprnt(2," UInt32 ");
3248 gaprnt(2,attr_name);
3249 gaprnt(2," ");
3250 for (i=0; i<attr_count; i++) {
3251 sprintf(pout,"%u ", uiattr_val[i]);
3252 gaprnt(2,pout);
3253 }
3254 gaprnt(2,"\n");
3255 }
3256 free(iattr_val);
3257 break;
3258 case (DFNT_FLOAT32): /* definition value 5 */
3259 fattr_val = malloc (attr_count * sizeof (float32));
3260 if (SDreadattr(sds_id, attr_index, fattr_val) == -1) {
3261 gaprnt(2,"SDreadattr failed for type FLOAT32\n");
3262 }
3263 else {
3264 gaprnt(2,abbrv);
3265 gaprnt(2," Float32 ");
3266 gaprnt(2,attr_name);
3267 gaprnt(2," ");
3268 for (i=0; i<attr_count; i++) {
3269 sprintf(pout,"%g ", fattr_val[i]);
3270 gaprnt(2,pout);
3271 }
3272 gaprnt(2,"\n");
3273 }
3274 free(fattr_val);
3275 break;
3276 case (DFNT_FLOAT64): /* definition value 6 */
3277 dattr_val = malloc (attr_count * sizeof (float64));
3278 if (SDreadattr(sds_id, attr_index, dattr_val) == -1) {
3279 gaprnt(2,"SDreadattr failed for type FLOAT64\n");
3280 }
3281 else {
3282 gaprnt(2,abbrv);
3283 gaprnt(2," Float64 ");
3284 gaprnt(2,attr_name);
3285 gaprnt(2," ");
3286 for (i=0; i<attr_count; i++) {
3287 sprintf(pout,"%g ", dattr_val[i]);
3288 gaprnt(2,pout);
3289 }
3290 gaprnt(2,"\n");
3291 }
3292 free(dattr_val);
3293 break;
3294 default:
3295 sprintf(pout,"Failed to retrieve attribute %d of type %d \n", attr_index, attr_dtype);
3296 gaprnt(2,pout);
3297 };
3298 } /* end of if-else statment following call to SDattrinfo */
3299 }
3300 return(n_atts);
3301 #endif
3302 #endif
3303 gaprnt(0,"Reading HDF-SDS attributes is not supported in this build\n");
3304 return(1);
3305 }
3306
3307 /* routine to print out a string attribute that may have carriage returns in it */
prntwrap(char * vname,char * aname,char * str)3308 void prntwrap(char *vname, char *aname, char *str ) {
3309 char *pos, *line;
3310 pos = line = str;
3311 while (*pos != '\0') {
3312 if (*pos == '\n') {
3313 *pos = '\0'; /* swap null for carriage return */
3314 gaprnt(2,line);
3315 sprintf(pout," \n%s String %s ",vname,aname); /* add varname, attr_type, and attr_name after carriage return */
3316 gaprnt(2,pout);
3317 *pos = '\n'; /* put the carriage return back in */
3318 line = pos+1;
3319 }
3320 pos++;
3321 }
3322 if (line < pos) { /* Print string that has no carriage returns in it */
3323 gaprnt(2,line);
3324 }
3325 gaprnt(2,"\n");
3326 }
3327