Pileus Git - ~andy/rsl/blob - toga_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 Darwin_Toga file and fill a Radar structure with the data.
  25  *  Darwin Type 1 data contains 4 fields: corrected & uncorrected refl,
  26  *  velocity, and spectrum width.
  27  *  Darwin type 19 data contains 1 field: uncorrected reflectivity.
  28  *  Ignore all other Darwin data types.
  29  *  Handle only PPI scans; ignore RHI scans.
  30  *
  31  * Object code from this file must be linked with the following
  32  * libraries:
  33  *    -ltg -lradar
  34  *
  35  *   "libtg.a" handles the reading and decompression of Darwin TOGA
  36  * data records from the input toga file.
  37  *
  38  *
  39  *   Kolander
  40  * Applied Research Corp.
  41  * NASA/GSFC
  42  *
  43  */
  44
  45 #include <unistd.h>
  46 #include <string.h>
  47 #include "toga.h"
  48 #include "rsl.h"
  49
  50 #define MAX_RAYS   512
  51 #define MAX_SWEEPS  20
  52 #define MISSING_VAL 0
  53
  54 extern int radar_verbose_flag;
  55 int tg_open(char *infile, tg_file_str *tg_file);
  56 int tg_read_ray(tg_file_str *tg_file);
  57
  58 /* Toga radar routines */
  59 void fill_ray(Ray *, tg_file_str *, int);
  60 void fill_ray_header(Ray *, tg_file_str *, int, int);
  61 void fill_sweep_header(Radar *, tg_map_head_str *, int, int);
  62 void fill_volume_header(Radar *, tg_map_head_str *);
  63 void fill_radar_header(Radar *, tg_map_head_str *);
  64
  65
  66
  67
  68 /********************************************************************/
  69 void fill_ray(Ray *ray, tg_file_str *tg_file, int datatype)
  70    /* transfer ray data values from the toga ray structure
  71           to the radar ray  */
  72    {
  73    int j;
  74
  75    if (ray == NULL) return;
  76    for (j=0; j<tg_file->ray.num_bins[datatype]; j++)
  77           {
  78           /* check to see if this is valid data */
  79           if (tg_file->ray.data[datatype][j] == TG_NO_DATA)
  80                  ray->range[j] = ray->h.invf(BADVAL);  /* bad data, store flag */
  81           else   /* good data, store into radar ray */
  82                  ray->range[j] = ray->h.invf(tg_file->ray.data[datatype][j]);
  83           }  /* end for */
  84    }
  85
  86
  87
  88 void fill_ray_header(Ray *ray, tg_file_str *tg_file, int elev_num,
  89                                          int datatype)
  90    /* Unresolved: Many of the following values */
  91    {
  92    if (ray == NULL) return;
  93    ray->h.month = tg_file->ray_head.mon;
  94    ray->h.day = tg_file->ray_head.day;
  95    ray->h.year = tg_file->ray_head.year;
  96    ray->h.hour = tg_file->ray_head.hour;
  97    ray->h.minute = tg_file->ray_head.min;
  98    ray->h.sec = tg_file->ray_head.hunsec/100.0;
  99    ray->h.unam_rng = MISSING_VAL;           /* ?? */
 100    ray->h.azimuth = tg_file->ray.azm;
 101    ray->h.ray_num = tg_file->ray_num;
 102    ray->h.elev = tg_file->ray.elev;
 103    ray->h.elev_num = elev_num;
 104    ray->h.range_bin1 = (int)(tg_file->ray.start_km[datatype]*1000.0);
 105    ray->h.gate_size = (int)(tg_file->ray.interval_km[datatype]*1000.0);
 106
 107    ray->h.vel_res = MISSING_VAL;           /* ?? */
 108    ray->h.sweep_rate = MISSING_VAL;        /* ?? */
 109    ray->h.prf = tg_file->map_head.prf;
 110    ray->h.azim_rate = MISSING_VAL;         /* ?? */
 111    ray->h.fix_angle = MISSING_VAL;         /* ?? */
 112    ray->h.pulse_count = MISSING_VAL;       /* ?? */
 113    ray->h.pulse_width = tg_file->map_head.pulsewd/100.0;
 114    ray->h.beam_width = 1.65;
 115    ray->h.frequency = 5625;
 116    ray->h.wavelength = (tg_file->map_head.wavelen/100.0)/100.0; /* m */
 117    ray->h.nyq_vel = (ray->h.wavelength * ray->h.prf)/4.0;  /* m/s */
 118    ray->h.nbins = tg_file->ray.num_bins[datatype];
 119
 120    if (datatype == TG_DM_IND)
 121           {
 122           ray->h.f = DZ_F;
 123           ray->h.invf = DZ_INVF;
 124           }
 125    else if (datatype == TG_DZ_IND)
 126           {
 127           ray->h.f = DZ_F;
 128           ray->h.invf = DZ_INVF;
 129           }
 130    else if (datatype == TG_VR_IND)
 131           {
 132           ray->h.f = VR_F;
 133           ray->h.invf = VR_INVF;
 134           }
 135    else if (datatype == TG_SW_IND)
 136           {
 137           ray->h.f = SW_F;
 138           ray->h.invf = SW_INVF;
 139           }
 140    }
 141
 142
 143
 144 void fill_sweep_header(Radar *radar, tg_map_head_str *map_head,
 145                                            int sweep_num, int nrays)
 146    /* Arrive here _after_ sweep ray data has been filled in. */
 147    {
 148    if (radar_verbose_flag)
 149       fprintf(stderr,"sweep_num:%02d  num_rays:%d\n",sweep_num, nrays);
 150
 151    if (map_head->data_set == 1)
 152           {
 153           radar->v[DZ_INDEX]->sweep[sweep_num]->h.f    = DZ_F;
 154           radar->v[DZ_INDEX]->sweep[sweep_num]->h.invf = DZ_INVF;
 155           radar->v[ZT_INDEX]->sweep[sweep_num]->h.f    = DZ_F;
 156           radar->v[ZT_INDEX]->sweep[sweep_num]->h.invf = DZ_INVF;
 157           radar->v[VR_INDEX]->sweep[sweep_num]->h.f    = VR_F;
 158           radar->v[VR_INDEX]->sweep[sweep_num]->h.invf = VR_INVF;
 159           radar->v[SW_INDEX]->sweep[sweep_num]->h.f    = SW_F;
 160           radar->v[SW_INDEX]->sweep[sweep_num]->h.invf = SW_INVF;
 161
 162           radar->v[DZ_INDEX]->sweep[sweep_num]->h.sweep_num = sweep_num;
 163           radar->v[ZT_INDEX]->sweep[sweep_num]->h.sweep_num = sweep_num;
 164           radar->v[VR_INDEX]->sweep[sweep_num]->h.sweep_num = sweep_num;
 165           radar->v[SW_INDEX]->sweep[sweep_num]->h.sweep_num = sweep_num;
 166
 167           radar->v[DZ_INDEX]->sweep[sweep_num]->h.elev =
 168                                       map_head->angfix[sweep_num]/10.0;
 169           radar->v[ZT_INDEX]->sweep[sweep_num]->h.elev =
 170                                       map_head->angfix[sweep_num]/10.0;
 171           radar->v[VR_INDEX]->sweep[sweep_num]->h.elev =
 172                                       map_head->angfix[sweep_num]/10.0;
 173           radar->v[SW_INDEX]->sweep[sweep_num]->h.elev =
 174                                       map_head->angfix[sweep_num]/10.0;
 175
 176           radar->v[DZ_INDEX]->sweep[sweep_num]->h.beam_width = 1.65;
 177           radar->v[ZT_INDEX]->sweep[sweep_num]->h.beam_width = 1.65;
 178           radar->v[VR_INDEX]->sweep[sweep_num]->h.beam_width = 1.65;
 179           radar->v[SW_INDEX]->sweep[sweep_num]->h.beam_width = 1.65;
 180
 181           radar->v[DZ_INDEX]->sweep[sweep_num]->h.horz_half_bw = 0.85;
 182           radar->v[ZT_INDEX]->sweep[sweep_num]->h.horz_half_bw = 0.85;
 183           radar->v[VR_INDEX]->sweep[sweep_num]->h.horz_half_bw = 0.85;
 184           radar->v[SW_INDEX]->sweep[sweep_num]->h.horz_half_bw = 0.85;
 185
 186           radar->v[DZ_INDEX]->sweep[sweep_num]->h.vert_half_bw = 0.85;
 187           radar->v[ZT_INDEX]->sweep[sweep_num]->h.vert_half_bw = 0.85;
 188           radar->v[VR_INDEX]->sweep[sweep_num]->h.vert_half_bw = 0.85;
 189           radar->v[SW_INDEX]->sweep[sweep_num]->h.vert_half_bw = 0.85;
 190
 191           radar->v[DZ_INDEX]->sweep[sweep_num]->h.nrays = nrays;
 192           radar->v[ZT_INDEX]->sweep[sweep_num]->h.nrays = nrays;
 193           radar->v[VR_INDEX]->sweep[sweep_num]->h.nrays = nrays;
 194           radar->v[SW_INDEX]->sweep[sweep_num]->h.nrays = nrays;
 195           }
 196    else if (map_head->data_set == 19)
 197           {
 198           radar->v[ZT_INDEX]->sweep[sweep_num]->h.sweep_num = sweep_num;
 199           radar->v[ZT_INDEX]->sweep[sweep_num]->h.elev =
 200                                       map_head->angfix[sweep_num]/10.0;
 201           radar->v[ZT_INDEX]->sweep[sweep_num]->h.beam_width = 1.65;
 202           radar->v[ZT_INDEX]->sweep[sweep_num]->h.horz_half_bw = 0.825;
 203           radar->v[ZT_INDEX]->sweep[sweep_num]->h.vert_half_bw = 0.825;
 204           radar->v[ZT_INDEX]->sweep[sweep_num]->h.nrays = nrays;
 205           radar->v[ZT_INDEX]->sweep[sweep_num]->h.f    = DZ_F;
 206           radar->v[ZT_INDEX]->sweep[sweep_num]->h.invf = DZ_INVF;
 207           }
 208    else  /* Another data type? */
 209           {
 210           }
 211    }
 212
 213
 214
 215 void fill_volume_header(Radar *radar, tg_map_head_str *map_head)
 216    {
 217    if (map_head->data_set == 1)
 218           {
 219           radar->v[DZ_INDEX] = RSL_new_volume(MAX_SWEEPS);
 220           radar->v[ZT_INDEX] = RSL_new_volume(MAX_SWEEPS);
 221           radar->v[VR_INDEX] = RSL_new_volume(MAX_SWEEPS);
 222           radar->v[SW_INDEX] = RSL_new_volume(MAX_SWEEPS);
 223           radar->v[DZ_INDEX]->h.type_str = strdup("Reflectivity");
 224           radar->v[ZT_INDEX]->h.type_str = strdup("Reflectivity");
 225           radar->v[VR_INDEX]->h.type_str = strdup("Velocity");
 226           radar->v[SW_INDEX]->h.type_str = strdup("Spectrum width");
 227           radar->v[DZ_INDEX]->h.nsweeps = map_head->numfix_ang;
 228           radar->v[ZT_INDEX]->h.nsweeps = map_head->numfix_ang;
 229           radar->v[VR_INDEX]->h.nsweeps = map_head->numfix_ang;
 230           radar->v[SW_INDEX]->h.nsweeps = map_head->numfix_ang;
 231           radar->v[DZ_INDEX]->h.f    = DZ_F;
 232           radar->v[DZ_INDEX]->h.invf = DZ_INVF;
 233           radar->v[ZT_INDEX]->h.f    = DZ_F;
 234           radar->v[ZT_INDEX]->h.invf = DZ_INVF;
 235           radar->v[VR_INDEX]->h.f    = VR_F;
 236           radar->v[VR_INDEX]->h.invf = VR_INVF;
 237           radar->v[SW_INDEX]->h.f    = SW_F;
 238           radar->v[SW_INDEX]->h.invf = SW_INVF;
 239           }
 240    else if (map_head->data_set == 19)
 241           {
 242           radar->v[ZT_INDEX] = RSL_new_volume(MAX_SWEEPS);
 243           radar->v[ZT_INDEX]->h.type_str = strdup("Reflectivity");
 244           radar->v[ZT_INDEX]->h.nsweeps = map_head->numfix_ang;
 245           radar->v[ZT_INDEX]->h.f    = DZ_F;
 246           radar->v[ZT_INDEX]->h.invf = DZ_INVF;
 247           }
 248    else  /* Another data type? */
 249           {
 250           }
 251    }
 252
 253
 254 void fill_radar_header(Radar *radar, tg_map_head_str *map_head)
 255    {
 256    radar->h.month = map_head->scan_mon;
 257    radar->h.day = map_head->scan_day;
 258    radar->h.year = map_head->scan_year;
 259    radar->h.hour = map_head->scan_hour;
 260    radar->h.minute = map_head->scan_min;
 261    radar->h.sec = map_head->scan_sec/100.0;
 262    strcpy(radar->h.radar_type, "toga");
 263    radar->h.number = 0;
 264    strcpy(radar->h.name, "");
 265    strcpy(radar->h.radar_name, "TOGA");
 266    strcpy(radar->h.city, "Darwin");
 267    strcpy(radar->h.state, "");
 268    /* The following lat, lon coordinates obtained from Dave Wolff */
 269    radar->h.latd = -12;
 270    radar->h.latm = -27;
 271    radar->h.lats = -26;
 272    radar->h.lond = 130;
 273    radar->h.lonm = 55;
 274    radar->h.lons = 31;
 275    radar->h.height = 37;
 276    radar->h.spulse = 500;  /* ns */
 277    radar->h.lpulse = 2000; /* ns */
 278    }
 279
 280
 281
 282
 283 Radar *RSL_toga_to_radar(char *infile)
 284    /* Ingest a Darwin_Toga file and fill a Radar structure with the data.
 285           Darwin Type 1 data contains 4 fields: corrected & uncorrected refl,
 286           velocity, and spectrum width.
 287           Darwin type 19 data contains 1 field: uncorrected reflectivity.
 288           Ignore all other Darwin data types.
 289           Handle only PPI scans; ignore RHI scans. */
 290    {
 291    int m, swp_num, num_bins;
 292    Radar *radar;
 293    Ray *new_ray;
 294    tg_file_str tg_file;
 295    extern int *rsl_qsweep; /* See RSL_read_these_sweeps in volume.c */
 296    extern int rsl_qsweep_max;
 297
 298    /* open the toga data file and read toga file header into the
 299           tg_file map_head_structure  */
 300    if (tg_open(infile,&tg_file) < 0) /* Understands NULL 'infile' as stdin. */
 301           {
 302           if (radar_verbose_flag)
 303              fprintf(stderr,"Error opening/reading data file\n");
 304           close(tg_file.fd);
 305           return NULL;
 306           }
 307    if (radar_verbose_flag)
 308           {
 309           fprintf(stderr,"Input file:  %s\n",infile);
 310           fprintf(stderr,"Scan_date: %.2d/%.2d/%d\n",tg_file.map_head.scan_mon,
 311                          tg_file.map_head.scan_day,tg_file.map_head.scan_year);
 312           fprintf(stderr,"Scan_time: %.2d:%.2d\n",tg_file.map_head.scan_hour,
 313                          tg_file.map_head.scan_min);
 314           }
 315    /* If this toga file does not contain a PPI (constant elevation)
 316           scan, do not bother with it. Just quit. */
 317    if (tg_file.map_head.scanmod != 1)
 318           {
 319           if (radar_verbose_flag)
 320              fprintf(stderr,"Darwtoga file %s does not contain a PPI scan.\n",infile);
 321           close(tg_file.fd);
 322           return NULL;
 323           }
 324    /* Can handle only datatypes 1 and 19 */
 325    if ((tg_file.map_head.data_set != 1) &&
 326            (tg_file.map_head.data_set != 19))
 327           {
 328           if (radar_verbose_flag)
 329              fprintf(stderr,"File %s does not contain Doppler or refl data\n",infile);
 330           close(tg_file.fd);
 331       return NULL;
 332           }
 333    if (radar_verbose_flag)
 334           {
 335           if (tg_file.map_head.data_set == 1) fprintf(stderr,"Type 1 data\n");
 336           else if (tg_file.map_head.data_set == 19) fprintf(stderr,"Type 19 data\n");
 337           else fprintf(stderr,"Unknown data type\n");
 338           }
 339
 340    radar = RSL_new_radar(MAX_RADAR_VOLUMES);
 341    fill_radar_header(radar, &tg_file.map_head);
 342    fill_volume_header(radar, &tg_file.map_head);
 343
 344    /* initialize counters */
 345    tg_file.ray_num = -1;
 346    swp_num = -1;
 347
 348    /* Main loop to read in a toga ray from the toga
 349           data file and store into radar structure */
 350    while ((m=tg_read_ray(&tg_file)) > 0)
 351           {
 352           /* check for end_of_sweep. Fill the sweep header _after_ we've
 353                  read in the sweep */
 354           if (swp_num < tg_file.ray_head.tilt - 1)  /* new sweep? */
 355                  {
 356                  if (swp_num >= 0)
 357                         {
 358                         /* fill the sweep header for the previous sweep */
 359                         fill_sweep_header(radar, &tg_file.map_head, swp_num,
 360                                                           tg_file.ray_num+1);
 361                         }
 362                  tg_file.ray_num = -1;  /* Reset ray_num. */
 363                  swp_num += 1;  /* increment sweep count */
 364                  if (rsl_qsweep != NULL) {
 365                    if (swp_num > rsl_qsweep_max) break;
 366                    if (rsl_qsweep[swp_num] == 0) continue;
 367                  }
 368                  /* Check for too many sweeps. */
 369                  if ((tg_file.map_head.numfix_ang < swp_num + 1) ||
 370                          (MAX_SWEEPS < swp_num + 1))
 371                         {
 372                         perror("toga_to_radar: Exceeded expected no. of sweeps");
 373                         close(tg_file.fd);
 374                         RSL_free_radar(radar);
 375                         return NULL;
 376                         }
 377                  /* Create new sweep structures. */
 378                  if (tg_file.map_head.data_set == 1)
 379                         {
 380                         radar->v[DZ_INDEX]->sweep[swp_num] = RSL_new_sweep(MAX_RAYS);
 381                         radar->v[ZT_INDEX]->sweep[swp_num] = RSL_new_sweep(MAX_RAYS);
 382                         radar->v[VR_INDEX]->sweep[swp_num] = RSL_new_sweep(MAX_RAYS);
 383                         radar->v[SW_INDEX]->sweep[swp_num] = RSL_new_sweep(MAX_RAYS);
 384                         }
 385                  else if (tg_file.map_head.data_set == 19)
 386                         {
 387                         radar->v[ZT_INDEX]->sweep[swp_num] = RSL_new_sweep(MAX_RAYS);
 388                         }
 389                  }  /* end if new sweep */
 390
 391           num_bins = tg_file.ray.num_bins[TG_DM_IND] + 4; /* types 1,19 data*/
 392           tg_file.ray_num += 1;
 393           /* Check for too many rays. */
 394           if (tg_file.ray_num > MAX_RAYS)
 395                  {
 396                  perror("toga_to_radar: Exceeded maximal no. of rays");
 397                  close(tg_file.fd);
 398                  RSL_free_radar(radar);
 399                  return NULL;
 400                  }
 401           /* check for uncorrected reflectivity data */
 402           if (tg_file.ray.da_inv[TG_DM_IND] == TRUE)
 403                  {
 404                  new_ray = RSL_new_ray(num_bins);
 405                  radar->v[ZT_INDEX]->sweep[swp_num]->ray[tg_file.ray_num] = new_ray;
 406                  fill_ray_header(new_ray, &tg_file, swp_num, TG_DM_IND);
 407                  fill_ray(new_ray, &tg_file, TG_DM_IND);
 408                  }
 409           /* check for corrected reflectivity data */
 410           if (tg_file.ray.da_inv[TG_DZ_IND] == TRUE)
 411                  {
 412                  new_ray = RSL_new_ray(num_bins);
 413                  radar->v[DZ_INDEX]->sweep[swp_num]->ray[tg_file.ray_num] = new_ray;
 414                  fill_ray_header(new_ray, &tg_file, swp_num, TG_DZ_IND);
 415                  fill_ray(new_ray, &tg_file, TG_DZ_IND);
 416                  }
 417           /* check for velocity data */
 418           if (tg_file.ray.da_inv[TG_VR_IND] == TRUE)
 419                  {
 420                  new_ray = RSL_new_ray(num_bins);
 421                  radar->v[VR_INDEX]->sweep[swp_num]->ray[tg_file.ray_num] = new_ray;
 422                  fill_ray_header(new_ray, &tg_file, swp_num, TG_VR_IND);
 423                  fill_ray(new_ray, &tg_file, TG_VR_IND);
 424                  }
 425           /* check for spectrum width data */
 426           if (tg_file.ray.da_inv[TG_SW_IND] == TRUE)
 427                  {
 428                  new_ray = RSL_new_ray(num_bins);
 429                  radar->v[SW_INDEX]->sweep[swp_num]->ray[tg_file.ray_num] = new_ray;
 430                  fill_ray_header(new_ray, &tg_file, swp_num, TG_SW_IND);
 431                  fill_ray(new_ray, &tg_file, TG_SW_IND);
 432                  }
 433           }  /* end while */
 434
 435    /* At this point, we are finished reading toga data records.
 436           Fill in the last sweep header. */
 437    fill_sweep_header(radar, &tg_file.map_head, swp_num,
 438                                          tg_file.ray_num+1);
 439
 440    /* Check which flag the TOGA read_record routines returned, and
 441           print appropriate terminating message */
 442    if (radar_verbose_flag)
 443           {
 444           switch (m)
 445                  {
 446            case TG_END_DATA:
 447                  fprintf(stderr,"Reached end of data file\n");
 448                  break;
 449            case TG_RAY_READ_ERR:
 450                  fprintf(stderr,"Error reading toga ray\n");
 451                  break;
 452            case TG_RAY_NOTYPE:
 453                  fprintf(stderr,"Found unrecognized toga ray type\n");
 454                  fprintf(stderr,"ray_type: %d\n",tg_file.ray_head.type);
 455                  break;
 456            case TG_REC_NOSEQ:
 457                  fprintf(stderr,"found out-of-sequence toga record \n");
 458                  break;
 459            case -1:
 460                  fprintf(stderr,"Error reading toga ray_header \n");
 461                  break;
 462            default:
 463                  fprintf(stderr,"Error reading toga file \n");
 464                  break;
 465                  }
 466           }  /* end if (radar_verbose_flag) */
 467
 468    close(tg_file.fd);
 469    if (m == TG_END_DATA) {
 470          radar = RSL_prune_radar(radar);
 471          return radar;
 472    } else
 473           {
 474           RSL_free_radar(radar);
 475           return NULL;
 476           }
 477
 478    }