Pileus Git - ~andy/rsl/blob - mcgill.c

   1 /*
   2     NASA/TRMM, Code 910.1.
   3     This is the TRMM Office Radar Software Library.
   4     Copyright (C) 1996, 1997
   5             Mike Kolander
   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 /******************************************************************
  25  * Opens a Mcgill format radar data file, and reads and
  26  * decompresses the contents.
  27  *
  28  * All Mcgill data structures referenced herein are
  29  * defined in the file "mcgill.h".
  30  *
  31  * There exist 3 functions in this file which are
  32  * accessible to application routines:
  33  * -----------------------------------------------------
  34  *   mcgFile_t *mcgFileOpen(int *code, char *filename);
  35  * Opens and verifies the content of a Mcgill format
  36  * data file, reads in the file header record, and
  37  * initializes the mcgFile_t structure.
  38  * -------------------------------------------------------
  39  *   int mcgFileClose(mcgFile_t *file);
  40  * Closes the Mcgill data file.
  41  * -------------------------------------------------------
  42  *   int mcgRayBuild(mcgRay_t *ray, mcgFile_t *file);
  43  * Returns in the mcgRay_t structure the reflectivity values from
  44  * one ray of data from a Mcgill format radar data file.
  45  * Successive calls return data from successive rays found
  46  * in the file.
  47  *
  48  *    Kolander
  49  *
  50  *******************************************************************/
  51
  52 #include <stdio.h>
  53 #include <stdlib.h>
  54 #include <string.h>
  55 #include "mcgill.h"
  56
  57 #define TRUE 1
  58 #define FALSE 0
  59 #define EMPTY -1.0
  60
  61
  62 /* Fixed elevation angles for Mcgill sweeps. */
  63 static float elev_angle[3][24] =
  64    {
  65    /* PAFB volume scan format types 3 and 4 (12 sweeps) */
  66           {
  67           0.6, 1.0, 1.6, 2.5, 3.6, 5.2, 7.5, 10.5, 14.7, 20.2,
  68           27.2, 35.6
  69           },
  70           /* PAFB volume scan format types 1 and 2 (24 sweeps) */
  71                  {
  72                  0.6, 0.79, 1.01, 1.27, 1.57, 1.93, 2.35, 2.84,
  73                  3.41, 4.08, 4.86, 5.78, 6.85, 8.09, 9.55, 11.23,
  74                  13.18, 15.42, 17.98, 20.89, 24.16, 27.78, 31.73, 35.97
  75                  },
  76                  /* Sao Paulo radar?? (20 sweeps) */
  77                         {
  78                         1.0, 1.3, 1.6, 2.0, 2.4, 2.9, 3.6, 4.2, 5.1,
  79                         6.1, 7.3, 8.7, 10.3, 12.2, 14.4, 16.9, 19.8,
  80                         23.1, 26.7, 30.8
  81                         }
  82    };
  83
  84 /* Fixed number_of_bins for Mcgill sweeps. */
  85 static int num_bins_normal[24] =
  86    {
  87    /* PAFB volume scan formats 1 and 3 (normal mode) */
  88    240, 180, 180, 180, 180, 180, 180, 180, 180, 180,
  89    180, 180, 180, 180, 180, 180, 180, 180, 180, 180,
  90    180, 180, 180, 180
  91    };
  92 static int num_bins_compressed[24] =
  93    {
  94    /* PAFB volume scan formats 2 and 4 (compressed mode) */
  95    180, 120, 120, 120, 120, 120, 120, 120, 120, 120,
  96    120, 120, 120, 120, 120, 120, 120, 120, 120, 120,
  97    120, 120, 120, 120
  98    };
  99
 100 FILE *uncompress_pipe (FILE *fp);
 101
 102 /**********************************************************************/
 103 mcgFile_t *mcgFileOpen(int *code, char *filename)
 104 /**********************************************************************/
 105    {
 106    /* Open a Mcgill format radar data file, and read and verify
 107       the content of the first (header) record from the file.
 108       This function returns one of the following coded integer values:
 109                  MCG_OK: Normal return.
 110          MCG_OPEN_FILE_ERR: Error occurred while attempting to open
 111                                     the data file.
 112          MCG_EOF: Unexpected End_of_File occurred reading from file.
 113          MCG_READ_ERR: Error occurred while reading from file.
 114          MCG_FORMAT_ERR: Format error encountered reading from file.
 115          SYS_NO_SPACE: Error attempting to allocate memory space.
 116    */
 117    mcgFile_t *file;
 118    char buffer[MCG_RECORD];
 119    char *csp_buffer;
 120
 121    /* Allocate space for the mcgFile structure. */
 122    if ((file = (mcgFile_t *)malloc(sizeof(mcgFile_t))) == NULL)
 123           {
 124           *code = SYS_NO_SPACE;
 125       return(NULL);
 126           }
 127
 128    /* Open Mcgill data file for reading */
 129    if ((file->fp = fopen(filename, "r")) == NULL)
 130           {
 131           *code = MCG_OPEN_FILE_ERR;
 132           return(NULL);
 133           }
 134    file->fp = uncompress_pipe(file->fp); /* Transparently, use gunzip. */
 135    /* Read first (header) record from data file into buffer */
 136    if (fread(buffer, sizeof(char), MCG_RECORD, file->fp) < MCG_RECORD)
 137           {
 138           *code = MCG_FORMAT_ERR;
 139           return(NULL);
 140           }
 141
 142    /* Move header values from buffer to the mcgHeader structure.*/
 143    memcpy(&file->head, buffer, sizeof(mcgHeader_t));
 144
 145    /* Check the leading header characters for the site ID string. */
 146    if ((strncmp((char *)&file->head, "PAFB RADAR VOLUME SCAN AT", 25) == 0) ||
 147            (strncmp((char *)&file->head, "P A B  RADAR VOLUME SCAN AT", 27) == 0))
 148           {
 149           file->site = MCG_PAFB_1_2;
 150           }
 151    /* Check if file from Sao Paulo.
 152    else if (strncmp((char *)&file->head, SAOP_ID_STRING) == 0)
 153           {
 154           file->site = MCG_SAOP;
 155           }
 156    */
 157    else  /* Can't identify site/format. */
 158           {
 159           *code = MCG_FORMAT_ERR;
 160           return(NULL);
 161           }
 162
 163    /* Check validity of some file header values. If
 164           not valid, we assume the file is garbled beyond legibility. */
 165    if ((file->head.vol_scan_format < 1) || (file->head.vol_scan_format > 4)
 166            || (file->head.day < 1) || (file->head.day > 31)
 167            || (file->head.month < 1) || (file->head.month > 12))
 168           {
 169           /* Invalid/unexpected file format */
 170           *code = MCG_FORMAT_ERR;
 171           return(NULL);
 172           }
 173
 174    /* Set the proper number_of_bins, depending on scan format. */
 175    switch (file->head.vol_scan_format)
 176           {
 177         case 1:
 178         case 3:
 179           file->num_bins = num_bins_normal;
 180           break;
 181         case 2:
 182         case 4:
 183           file->num_bins = num_bins_compressed;
 184           break;
 185         default:
 186           break;
 187           }  /* end switch */
 188
 189    /* I've never seen scan_formats 3 or 4 from PAFB, but I make
 190           allowance for them here. */
 191    if (file->site == MCG_PAFB_1_2)
 192           {
 193       if ((file->head.vol_scan_format == 3) ||
 194                   (file->head.vol_scan_format == 4))
 195                  {
 196                  file->site = MCG_PAFB_3_4;
 197                  }
 198           }
 199
 200    /* If a Mcgill Csp block is contained in this file, it's next.
 201           Skip past it. It's useless for our purposes. */
 202    if (file->head.csp_rec == 1)
 203           {
 204           /* if (fseek(file->fp, MCG_CSP, SEEK_CUR) != 0)
 205                  {
 206                  *code = MCG_READ_ERR;
 207                  return(NULL);
 208                  }
 209           */
 210           /* Allocate buffer space for the csp block, read the block
 211                  from the disk file into the buffer, then free the buffer.
 212                  An ugly way to discard the block, but the fseek solution
 213                  above won't work with stdin capabilities.
 214           */
 215           if ((csp_buffer=(char *)malloc(MCG_CSP)) == NULL)
 216                  {
 217                  *code = SYS_NO_SPACE;
 218                  return(NULL);
 219                  }
 220           if (fread(csp_buffer, sizeof(char), MCG_CSP, file->fp) < MCG_CSP)
 221                  {
 222                  *code = MCG_READ_ERR;
 223                  return(NULL);
 224                  }
 225           free(csp_buffer);
 226           }
 227
 228    /* File is open and properly initialized. */
 229    *code = MCG_OK;
 230    return(file);
 231    }
 232
 233
 234 int rsl_pclose( FILE *stream);
 235 /**********************************************************************/
 236 int mcgFileClose(mcgFile_t *file)
 237 /**********************************************************************/
 238    {
 239    /* Close a Mcgill format radar data file. Returns one of the
 240           following coded integer values:
 241              MCG_OK: File closed successfully.
 242                  MCG_CLOSE_FILE_ERR: System error occurred during file close
 243                                      operation.
 244    */
 245
 246
 247    if (rsl_pclose(file->fp) == 0)
 248           {
 249           file->fp = NULL;
 250           return(MCG_OK);
 251           }
 252    else
 253           return(MCG_CLOSE_FILE_ERR);
 254    }
 255
 256
 257
 258 /************************************************************************/
 259 int mcgRecordRead(mcgRecord_t *record, mcgFile_t *file)
 260 /************************************************************************/
 261    {
 262    /* Reads a Mcgill logical record (2048 bytes) from the data file
 263           into the mcgRecord_t structure. Returns one of the following coded
 264           integer values:
 265              MCG_OK: Successfully read a record.
 266                  MCG_EOF: End_of_File condition occurred while reading a record.
 267                  MCG_READ_ERR: Read error occurred.
 268    */
 269
 270    /* Read data from file directly into the record structure. */
 271    if (fread(record, sizeof(char), MCG_RECORD, file->fp) < MCG_RECORD)
 272           {
 273           if (feof(file->fp))
 274              return(MCG_EOF);
 275           else
 276              return(MCG_READ_ERR);
 277           }
 278
 279    return(MCG_OK);
 280    }
 281
 282
 283
 284 /************************************************************************/
 285 mcgSegmentID mcgSegmentKeyIdentify(int key)
 286 /************************************************************************/
 287    {
 288    /* Decodes the first byte from a Mcgill file segment.
 289           The first byte of the segment identifies the type of the segment.
 290           Mcgill segment types are data, elevation, and end_of_data.
 291    */
 292    if ((key > 0) && (key < 16))
 293       return(MCG_DATA_SEG);
 294    else if ((key > 31) && (key < 63))
 295       return(MCG_ELEV_SEG);
 296    else if (key == 63)
 297       return(MCG_EOD_SEG);
 298    else
 299       return(MCG_UNDEFINED_SEG);
 300    }
 301
 302
 303
 304 /*************************************************************************/
 305 int mcgRayBuild(mcgRay_t *ray, mcgFile_t *file)
 306 /*************************************************************************/
 307    {
 308    /* Returns in the ray structure the reflectivity values from one
 309           ray of data from a Mcgill format radar data file. Successive calls
 310           to this routine return data from successive rays found in the file.
 311           This function returns the following coded integer values:
 312              MCG_OK : Ray structure successfully filled with data values.
 313                  MCG_EOD: Empty ray structure, No more data in file.
 314                  MCG_EOF: Premature End_of_File condition encountered.
 315                  MCG_READ_ERR: System error occurred while reading from file.
 316                  MCG_FORMAT_ERR: Format error encountered while reading from file.
 317    */
 318
 319    /* This function is typically called about 9000 times while reading a
 320           Mcgill data file. The following static variables must retain their
 321           values between successive calls. */
 322    static int seg_num=MCG_MAX_SEG_NUM;
 323    static int eod_found = FALSE;
 324    static int sweep_num=0;
 325    static float elev=0.0;
 326    static float azm;
 327    int base, j, n;
 328    mcgSegmentID seg_type;
 329    static mcgRecord_t record;
 330
 331
 332    /* If we've previously found the end_of_data (eod) data segment in
 333           the last Mcgill record, we're done. Do this for robustness, in case
 334           the calling application routine screws up. */
 335    if (eod_found)
 336           {
 337           return(MCG_EOD);
 338           }
 339
 340    /* Initialize the ray data values. */
 341    ray->elev = elev;
 342    ray->azm = EMPTY;
 343    ray->sweep_num = sweep_num;
 344    memset(ray->data, 0, 240);
 345    ray->num_bins = file->num_bins[sweep_num];
 346
 347    /* Loop to fill the ray bins with data values from consecutive
 348           Mcgill data segments. */
 349    while (1)
 350           {
 351           /* Read in from the Mcgill file a new record, if necessary. */
 352           if (seg_num == MCG_MAX_SEG_NUM)
 353                  {
 354                  if ((n=mcgRecordRead(&record, file)) < MCG_OK)
 355                     return(n);
 356                  else
 357                         {
 358                           /* record_empty = FALSE; */
 359                         seg_num = 0;
 360                         }
 361                  }
 362
 363           /* The first byte of segment is used to identify segment type. */
 364           seg_type = mcgSegmentKeyIdentify(record.segment[seg_num][0]);
 365           switch (seg_type)
 366                  {
 367            case MCG_DATA_SEG:
 368                  /* Determine the azimuth angle for this data segment */
 369                  azm = record.segment[seg_num][1] * 64.0 +
 370                        record.segment[seg_num][2] - 1.0;
 371                  /* Check if this data segment contains data belonging to
 372                     the current ray. If so, add the segment bin values to
 373                         the ray structure. */
 374                  if ((azm == ray->azm) || (ray->azm == EMPTY))
 375                         {
 376                         ray->azm = azm;
 377                         /* Compute the ray bin base_address for the 16 bin values
 378                            from this data segment. */
 379                         base = (record.segment[seg_num][0] - 1)*16;
 380                         /* Move the 16 bin values from segment to ray structure */
 381                         for (j=0; j<16; j++)
 382                            {
 383                            ray->data[base + j] = record.segment[seg_num][3+j];
 384                            }
 385                         }
 386                  else  /* This data segment belongs to the next ray. Return the
 387                                   current ray. We'll start with this segment the next
 388                                   time mcgRayBuild() is called. */
 389                         {
 390                         return(MCG_OK);
 391                         }
 392                  break;
 393
 394            case MCG_ELEV_SEG:
 395                  sweep_num = record.segment[seg_num][0] - 31 - 1;
 396                  elev = elev_angle[file->site][sweep_num];
 397                  /* If ray structure is not empty, return it, since
 398                         we've found a Mcgill elevation segment, which initiates
 399                         a new sweep. */
 400                  if (ray->azm != EMPTY)
 401                         {
 402                         seg_num += 1;
 403                         if (seg_num == MCG_MAX_SEG_NUM)
 404                           /*                       record_empty = TRUE; */
 405                         return(MCG_OK);
 406                         }
 407                  else   /* Ray structure is empty. Continue on to read the
 408                                    next Mcgill segment.*/
 409                         {
 410                         ray->elev = elev;
 411                         ray->sweep_num = sweep_num;
 412                         }
 413                  break;
 414
 415            case MCG_EOD_SEG:
 416                  eod_found = TRUE;
 417                  if (ray->azm == EMPTY)
 418                         return(MCG_EOD);
 419                  else  /* Return the ray structure we've filled. */
 420                         return(MCG_OK);
 421                  break;
 422
 423            default:
 424                  return(MCG_FORMAT_ERR);
 425                  break;
 426                  }  /* end switch */
 427
 428           seg_num += 1;
 429           }  /* end while (1) */
 430    }