Pileus Git - ~andy/rsl/blob - mcgill_to_radar.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  *  Ingest a Mcgill format file and fill a RSL Radar structure
  25  *  with the data.
  26  *
  27  * Object code from this file must be linked with the following libraries:
  28  *     librsl  libmcg
  29  *
  30  *-----------------------------------------------------------------
  31  * NOTE: Adjust Mcgill dbz bias in function RayFill() .
  32  * Presently, if (Mcgill_dbz_value != 0)
  33  *            then Mcgill_dbz_value + 16.5 is stored in RSL Ray.
  34  *            if (Mcgill dbz value == 0)
  35  *            then 0 is stored in RSL Ray.
  36  * I don't find the Mcgill documentation very clear in this matter.
  37  *-----------------------------------------------------------------
  38  * Presently sets radar_head->radar_type = DARWIN_TOGA_SIGMET
  39  * (Need a radar_type defined in rsl.h for PAFB or MCGILL)
  40  *----------------------------------------------------------------
  41  * Minor changes will be required in the associated Mcgill library to
  42  * accomodate Mcgill data from Sao Paulo. Since I don't have any files
  43  * from Sao Paulo to test, I haven't incorporated any such changes.
  44  *-------------------------------------------------------------------
  45  * Functions defined in this file:
  46  *
  47  * void RayFill(Ray *rsl_ray, mcgRay_t *mcg_ray);
  48  * void Ray_headerInit(Ray *ray, mcgHeader_t *head,
  49  *                    mcgRay_t *mcg_ray, int ray_num, int num_bins_rsl);
  50  * void Sweep_headerInit(Sweep *sweep, mcgRay_t *mcg_ray,
  51  *                       int nrays);
  52  * void Volume_headerInit(Volume *volume, short vol_scan_format);
  53  * void Radar_headerInit(Radar *radar, mcgHeader_t *mcg_head);
  54  * Radar *RSL_mcgill_to_radar(char *infile);
  55  *
  56  *   Kolander
  57  *   09 May 95
  58  *
  59  *******************************************************************/
  60
  61 #include <string.h>
  62 #include <stdlib.h>
  63 #include "mcgill.h"
  64 #include "rsl.h"
  65
  66 #define MAX_RAYS   512
  67 #define MAX_SWEEPS  32
  68 #define MISSING_VAL 0
  69 #define MCG_DBZ_BIAS 16.5
  70 #define MCG_NOISE_BIAS 0.0
  71
  72 extern int radar_verbose_flag;
  73
  74 /*********************** Function Prototypes ***************************/
  75 static void RayFill(Ray *rsl_ray, mcgRay_t *mcg_ray);
  76 static void Ray_headerInit(Ray *ray, mcgHeader_t *head,
  77                     mcgRay_t *mcg_ray, int ray_num, int num_bins_rsl);
  78 static void Sweep_headerInit(Sweep *sweep, mcgRay_t *mcg_ray, int nrays);
  79 static void Volume_headerInit(Volume *volume, short vol_scan_format);
  80 static void Radar_headerInit(Radar *radar, mcgHeader_t *mcg_head);
  81 Radar *RSL_mcgill_to_radar(char *infile);
  82
  83 static float (*f)(Range x);
  84 static Range (*invf)(float x);
  85
  86 /* Fixed number_of_RSL_bins for Mcgill sweeps. Note that
  87    the number of bins in the RSL structure differs from
  88    the number of bins in the Mcgill ray stucture, since
  89    the km/bin spacing differs. SEE MCGILL DOCUMENTATION */
  90 static int num_bins_normal[24] =
  91    {
  92    /* PAFB volume scan formats 1 and 3 (normal mode) */
  93    480, 240, 240, 240, 240, 240, 240, 240, 240, 240,
  94    240, 240, 240, 240, 240, 240, 240, 240, 240, 240,
  95    240, 240, 240, 240
  96    };
  97 static int num_bins_compressed[24] =
  98    {
  99    /* PAFB volume scan formats 2 and 4 (compressed mode) */
 100    240, 120, 120, 120, 120, 120, 120, 120, 120, 120,
 101    120, 120, 120, 120, 120, 120, 120, 120, 120, 120,
 102    120, 120, 120, 120
 103    };
 104
 105
 106 /********************************************************************/
 107 void RayFill(Ray *rsl_ray, mcgRay_t *mcg_ray)
 108 /********************************************************************/
 109    /* transfer ray data values from the mcgRay structure
 110           to the radar->v[]->sweep[]->ray[] structure  */
 111    {
 112    int j, mcg_bin;
 113    int rsl_bin=0;
 114    int index=1;    /* Used to accomodate varying Mcgill bin spacing */
 115
 116    if (rsl_ray == NULL) return;
 117
 118    /* Note that the number of bins in the Mcgill ray structure
 119           is different from the number of bins in the RSL ray structure,
 120           since the Mcgill bin spacing changes with range. RSL ray bin
 121           spacing is constant.
 122           ---------------------------------------------------------
 123       MCGILL NORMAL FORMAT:
 124                  The first 120 mcg_bins have a length of 1 km, and span the
 125                  the range [0, 120] km.
 126                  Mcg_bins 120 to 180 have a length of 2 km, and span the
 127                  range [120, 240] km.
 128                  Mcg_bins 180 to 240 have a length of 4 km, and span the
 129                  range [240, 480] km.
 130           ---------------------------------------------------------
 131           MCGILL COMPRESSED FORMAT:
 132                  The first 120 mcg_bins have a length of 2 km, and span the
 133                  range [0, 240] km.
 134                  Mcg_bins 120 to 180 have a length of 4 km, and span the
 135                  range [240, 480] km.
 136    */
 137    for (mcg_bin=0; mcg_bin<mcg_ray->num_bins; mcg_bin++)
 138           {
 139           /* I do the case dbz_value of zero separately, so produced
 140                  color images have a black background. Otherwise, if >>all<<
 141                  dbz values are biased by 16.5, we get a colored background
 142                  where there is no precipitation. Looks peculiar. Somebody
 143                  knowledgable should look into this matter. */
 144           if (mcg_ray->data[mcg_bin] == 0)
 145                  {
 146                  for (j=0; j<index; j++)
 147                         {
 148                         rsl_ray->range[rsl_bin] = invf(MCG_NOISE_BIAS);
 149                         rsl_bin += 1;
 150                         }
 151                  }
 152           else
 153                  {
 154                  for (j=0; j<index; j++)
 155                         {
 156                         rsl_ray->range[rsl_bin] = invf(mcg_ray->data[mcg_bin] +
 157                                                                                    MCG_DBZ_BIAS);
 158                         rsl_bin += 1;
 159                         }
 160                  }  /* end else */
 161           /* Mcgill bin spacing changes at bins 120 and 180 */
 162           if ((mcg_bin == 119) || (mcg_bin == 179))
 163                  {
 164                  /* In normal format, put range rings at 120km and 240km.
 165                         In compressed format, put range rings at 240km and 480km */
 166 /*               rsl_ray->range[rsl_bin-1] = RSL_INVF(60.0); */
 167                  index = index*2;  /* Set up new bin spacing index */
 168                  }
 169           }  /* end for (mcg_bin=0...*/
 170    }
 171
 172
 173 /*************************************************************************/
 174 void Ray_headerInit(Ray *ray, mcgHeader_t *head, mcgRay_t *mcg_ray,
 175                                         int ray_num, int num_bins_rsl)
 176 /*************************************************************************/
 177    {
 178    if (ray == NULL) return;
 179    ray->h.month = (int)head->month;
 180    ray->h.day = (int)head->day;
 181    ray->h.year = 1900 + (int)head->year;
 182    if (ray->h.year < 1980) ray->h.year += 100; /* Year > 2000 */
 183    ray->h.hour = (int)head->hour;
 184    ray->h.minute = (int)head->min;
 185    ray->h.sec = (float)head->sec;
 186    ray->h.unam_rng = MISSING_VAL;           /***** ?? ******/
 187    ray->h.azimuth = mcg_ray->azm;
 188    ray->h.ray_num = ray_num;
 189    ray->h.elev = mcg_ray->elev;
 190    ray->h.elev_num = mcg_ray->sweep_num;
 191    ray->h.range_bin1 = 0;
 192    switch (head->vol_scan_format)
 193           {
 194         case 1:
 195           ray->h.gate_size = 1000;  /* 1 km/bin */
 196           ray->h.fix_angle = 24;    /* 24 sweeps */
 197           break;
 198         case 2:
 199           ray->h.gate_size = 2000;  /* 2 km/bin */
 200           ray->h.fix_angle = 24;    /* 24 sweeps */
 201           break;
 202         case 3:
 203           ray->h.gate_size = 1000;  /* 1 km/bin */
 204           ray->h.fix_angle = 12;    /* 12 sweeps */
 205           break;
 206         case 4:
 207           ray->h.gate_size = 2000;  /* 2 km/bin */
 208           ray->h.fix_angle = 12;    /* 12 sweeps */
 209           break;
 210         default:
 211           break;
 212           }  /* end switch */
 213    ray->h.vel_res = MISSING_VAL;           /* ?? */
 214    ray->h.sweep_rate = MISSING_VAL;        /* ?? */
 215    ray->h.prf = MISSING_VAL;
 216    ray->h.azim_rate = MISSING_VAL;         /* ?? */
 217    ray->h.pulse_count = MISSING_VAL;       /* ?? */
 218    ray->h.pulse_width = MISSING_VAL;       /* ?? */
 219    ray->h.frequency = MISSING_VAL;         /* ?? */
 220    ray->h.wavelength = MISSING_VAL;        /* ?? */
 221    ray->h.nyq_vel = MISSING_VAL;           /* ?? */
 222    ray->h.nbins = num_bins_rsl;
 223    ray->h.f = f;
 224    ray->h.invf = invf;
 225
 226    }
 227
 228
 229
 230 /*********************************************************************/
 231 void Sweep_headerInit(Sweep *sweep, mcgRay_t *mcg_ray, int nrays)
 232 /*********************************************************************/
 233    /* Arrive here _after_ sweep ray data has been filled in. */
 234    {
 235    if (radar_verbose_flag)
 236       fprintf(stderr,"sweep_num:%02d  num_rays:%d\n",mcg_ray->sweep_num, nrays);
 237    if (sweep == NULL) return;
 238    sweep->h.sweep_num = mcg_ray->sweep_num;
 239    sweep->h.elev = mcg_ray->elev;
 240    sweep->h.beam_width = 2.0;  /* From Dennis' function mcg2uf().
 241                                                                      I don't know where he got it. */
 242    sweep->h.horz_half_bw = 1.0;
 243    sweep->h.vert_half_bw = 1.0;
 244    sweep->h.nrays = nrays;
 245    sweep->h.f = f;
 246    sweep->h.invf = invf;
 247    }
 248
 249
 250
 251 /************************************************************************/
 252 void Volume_headerInit(Volume *volume, short vol_scan_format)
 253 /************************************************************************/
 254    {
 255    if (volume == NULL) return;
 256    volume->h.type_str = strdup("Reflectivity");
 257    volume->h.f = f;
 258    volume->h.invf = invf;
 259
 260    switch (vol_scan_format)
 261           {
 262         case 1: case 2:
 263           volume->h.nsweeps = 24;
 264           break;
 265         case 3: case 4:
 266           volume->h.nsweeps = 12;
 267           break;
 268         default:
 269           break;
 270           }
 271    }
 272
 273
 274 /***********************************************************************/
 275 void Radar_headerInit(Radar *radar, mcgHeader_t *mcg_head)
 276 /***********************************************************************/
 277    {
 278    radar->h.month = (int)mcg_head->month;
 279    radar->h.day = (int)mcg_head->day;
 280    radar->h.year = 1900 + (int)mcg_head->year;
 281    radar->h.hour = (int)mcg_head->hour;
 282    radar->h.minute = (int)mcg_head->min;
 283    radar->h.sec = (float)mcg_head->sec;
 284    strcpy(radar->h.radar_type, "mcgill");
 285    radar->h.nvolumes = 1;  /* Mcgill contains only refl data */
 286    radar->h.number = MISSING_VAL;         /****************/
 287    strcpy(radar->h.name, "PAFB");
 288    strcpy(radar->h.radar_name, "MCGILL");
 289    strcpy(radar->h.city, "MELB");
 290    strcpy(radar->h.state, "FL");
 291    radar->h.latd = 28;
 292    radar->h.latm = 15;
 293    radar->h.lats = 19;
 294    radar->h.lond = -80;
 295    radar->h.lonm = -36;
 296    radar->h.lons = -21;
 297    radar->h.height = 0;
 298    radar->h.spulse = MISSING_VAL;  /* ns */    /*************/
 299    radar->h.lpulse = MISSING_VAL;  /* ns */    /*************/
 300    }
 301
 302
 303
 304 /*********************************************************************/
 305 Radar *RSL_mcgill_to_radar(char *infile)
 306 /*********************************************************************/
 307    {
 308    /* Ingest a Mcgill format radar data file and fill a Radar RSL
 309           structure with the data.
 310    */
 311    int ray_num, code;
 312    int *num_bins_rsl;
 313    mcgFile_t *file;
 314    Radar *radar;
 315    mcgRay_t *mcg_ray, *mcg_ray_last, *swap;
 316    extern int *rsl_qsweep; /* See RSL_read_these_sweeps in volume.c */
 317    extern int rsl_qsweep_max;
 318
 319    /* If no filename has been passed, there's nothing to do. */
 320    if (infile == NULL)
 321       return NULL;
 322
 323    /* Default conversion functions. */
 324    f = DZ_F;
 325    invf = DZ_INVF;
 326
 327    /* Create a structure of type Radar */
 328    radar = (Radar *)RSL_new_radar(MAX_RADAR_VOLUMES);
 329    if (radar == NULL) {
 330          perror("RSL_mcgill_to_radar: RSL_new_radar:");
 331          return NULL;
 332    }
 333
 334
 335    /* Open the Mcgill data file and read Mcgill file header into the
 336           mcgFile.head structure */
 337    file = (mcgFile_t *)mcgFileOpen(&code, infile);
 338    if (file == NULL)
 339       goto quit;
 340
 341    if (radar_verbose_flag)
 342           {
 343           fprintf(stderr,"Input file:  %s\n", infile);
 344           fprintf(stderr,"Scan_date: %d/%d/%d\n", file->head.month, file->head.day,
 345                          file->head.year);
 346           fprintf(stderr,"Scan_time: %d:%d:%d\n", file->head.hour, file->head.min,
 347                          file->head.sec);
 348           fprintf(stderr,"num_recs:%d  format:%d\n\n",
 349                          file->head.num_records, file->head.vol_scan_format);
 350           }
 351
 352    /* Initialize num_bins_rsl to point to the proper bin count array for
 353           the scan format. */
 354    if ((file->head.vol_scan_format == 1) ||
 355            (file->head.vol_scan_format == 3)) /* Normal mode */
 356           {
 357           num_bins_rsl = num_bins_normal;
 358           }
 359    else  /* Compressed mode */
 360       num_bins_rsl = num_bins_compressed;
 361
 362    /* Allocate space for two mcg_rays. mcg_ray will contain the current
 363           ray, and mcg_ray_last will contain the previous ray. The previous ray
 364           is needed only at sweep increments. */
 365    mcg_ray = (mcgRay_t *)malloc(sizeof(mcgRay_t));
 366    mcg_ray_last = (mcgRay_t *)malloc(sizeof(mcgRay_t));
 367
 368    /* Initialize the values of the Radar_header structure */
 369    Radar_headerInit(radar, &file->head);
 370    /* Create a structure of type Volume */
 371    radar->v[DZ_INDEX] = RSL_new_volume(MAX_SWEEPS);
 372    if (radar->v[DZ_INDEX] == NULL) {
 373          perror("mcgill_to_radar: RSL_new_volume");
 374          return radar;
 375    }
 376    /* Initialize the values of the Volume_header structure. */
 377    Volume_headerInit(radar->v[DZ_INDEX], file->head.vol_scan_format);
 378
 379    /* initialize counters */
 380    ray_num = -1;
 381    mcg_ray_last->sweep_num = -1;
 382
 383    /* Main loop to read in a ray from the Mcgill data file and store
 384           into RSL radar structure */
 385    while ((code = mcgRayBuild(mcg_ray, file)) == MCG_OK)
 386           {
 387           /* Discard rays with bogus azimuth values. */
 388           if (mcg_ray->azm > 360.0)
 389                  {
 390                  if (radar_verbose_flag)
 391                     fprintf(stderr,"**** Bogus azm:%.1f, discarding ray.\n", mcg_ray->azm);
 392              continue;
 393                  }
 394           /* Check for end_of_sweep. Fill the sweep header _after_ we've
 395                  read in the sweep, since we don't know the number of
 396                  rays in the sweep until we find the start of the next sweep*/
 397           if (mcg_ray_last->sweep_num < mcg_ray->sweep_num)  /* new sweep? */
 398                  {
 399                  if (mcg_ray_last->sweep_num >= 0)
 400                         {
 401                         /* fill the sweep header for the previous sweep */
 402                         Sweep_headerInit(radar->v[DZ_INDEX]->sweep[mcg_ray_last->sweep_num],
 403                                                           mcg_ray_last, ray_num+1);
 404                         }
 405                  ray_num = -1;  /* Reset ray counter at start of each sweep */
 406
 407                  /* Check for too many sweeps. If we've exceeded MAX_SWEEPS,
 408                     we're lost. */
 409                  if (MAX_SWEEPS < mcg_ray->sweep_num + 1)
 410                         {
 411                         perror("RSL_mcgill_to_radar: Exceeded expected no. of sweeps");
 412                         mcgFileClose(file);
 413                         RSL_free_radar(radar);
 414                         return NULL;
 415                         }
 416                  if (rsl_qsweep != NULL) {
 417                    if (mcg_ray->sweep_num > rsl_qsweep_max) break;
 418                    if (rsl_qsweep[mcg_ray->sweep_num] == 0) continue;
 419                  }
 420
 421                  /* Create new sweep structure. */
 422                  radar->v[DZ_INDEX]->sweep[mcg_ray->sweep_num] = RSL_new_sweep(MAX_RAYS);
 423                  if (radar->v[DZ_INDEX]->sweep[mcg_ray->sweep_num] == NULL) {
 424                    perror("mcgill_to_radar: RSL_new_sweep");
 425                    return radar;
 426                  }
 427                  }  /* end if new sweep */
 428
 429           ray_num += 1;  /* Increment ray counter */
 430           /* Check for too many rays. */
 431           if (ray_num > MAX_RAYS)
 432                  {
 433                  perror("mcgill_to_radar: Exceeded maximal no. of rays");
 434                  mcgFileClose(file);
 435                  RSL_free_radar(radar);
 436                  return NULL;
 437                  }
 438
 439           /* Create a new RSL ray containing correct # of bins, and fill it. */
 440           radar->v[DZ_INDEX]->sweep[mcg_ray->sweep_num]->ray[ray_num] =
 441                                      RSL_new_ray(num_bins_rsl[mcg_ray->sweep_num] + 8);
 442           if (radar->v[DZ_INDEX]->sweep[mcg_ray->sweep_num]->ray[ray_num] == NULL ) {
 443                 perror("mcgill_to_radar: RSL_new_ray");
 444                 return radar;
 445           }
 446
 447           Ray_headerInit(radar->v[DZ_INDEX]->sweep[mcg_ray->sweep_num]->ray[ray_num],
 448                                          &file->head, mcg_ray, ray_num, num_bins_rsl[mcg_ray->sweep_num]);
 449           RayFill(radar->v[DZ_INDEX]->sweep[mcg_ray->sweep_num]->ray[ray_num],
 450                           mcg_ray);
 451
 452           /* Swap mcg_ray pointers so that the current mcg_ray, which we've
 453                  finished loading into an RSL_ray, becomes mcg_ray_last . */
 454           swap = (mcgRay_t *)mcg_ray_last;
 455           mcg_ray_last = (mcgRay_t *)mcg_ray;
 456           mcg_ray = (mcgRay_t *)swap;
 457           }  /* end while */
 458
 459    /* At this point, we're finished reading file records.
 460           Fill in the last sweep header. */
 461    Sweep_headerInit(radar->v[DZ_INDEX]->sweep[mcg_ray_last->sweep_num],
 462                                         mcg_ray_last, ray_num+1);
 463
 464    /* Check which flag the mcgill routines returned, and
 465           print appropriate terminating message */
 466 quit: if (radar_verbose_flag)
 467           {
 468           switch (code)
 469                  {
 470            case MCG_EOD:
 471                  fprintf(stderr,"MCG_EOD: Reached end of data file\n");
 472                  break;
 473            case MCG_OPEN_FILE_ERR:
 474                  fprintf(stderr,"MCG_OPEN_FILE_ERR: Error opening Mcgill data file\n");
 475                  break;
 476            case MCG_EOF:
 477                  fprintf(stderr,"MCG_EOF: End of data file reached prematurely\n");
 478                  break;
 479            case MCG_READ_ERR:
 480                  fprintf(stderr,"MCG_READ_ERR: Error occurred while reading from data file\n");
 481                  break;
 482            case MCG_FORMAT_ERR:
 483                  fprintf(stderr,"MCG_FORMAT_ERR: Format error encountered in data file\n");
 484                  break;
 485            default:
 486                  fprintf(stderr,"Error reading data file \n");
 487                  break;
 488                  }
 489           }  /* end if (radar_verbose_flag) */
 490
 491    if (code == MCG_EOD)  /* Successfully read in Mcgill file? */
 492           {
 493           mcgFileClose(file);
 494           radar = RSL_prune_radar(radar);
 495       return(radar);
 496           }
 497    else  /* Fatal error occurred */
 498           {
 499           RSL_free_radar(radar);
 500           return(NULL);
 501           }
 502    }
 503