]> Pileus Git - ~andy/rsl/blob - nsig.c
038b7de6b16e79fb0056b61393fafd2a751c0128
[~andy/rsl] / nsig.c
1 /*
2     NASA/TRMM, Code 910.1.
3     This is the TRMM Office Radar Software Library.
4     Copyright (C) 1996, 1997
5             John H. Merritt
6             Space Applications Corporation
7             Vienna, Virginia
8
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.
13
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.
18
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.
22 */
23 /* 
24  * Read SIGMET version 1 and version 2 formatted files.
25  *
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.
29  *
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.
33  *
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
37  *                       end-of-file.
38  *----------------------------------------------------------------------
39  * 8/13/96
40  *
41  * John H. Merritt
42  * Space Applications Corp.
43  * NASA/GSFC Code 910.1 
44  *
45  * Copyright 1996
46  */
47
48 #include <stdio.h>
49 #include <string.h>
50 #include <stdlib.h>
51 #include <unistd.h>
52
53 #include "nsig.h"
54
55 FILE *uncompress_pipe(FILE *fp);
56 int big_endian(void);
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);
61
62 /*********************************************************************
63  * Open a file and possibly setup a gunzip pipe                      *
64  *********************************************************************/
65 FILE *nsig_open(char *file_name)
66 {
67   FILE *fp;
68   int save_fd;
69
70    /* Open input file */
71   if (file_name == NULL) { /* Use stdin */
72         save_fd = dup(0);
73         fp = fdopen(save_fd, "r");
74   } else if((fp = fopen(file_name,"r")) == NULL) {
75         perror(file_name);
76         return fp;
77   }
78
79   fp = uncompress_pipe(fp); /* Transparently gunzip. */
80   return fp;
81 }
82
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)
88 {
89    int n;
90    int nbytes;
91    char *buf;
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.
95     */
96    /* Return the number of bytes read. */
97    buf = (char *)nsig_rec;
98
99    if (feof(fp)) return -1;
100    nbytes = 0;
101    while((n = fread(&buf[nbytes], sizeof(char),NSIG_BLOCK-nbytes, fp)) > 0) {
102          nbytes += n;
103    }
104    return nbytes;
105 }
106
107
108 /*********************************************************
109  * Close nsig file                                       *
110  *********************************************************/
111 void nsig_close(FILE *fp)
112    {
113          rsl_pclose(fp);
114    }
115
116 static int do_swap;
117
118 int nsig_endianess(NSIG_Record1 *rec1)
119 {
120   /*
121    * If NSIG is version1 and on big-endian, then swap.
122    * If NSIG is version2 and on little-endian, then swap.
123    */
124   /*
125   printf("id = %d %d\n", (int)rec1->struct_head.id[0], (int)rec1->struct_head.id[1]);
126   */
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();
131         else
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();
137   }
138   /*
139   printf("DO SWAP = %d\n", do_swap);
140   */
141   return do_swap;
142 }
143
144 short NSIG_I2 (twob x)
145 {/* x is already a pointer. */
146   short s;
147   memmove(&s, x, sizeof(twob));
148   if (do_swap) swap_2_bytes(&s);
149   return s;
150 }
151
152 int NSIG_I4 (fourb x)
153 { /* x is already a pointer. */
154   int i;
155   memmove(&i, x, sizeof(fourb));
156   if (do_swap) swap_4_bytes(&i);
157   return i;
158 }
159
160
161 void nsig_free_ray(NSIG_Ray *r)
162 {
163   if (r == NULL) return;
164   free(r->range);
165   free(r);
166 }
167
168 void nsig_free_sweep(NSIG_Sweep **s)
169 {
170   int i=0,itype;
171   if (s == NULL) return;
172   for (itype=0; itype<s[0]->nparams; itype++) {
173         if (s[itype] == NULL) continue;
174         
175         if (s[itype]->idh.data_type == NSIG_DTB_EXH)
176           free(s[itype]->ray[i]);
177         else {
178           for (i=0; i<NSIG_I2(s[itype]->idh.num_rays_act); i++)
179                 nsig_free_ray(s[itype]->ray[i]);
180         }
181         free(s[itype]->ray);
182         free(s[itype]);
183   }
184   free(s);
185 }
186
187 static int ipos = 0;  /* Current position in the data buffer. */
188 static NSIG_Data_record data;
189
190 int nsig_read_chunk(FILE *fp, char *chunk)
191 {
192   int i, n;
193   int the_code;
194   int nwords, end_nwords;
195
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.
200    *
201    * Return number of bytes in 'chunk'  -- this variable is 'i'.
202    */
203
204   i = n = 0;
205   the_code = 0;
206
207 #define Vprint
208 #undef  Vprint
209   while(the_code != 1) {
210         if (feof(fp)) return -1;
211         if (ipos == sizeof(data)) { /* the_code is in the next chunk */
212 #ifdef Vprint
213           printf("Exceeded block size looking for the_code. Get it from next buffer.\n");
214 #endif
215           n = nsig_read_record(fp, (char *)data);
216           if (n <= 0) return n;  /* Problem. */
217
218 #ifdef Vprint
219           printf("Read %d bytes.\n", n);
220           printf("Resetting ipos\n");
221 #endif
222           ipos = sizeof(NSIG_Raw_prod_bhdr);
223         }
224           
225 #ifdef Vprint
226         printf("ipos = %d -- ", ipos);
227 #endif
228         the_code = NSIG_I2(&data[ipos]);
229 #ifdef Vprint
230         printf("the_code = %d (%d) -- ", the_code, (unsigned short)the_code);
231 #endif
232         ipos += sizeof(twob);
233         
234         if (the_code < 0) {  /* THIS IS DATA */
235           nwords = the_code & 0x7fff;
236 #ifdef Vprint
237           printf("#data words (2-bytes) is %d\n", nwords);
238 #endif
239           if (ipos + sizeof(twob)*nwords > sizeof(data)) {
240 #ifdef Vprint
241                 printf("Exceeded block size... transferring and reading new (i=%d).\n", i);
242 #endif
243                 /* Need another phyical block. */
244                         
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);
251
252                 n = nsig_read_record(fp, (char *)data);
253                 if (n <= 0) return n;  /* Problem. */
254                 /* New ipos */
255                 nwords -= end_nwords;
256                 ipos = sizeof(NSIG_Raw_prod_bhdr);
257                         
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);
263                 
264 #ifdef Vprint
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);
268 #endif
269           } else { /* Normal situation.  Position to end of data.
270                                 * But, first transfer it to the chunk.
271                           */
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;
276           }
277                   
278         } else if (the_code == 1) {  /* END OF THE RAY. */
279 #ifdef Vprint
280           printf("------------------------------> Reached end of ray.\n");
281 #endif
282           break; /* or continue; */
283
284         } else if (the_code == 0) {  /* UNKNOWN */
285           break;
286         } else {  /* NUMBER OF ZERO's */
287 #ifdef Vprint
288           printf("#000000000000 to skip is %d (i=%d)\n", the_code, i);
289 #endif
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);
293         }
294                 
295         if (ipos >= sizeof(data)) {
296 #ifdef Vprint
297           printf("Exceeded block size ... ipos = %d\n", ipos);
298           printf("This should be right at the end of the block.\n");
299 #endif
300           n = nsig_read_record(fp, (char *)data);
301           if (n <= 0) return n; /* Problem. */
302           ipos = sizeof(NSIG_Raw_prod_bhdr);
303         }
304         
305   } /* End while. */                                                     
306   return i;
307 }
308
309 NSIG_Ext_header_ver0 *nsig_read_ext_header_ver0(FILE *fp)
310 {
311   NSIG_Data_record chunk;
312   int n;
313   NSIG_Ext_header_ver0 *xh0;
314   xh0 = NULL;
315   n = nsig_read_chunk(fp, (char *)chunk);
316   if (n <= 0) return xh0;
317 #ifdef Vprint
318   printf("Ver0 x-header.  %d bytes found.\n", n);
319 #endif
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);
323   return xh0;
324 }
325
326
327 NSIG_Ext_header_ver1 *nsig_read_ext_header_ver1(FILE *fp)
328 {
329   NSIG_Data_record chunk;
330   int n;
331   NSIG_Ext_header_ver1 *xh1;
332   xh1 = NULL;
333   n = nsig_read_chunk(fp, (char *)chunk);
334   if (n <= 0) return xh1;
335 #ifdef Vprint
336   printf("Ver1 x-header.  %d bytes found.\n", n);
337 #endif
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);
341   return xh1;
342 }
343
344 NSIG_Ray *nsig_read_ray(FILE *fp)
345 {
346   int n, nbins;
347   NSIG_Ray_header rayh;
348   static NSIG_Data_record chunk;
349   NSIG_Ray *ray;
350   
351   n = nsig_read_chunk(fp, (char *)chunk);
352   /* Size of chunk is n */
353
354   if (n == 0) return NULL; /* Silent error. */
355
356   if (n < 0) {
357         fprintf(stderr, "nsig_read_ray: chunk return code = %d.\n", n);
358         return NULL;
359   }
360
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);
364         return NULL;
365   }
366
367   ray = (NSIG_Ray *) calloc(1, sizeof(NSIG_Ray));
368   memcpy(&ray->h, chunk, sizeof(NSIG_Ray_header));
369   n -= sizeof(NSIG_Ray_header);
370 #ifdef Vprint
371   printf("nsig_read_ray: allocating %d bytes for range\n", n);
372 #endif
373   memcpy(&rayh, chunk, sizeof(NSIG_Ray_header));
374   nbins = NSIG_I2(rayh.num_bins);
375   if (nbins <= 0) return NULL;
376 #ifdef Vprint
377   printf("               rayh.num_bins = %d (nbins %d, n %d)\n", NSIG_I2(rayh.num_bins), nbins, n);
378 #endif
379   ray->range = (unsigned char *)calloc(nbins, sizeof(unsigned char));
380   memmove(ray->range, &chunk[sizeof(NSIG_Ray_header)], nbins);
381   
382   return ray;
383 }
384
385
386 NSIG_Sweep **nsig_read_sweep(FILE *fp, NSIG_Product_file *prod_file)
387 {
388   NSIG_Sweep **s;
389   int i, n;
390   static NSIG_Ingest_data_header **idh = NULL;
391   static NSIG_Raw_prod_bhdr *bhdr = NULL;
392   NSIG_Ray *nsig_ray;
393   int data_mask, iray, nrays[12], max_rays;
394   int nparams;
395   int is_new_ray;
396   int idtype[12];
397   int is_new_sweep;
398   int xh_size;
399   NSIG_Ext_header_ver0 *exh0;
400   NSIG_Ext_header_ver1 *exh1;
401
402   /*
403    * The organization of a RAW PRODUCT FILE:  (page III-38)
404    *
405    *  Record #1   { <Product_hdr> 0, 0, 0... }
406    *  Record #2   { <Ingest_Summary> 0, 0, 0... }
407    *  Record #3   { <Raw_Prod_BHdr> <ingest_data_header> Data...} \
408    *  Record #4   { <Raw_Prod_BHdr> Data...}                       \
409    *      .               .           .                            | Data for
410    *      .               .           .                            /  Sweep
411    *  Record #N   { <Raw_Prod_BHdr> Data 0...}                    /     #1
412    *  Record #N+1 { <Raw_Prod_BHdr> <ingest_data_header> Data...} \
413    *  Record #N+2 { <Raw_Prod_BHdr> Data...}                       \
414    *      .               .           .                            | Data for
415    *      .               .           .                            /  Sweep
416    *  Record #M   { <Raw_Prod_BHdr> Data 0...}                    /     #2
417    *
418    *  What about the order of info in 'Data'?
419    *    Data, when it begins a sweep:
420    *         a. Raw Product Bhdr
421    *         b. Ingest data header for param 1
422    *                .
423    *                .
424    *            Ingest data header for param n+1
425    *         c. Ray header
426    *         d. Ray data
427    *
428    *    Ray header and Ray data are encoded with the compression algorithm.
429    *    If Ray data spans more than one physical NSIG BLOCK (6144 bytes),
430    *    then the 'Data' consists of:
431    *         a. Raw Product Bhdr
432    *         b. Ray header
433    *         c. Ray data
434    *
435    *    It is just missing all the Ingest data header fields.
436    */
437
438 #define Vprint
439 #undef  Vprint
440   /* Determine if we need to byte-swap values. */
441   (void)nsig_endianess(&prod_file->rec1);
442   
443   /* Setup the array of ingest data headers [0..nparams-1] */
444   memmove(&data_mask, prod_file->rec2.task_config.dsp_info.data_mask, sizeof(fourb));
445 #ifdef Vprint
446   printf("data_mask %x\n", data_mask);
447 #endif
448   for (nparams=i=0; i<32; i++)
449         nparams += (data_mask >> i) & 0x1;
450   /* Number of sweeps */
451 #ifdef Vprint
452   {int nsweeps;
453   nsweeps = NSIG_I2(prod_file->rec2.task_config.scan_info.num_swp);
454   printf("nsig2.c:::nparams = %d, nsweeps = %d\n", nparams, nsweeps);
455   }
456 #endif
457
458
459   if (idh == NULL) {
460         
461         idh = (NSIG_Ingest_data_header **)calloc(nparams, sizeof(NSIG_Ingest_data_header *));
462         ipos = 0;
463         for (i=0; i<nparams; i++) {
464           idh[i] = (NSIG_Ingest_data_header *)&data[sizeof(NSIG_Raw_prod_bhdr) + i*sizeof(NSIG_Ingest_data_header)];
465         }
466         bhdr = (NSIG_Raw_prod_bhdr *)prod_file->data;
467   }
468
469   xh_size = NSIG_I2(prod_file->rec2.ingest_head.size_ext_ray_headers);
470 #ifdef Vprint
471   {int rh_size;
472   rh_size = NSIG_I2(prod_file->rec2.ingest_head.size_ray_headers);
473   printf("Extended header is %d bytes long.\n", xh_size);
474   printf("     Ray header is %d bytes long.\n", rh_size);
475   }
476 #endif
477   is_new_ray = 1;
478   is_new_sweep = 1;
479   max_rays  = NSIG_I2(prod_file->rec2.ingest_head.num_rays);
480
481   /* Ingest initial block for the sweep.  All remaining I/O will
482    * be performed in the de-compression loop.
483    */
484   if (feof(fp)) return NULL;
485   n = nsig_read_record(fp, (char *)data);
486   if (n <= 0) return NULL;
487 #ifdef Vprint
488   printf("Read %d bytes for data.\n", n);
489 #endif
490
491
492   /* This is a NEW sweep. */
493   iray = 0;
494 #ifdef Vprint
495   {int isweep;
496   isweep = NSIG_I2(idh[0]->sweep_num);
497   printf("Number of rays in sweep %d is %d\n", isweep, max_rays);
498   }
499 #endif
500   /* Allocate memory for sweep. */
501   s = (NSIG_Sweep **) calloc (nparams, sizeof(NSIG_Sweep*));
502
503   /* Now pointers to all possible rays. */
504   for (i=0; i<nparams; i++) {
505   /* Load sweep headers */
506         s[i] = (NSIG_Sweep *) calloc (nparams, sizeof(NSIG_Sweep));
507         s[i]->nparams = nparams;
508         memmove(&s[i]->bhdr, &bhdr, sizeof(NSIG_Raw_prod_bhdr));
509         memmove(&s[i]->idh, idh[i], sizeof(NSIG_Ingest_data_header));
510         s[i]->ray = (NSIG_Ray **) calloc (max_rays, sizeof(NSIG_Ray *));
511   }     
512
513   /* Process this sweep. Keep track of the end of the ray. */
514   ipos = sizeof(NSIG_Raw_prod_bhdr); /* Position in the 'data' array */
515
516   max_rays = 0;
517   for (i=0; i<nparams; i++) {
518         idtype[i] = NSIG_I2(idh[i]->data_type);
519         nrays[i] = (int)NSIG_I2(idh[i]->num_rays_act);
520         if (nrays[i] > max_rays) max_rays = nrays[i];
521 #ifdef Vprint
522         printf("New ray: parameter %d has idtype=%d\n", i, idtype[i]);
523         printf("Number of expected rays in sweep %d is %d\n", isweep, (int)NSIG_I2(idh[i]->num_rays_exp));
524         printf("Number of actual   rays in sweep %d is %d\n", isweep,  (int)NSIG_I2(idh[i]->num_rays_act));
525 #endif
526
527   }
528   if (is_new_sweep) 
529         ipos += nparams * sizeof(NSIG_Ingest_data_header);
530   
531   /*    ipos = sizeof(NSIG_Raw_prod_bhdr) + nparams*sizeof(NSIG_Ingest_data_header); */
532   /* 'iray' is the true ray index into 's', whereas, 'nsig_iray' is what
533    * the NSIG file says it is.  I'll trust 'iray'
534    *
535    * I have a cursor into the 'data' buffer representing my current
536    * position for processing rays.  This cursor will dictate if I read
537    * a new NSIG block.  The cursor is call 'ipos'.  It is initialized
538    * each time a new ray is encountered.
539    */
540
541 #ifdef Vprint
542   { int ioff, nsig_iray;
543   /* Check that all idh pointers 'id' is Ingest data header. */
544   ioff = NSIG_I2(bhdr->ray_loc);
545   nsig_iray = NSIG_I2(bhdr->ray_num);
546   printf("Offset to begining of ray %d is %d, iray=%d\n", nsig_iray, ioff,iray);
547   }
548 #endif
549   
550   /* DECODE THE DATA HERE */
551   /* From III-39 */
552   /*
553    *       Table 3.5-5
554    *   Compression Code Meanings
555    *
556    * MSB LOW-bits  Meaning
557    *  0     0      <unused>
558    *  0     1      End of ray.
559    *  0     2      <unused>
560    *  0   3-32767  3 to 32767 zeros skipped.
561    *  1     0      <unused>
562    *  1   1-32767  1 to 32767 data words follow.
563    */
564   
565   do {
566 #ifdef Vprint
567         printf("---------------------- New Ray <%d> --------------------\n", iray);
568 #endif
569         if (feof(fp)) { /* Premature eof */
570           return NULL; /* This will have to do. */
571         }
572         /* For all parameters present. */
573         is_new_ray = 0;
574         for (i=0; i<nparams; i++) {
575
576           /* Keep track of the cursor within the buffer, and, possibly
577        * read another buffer when a ray is split across two NSIG blocks
578            */
579           nsig_ray = NULL;
580           if (idtype[i] != 0) { /* Not an extended header. */
581                 nsig_ray = nsig_read_ray(fp);
582
583           } else { /* Check extended header version. */
584                 if (xh_size <= 20) {
585                   exh0 = nsig_read_ext_header_ver0(fp);
586                   if (exh0) {
587                         nsig_ray = (NSIG_Ray *)calloc(1, sizeof(NSIG_Ray));
588                         nsig_ray->range = (unsigned char *)exh0;
589                   }
590                 } else {
591                   exh1 = nsig_read_ext_header_ver1(fp);
592                   if (exh1) {
593                         nsig_ray = (NSIG_Ray *)calloc(1, sizeof(NSIG_Ray));
594                         nsig_ray->range = (unsigned char *)exh1;
595                   }
596                 }
597           }
598           if (nsig_ray) is_new_ray = 1;
599           if (iray > nrays[i]) break;
600           s[i]->ray[iray] = nsig_ray;
601
602         } /* End for */
603         if (is_new_ray) iray++;
604         
605   } while (iray < max_rays);
606 #ifdef Vprint
607   printf("iray = %d\n", iray);
608 #endif
609   return s;
610   
611 }
612
613 /**************************************************
614  *  Convert 2 byte binary angle to floating point *
615  **************************************************/
616 float nsig_from_bang(bang in)
617    {
618    float result,maxval = 65536.0;
619    unsigned short bi_ang;
620
621    memmove(&bi_ang, in, sizeof(bang));
622    if (do_swap) swap_2_bytes(&bi_ang);
623    result = ((float)(bi_ang)/maxval) * 360.0;
624    
625    return (result);
626    }
627
628 /*************************************************
629  * convert 4 byte binary angle to floating point *
630  *************************************************/
631 float nsig_from_fourb_ang(fourb ang)
632    {
633    double result,maxval;
634    unsigned int bi_ang;
635
636    maxval = 4294967296.0;
637
638    memmove(&bi_ang, ang, sizeof(fourb));
639    if (do_swap) swap_4_bytes(&bi_ang);
640    
641    result = ((double)(bi_ang)/maxval) * 360.0;
642    
643    return ((float)result);
644    }