3 This is the TRMM Office Radar Software Library.
4 Copyright (C) 1996, 1997
6 Space Applications Corporation
9 This library is free software; you can redistribute it and/or
10 modify it under the terms of the GNU Library General Public
11 License as published by the Free Software Foundation; either
12 version 2 of the License, or (at your option) any later version.
14 This library is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 Library General Public License for more details.
19 You should have received a copy of the GNU Library General Public
20 License along with this library; if not, write to the Free
21 Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
24 * Read SIGMET version 1 and version 2 formatted files.
26 * Data is written in little-endian format for version 1 files.
27 * This means that on big endian machines, bytes must be swapped.
28 * This is auto-detected and all swapping is automatic.
30 * Note that this is different in SIGMET version 2 data files. There
31 * the data is written in big-endian format (written w/ an IRIS computer).
32 * For that case, the byte swapping logic is reversed.
34 * The highest level functions provided is:
35 * nsig_read_sweep -- Read an entire sweep including all fields.
36 * Call it until NULL is returned. That indicates
38 *----------------------------------------------------------------------
42 * Space Applications Corp.
43 * NASA/GSFC Code 910.1
55 FILE *uncompress_pipe(FILE *fp);
57 int little_endian(void);
58 void swap_4_bytes(void *word);
59 void swap_2_bytes(void *word);
60 int rsl_pclose(FILE *fp);
62 /*********************************************************************
63 * Open a file and possibly setup a gunzip pipe *
64 *********************************************************************/
65 FILE *nsig_open(char *file_name)
71 if (file_name == NULL) { /* Use stdin */
73 fp = fdopen(save_fd, "r");
74 } else if((fp = fopen(file_name,"r")) == NULL) {
79 fp = uncompress_pipe(fp); /* Transparently gunzip. */
83 /**********************************************************************
84 * Given an opened file stream read in the headers and fill in *
85 * the nsig_file data structure. *
86 **********************************************************************/
87 int nsig_read_record(FILE *fp, char *nsig_rec)
92 /* Input could be a chain of pipes. So read until no more data.
93 * For instance, a gzip pipe will write chunks of 4096 bytes,
94 * even when there is really more that we want to read.
96 /* Return the number of bytes read. */
97 buf = (char *)nsig_rec;
99 if (feof(fp)) return -1;
101 while((n = fread(&buf[nbytes], sizeof(char), NSIG_BLOCK-nbytes, fp)) > 0) {
108 /*********************************************************
110 *********************************************************/
111 void nsig_close(FILE *fp)
118 int nsig_endianess(NSIG_Record1 *rec1)
121 * If NSIG is version1 and on big-endian, then swap.
122 * If NSIG is version2 and on little-endian, then swap.
125 printf("id = %d %d\n", (int)rec1->struct_head.id[0], (int)rec1->struct_head.id[1]);
127 if (rec1->struct_head.id[0] == 0) { /* Possible little-endian */
128 if (rec1->struct_head.id[1] >= 20)
129 /* This is a big-endian file. Typically, version 2. */
130 do_swap = little_endian();
132 do_swap = big_endian();
133 } else if ((rec1->struct_head.id[1] == 0)) { /* Possible big-endian */
134 if (rec1->struct_head.id[0] <= 7)
135 /* This is a little-endian file. Version 1. */
136 do_swap = big_endian();
139 printf("DO SWAP = %d\n", do_swap);
144 short NSIG_I2 (twob x)
145 {/* x is already a pointer. */
147 memmove(&s, x, sizeof(twob));
148 if (do_swap) swap_2_bytes(&s);
152 int NSIG_I4 (fourb x)
153 { /* x is already a pointer. */
155 memmove(&i, x, sizeof(fourb));
156 if (do_swap) swap_4_bytes(&i);
161 void nsig_free_ray(NSIG_Ray *r)
163 if (r == NULL) return;
168 void nsig_free_sweep(NSIG_Sweep **s)
171 if (s == NULL) return;
172 for (itype=0; itype<s[0]->nparams; itype++) {
173 if (s[itype] == NULL) continue;
175 if (s[itype]->idh.data_type == NSIG_DTB_EXH)
176 free(s[itype]->ray[i]);
178 for (i=0; i<NSIG_I2(s[itype]->idh.num_rays_act); i++)
179 nsig_free_ray(s[itype]->ray[i]);
187 static int ipos = 0; /* Current position in the data buffer. */
188 static NSIG_Data_record data;
190 int nsig_read_chunk(FILE *fp, char *chunk)
194 int nwords, end_nwords;
196 /* The information saved from call to call is 'data' and represents
197 * the data that has been buffered. When new data is needed
198 * it will be read from the file. 'ipos' is global and initially
199 * set in 'nsig_read_sweep'. Assumptions: chunk is big enough.
201 * Return number of bytes in 'chunk' -- this variable is 'i'.
209 while(the_code != 1) {
210 if (feof(fp)) return -1;
211 if (ipos == sizeof(data)) { /* the_code is in the next chunk */
213 printf("Exceeded block size looking for the_code. Get it from next buffer.\n");
215 n = nsig_read_record(fp, (char *)data);
216 if (n <= 0) return n; /* Problem. */
219 printf("Read %d bytes.\n", n);
220 printf("Resetting ipos\n");
222 ipos = sizeof(NSIG_Raw_prod_bhdr);
226 printf("ipos = %d -- ", ipos);
228 the_code = NSIG_I2(&data[ipos]);
230 printf("the_code = %d (%d) -- ", the_code, (unsigned short)the_code);
232 ipos += sizeof(twob);
234 if (the_code < 0) { /* THIS IS DATA */
235 nwords = the_code & 0x7fff;
237 printf("#data words (2-bytes) is %d\n", nwords);
239 if (ipos + sizeof(twob)*nwords > sizeof(data)) {
241 printf("Exceeded block size... transferring and reading new (i=%d).\n", i);
243 /* Need another phyical block. */
245 /* But, first transfer the remainder to the output chunk. */
246 /* And, transfer begining of next block */
247 /* Transfer end of current buffer. */
248 end_nwords = (NSIG_BLOCK - ipos)/sizeof(twob);
249 memmove(&chunk[i], &data[ipos], sizeof(twob)*end_nwords);
250 i += end_nwords * sizeof(twob);
252 n = nsig_read_record(fp, (char *)data);
253 if (n <= 0) return n; /* Problem. */
255 nwords -= end_nwords;
256 ipos = sizeof(NSIG_Raw_prod_bhdr);
258 /* Transfer beginning of new buffer */
259 if (i+nwords * sizeof(twob) > NSIG_BLOCK) return -1;
260 memmove(&chunk[i], &data[ipos], sizeof(twob) * nwords);
261 i += nwords * sizeof(twob);
262 ipos += nwords * sizeof(twob);
265 printf("Words to transfer (at end of block) is %d\n", end_nwords);
266 printf("Transfer %d words from beginning of next buffer.\n", nwords);
267 printf("ipos in new buffer is %d\n", ipos);
269 } else { /* Normal situation. Position to end of data.
270 * But, first transfer it to the chunk.
272 if (i+nwords * sizeof(twob) > NSIG_BLOCK) return -1;
273 memmove(&chunk[i], &data[ipos], sizeof(twob) * nwords);
274 i += nwords * sizeof(twob);
275 ipos += sizeof(twob) * nwords;
278 } else if (the_code == 1) { /* END OF THE RAY. */
280 printf("------------------------------> Reached end of ray.\n");
282 break; /* or continue; */
284 } else if (the_code == 0) { /* UNKNOWN */
286 } else { /* NUMBER OF ZERO's */
288 printf("#000000000000 to skip is %d (i=%d)\n", the_code, i);
290 if (i+the_code * sizeof(twob) > NSIG_BLOCK) return -1;
291 memset(&chunk[i], 0, the_code*sizeof(twob));
292 i += the_code * sizeof(twob);
295 if (ipos >= sizeof(data)) {
297 printf("Exceeded block size ... ipos = %d\n", ipos);
298 printf("This should be right at the end of the block.\n");
300 n = nsig_read_record(fp, (char *)data);
301 if (n <= 0) return n; /* Problem. */
302 ipos = sizeof(NSIG_Raw_prod_bhdr);
309 NSIG_Ext_header_ver0 *nsig_read_ext_header_ver0(FILE *fp)
311 NSIG_Data_record chunk;
313 NSIG_Ext_header_ver0 *xh0;
315 n = nsig_read_chunk(fp, (char *)chunk);
316 if (n <= 0) return xh0;
318 printf("Ver0 x-header. %d bytes found.\n", n);
320 n -= sizeof(NSIG_Ray_header);
321 xh0 = (NSIG_Ext_header_ver0 *)calloc(1, n);
322 if (xh0) memmove(xh0, &chunk[sizeof(NSIG_Ray_header)], n);
327 NSIG_Ext_header_ver1 *nsig_read_ext_header_ver1(FILE *fp)
329 NSIG_Data_record chunk;
331 NSIG_Ext_header_ver1 *xh1;
333 n = nsig_read_chunk(fp, (char *)chunk);
334 if (n <= 0) return xh1;
336 printf("Ver1 x-header. %d bytes found.\n", n);
338 n -= sizeof(NSIG_Ray_header);
339 xh1 = (NSIG_Ext_header_ver1 *)calloc(1, n);
340 if (xh1) memmove(xh1, &chunk[sizeof(NSIG_Ray_header)], n);
344 NSIG_Ray *nsig_read_ray(FILE *fp)
347 NSIG_Ray_header rayh;
348 static NSIG_Data_record chunk;
351 n = nsig_read_chunk(fp, (char *)chunk);
352 /* Size of chunk is n */
354 if (n == 0) return NULL; /* Silent error. */
357 fprintf(stderr, "nsig_read_ray: chunk return code = %d.\n", n);
361 if (n > NSIG_BLOCK) { /* Whoa! */
362 fprintf(stderr, "nsig_read_ray: chunk bigger than buffer. n = %d,\
363 maximum block size allowed is %d\n", n, NSIG_BLOCK);
367 ray = (NSIG_Ray *) calloc(1, sizeof(NSIG_Ray));
368 memcpy(&ray->h, chunk, sizeof(NSIG_Ray_header));
369 n -= sizeof(NSIG_Ray_header);
371 printf("nsig_read_ray: allocating %d bytes for range\n", n);
373 memcpy(&rayh, chunk, sizeof(NSIG_Ray_header));
374 nbins = NSIG_I2(rayh.num_bins);
375 if (nbins <= 0) return NULL;
377 printf(" rayh.num_bins = %d (nbins %d, n %d)\n", NSIG_I2(rayh.num_bins), nbins, n);
379 ray->range = (unsigned char *)calloc(n, sizeof(unsigned char));
380 /* Changed calloc nbins to calloc n for 2-byte data.
381 ray->range = (unsigned char *)calloc(nbins, sizeof(unsigned char));
382 memmove(ray->range, &chunk[sizeof(NSIG_Ray_header)], nbins);
383 Can remove this commented-out code once we know changes work.
385 memmove(ray->range, &chunk[sizeof(NSIG_Ray_header)], n);
391 NSIG_Sweep **nsig_read_sweep(FILE *fp, NSIG_Product_file *prod_file)
395 static NSIG_Ingest_data_header **idh = NULL;
396 static NSIG_Raw_prod_bhdr *bhdr = NULL;
398 int data_mask, iray, nrays[12], max_rays;
405 NSIG_Ext_header_ver0 *exh0;
406 NSIG_Ext_header_ver1 *exh1;
409 * The organization of a RAW PRODUCT FILE: (page III-38)
411 * Record #1 { <Product_hdr> 0, 0, 0... }
412 * Record #2 { <Ingest_Summary> 0, 0, 0... }
413 * Record #3 { <Raw_Prod_BHdr> <ingest_data_header> Data...} \
414 * Record #4 { <Raw_Prod_BHdr> Data...} \
417 * Record #N { <Raw_Prod_BHdr> Data 0...} / #1
418 * Record #N+1 { <Raw_Prod_BHdr> <ingest_data_header> Data...} \
419 * Record #N+2 { <Raw_Prod_BHdr> Data...} \
422 * Record #M { <Raw_Prod_BHdr> Data 0...} / #2
424 * What about the order of info in 'Data'?
425 * Data, when it begins a sweep:
426 * a. Raw Product Bhdr
427 * b. Ingest data header for param 1
430 * Ingest data header for param n+1
434 * Ray header and Ray data are encoded with the compression algorithm.
435 * If Ray data spans more than one physical NSIG BLOCK (6144 bytes),
436 * then the 'Data' consists of:
437 * a. Raw Product Bhdr
441 * It is just missing all the Ingest data header fields.
446 /* Determine if we need to byte-swap values. */
447 (void)nsig_endianess(&prod_file->rec1);
449 /* Setup the array of ingest data headers [0..nparams-1] */
451 memmove(&masks[0], prod_file->rec2.task_config.dsp_info.data_mask_cur.mask_word_0,
453 memmove(&masks[1], &prod_file->rec2.task_config.dsp_info.data_mask_cur.mask_word_1,
456 for (j=0; j < 5; j++) {
457 data_mask = masks[j];
459 nparams += (data_mask >> i) & 0x1;
462 memmove(&data_mask, prod_file->rec2.task_config.dsp_info.data_mask, sizeof(fourb));
463 for (nparams=i=0; i<32; i++)
464 nparams += (data_mask >> i) & 0x1;
466 printf("data_mask %x\n", data_mask);
469 /* Number of sweeps */
472 nsweeps = NSIG_I2(prod_file->rec2.task_config.scan_info.num_swp);
473 printf("nsig2.c:::nparams = %d, nsweeps = %d\n", nparams, nsweeps);
480 idh = (NSIG_Ingest_data_header **)calloc(nparams, sizeof(NSIG_Ingest_data_header *));
482 for (i=0; i<nparams; i++) {
483 idh[i] = (NSIG_Ingest_data_header *)&data[sizeof(NSIG_Raw_prod_bhdr) + i*sizeof(NSIG_Ingest_data_header)];
485 bhdr = (NSIG_Raw_prod_bhdr *)prod_file->data;
488 xh_size = NSIG_I2(prod_file->rec2.ingest_head.size_ext_ray_headers);
491 rh_size = NSIG_I2(prod_file->rec2.ingest_head.size_ray_headers);
492 printf("Extended header is %d bytes long.\n", xh_size);
493 printf(" Ray header is %d bytes long.\n", rh_size);
498 max_rays = NSIG_I2(prod_file->rec2.ingest_head.num_rays);
500 /* Ingest initial block for the sweep. All remaining I/O will
501 * be performed in the de-compression loop.
503 if (feof(fp)) return NULL;
504 n = nsig_read_record(fp, (char *)data);
505 if (n <= 0) return NULL;
507 printf("Read %d bytes for data.\n", n);
511 /* This is a NEW sweep. */
515 isweep = NSIG_I2(idh[0]->sweep_num);
516 printf("Number of rays in sweep %d is %d\n", isweep, max_rays);
519 /* Allocate memory for sweep. */
520 s = (NSIG_Sweep **) calloc (nparams, sizeof(NSIG_Sweep*));
522 /* Now pointers to all possible rays. */
523 for (i=0; i<nparams; i++) {
524 /* Load sweep headers */
525 s[i] = (NSIG_Sweep *) calloc (nparams, sizeof(NSIG_Sweep));
526 s[i]->nparams = nparams;
527 memmove(&s[i]->bhdr, &bhdr, sizeof(NSIG_Raw_prod_bhdr));
528 memmove(&s[i]->idh, idh[i], sizeof(NSIG_Ingest_data_header));
529 s[i]->ray = (NSIG_Ray **) calloc (max_rays, sizeof(NSIG_Ray *));
532 /* Process this sweep. Keep track of the end of the ray. */
533 ipos = sizeof(NSIG_Raw_prod_bhdr); /* Position in the 'data' array */
536 for (i=0; i<nparams; i++) {
537 idtype[i] = NSIG_I2(idh[i]->data_type);
538 nrays[i] = (int)NSIG_I2(idh[i]->num_rays_act);
539 if (nrays[i] > max_rays) max_rays = nrays[i];
541 printf("New ray: parameter %d has idtype=%d\n", i, idtype[i]);
542 printf("Number of expected rays in sweep %d is %d\n", isweep, (int)NSIG_I2(idh[i]->num_rays_exp));
543 printf("Number of actual rays in sweep %d is %d\n", isweep, (int)NSIG_I2(idh[i]->num_rays_act));
548 ipos += nparams * sizeof(NSIG_Ingest_data_header);
550 /* ipos = sizeof(NSIG_Raw_prod_bhdr) + nparams*sizeof(NSIG_Ingest_data_header); */
551 /* 'iray' is the true ray index into 's', whereas, 'nsig_iray' is what
552 * the NSIG file says it is. I'll trust 'iray'
554 * I have a cursor into the 'data' buffer representing my current
555 * position for processing rays. This cursor will dictate if I read
556 * a new NSIG block. The cursor is call 'ipos'. It is initialized
557 * each time a new ray is encountered.
561 { int ioff, nsig_iray;
562 /* Check that all idh pointers 'id' is Ingest data header. */
563 ioff = NSIG_I2(bhdr->ray_loc);
564 nsig_iray = NSIG_I2(bhdr->ray_num);
565 printf("Offset to begining of ray %d is %d, iray=%d\n", nsig_iray, ioff,iray);
569 /* DECODE THE DATA HERE */
573 * Compression Code Meanings
575 * MSB LOW-bits Meaning
579 * 0 3-32767 3 to 32767 zeros skipped.
581 * 1 1-32767 1 to 32767 data words follow.
586 printf("---------------------- New Ray <%d> --------------------\n", iray);
588 if (feof(fp)) { /* Premature eof */
589 return NULL; /* This will have to do. */
591 /* For all parameters present. */
593 for (i=0; i<nparams; i++) {
595 /* Keep track of the cursor within the buffer, and, possibly
596 * read another buffer when a ray is split across two NSIG blocks
599 if (idtype[i] != 0) { /* Not an extended header. */
600 nsig_ray = nsig_read_ray(fp);
602 } else { /* Check extended header version. */
604 exh0 = nsig_read_ext_header_ver0(fp);
606 nsig_ray = (NSIG_Ray *)calloc(1, sizeof(NSIG_Ray));
607 nsig_ray->range = (unsigned char *)exh0;
610 exh1 = nsig_read_ext_header_ver1(fp);
612 nsig_ray = (NSIG_Ray *)calloc(1, sizeof(NSIG_Ray));
613 nsig_ray->range = (unsigned char *)exh1;
617 if (nsig_ray) is_new_ray = 1;
618 if (iray > nrays[i]) break;
619 s[i]->ray[iray] = nsig_ray;
622 if (is_new_ray) iray++;
624 } while (iray < max_rays);
626 printf("iray = %d\n", iray);
632 /**************************************************
633 * Convert 2 byte binary angle to floating point *
634 **************************************************/
635 float nsig_from_bang(bang in)
637 float result,maxval = 65536.0;
638 unsigned short bi_ang;
640 memmove(&bi_ang, in, sizeof(bang));
641 if (do_swap) swap_2_bytes(&bi_ang);
642 result = ((float)(bi_ang)/maxval) * 360.0;
647 /*************************************************
648 * convert 4 byte binary angle to floating point *
649 *************************************************/
650 float nsig_from_fourb_ang(fourb ang)
652 double result,maxval;
655 maxval = 4294967296.0;
657 memmove(&bi_ang, ang, sizeof(fourb));
658 if (do_swap) swap_4_bytes(&bi_ang);
660 result = ((double)(bi_ang)/maxval) * 360.0;
662 return ((float)result);