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

   1 /*
   2     NASA/TRMM, Code 910.1.
   3     This is the TRMM Office Radar Software Library.
   4     Copyright (C) 1996-1999
   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 #include <stdio.h>
  25 #include <stdlib.h>
  26 #include <unistd.h>
  27 #include <string.h>
  28 #define USE_RSL_VARS
  29 #include "rsl.h"
  30 #include "dorade.h"
  31
  32 extern int radar_verbose_flag;
  33
  34 /********************************************************************/
  35 /*                                                                  */
  36 /*                    find_rsl_field_index                          */
  37 /*                                                                  */
  38 /********************************************************************/
  39 int find_rsl_field_index(char *dorade_field_name)
  40 {
  41   /*
  42    * Dorade: VE, DM, SW, DBZ, ZDR, PHI, RHO, LDR,  DX,  CH,  AH,  CV,  AV
  43    *    RSL: VR, DM, SW, DZ,  ZD,  PH,  RH,  LR,  *DX, *CH, *AH, *CV, *AV.
  44    */
  45   if (strncasecmp(dorade_field_name, "ve", 2) == 0)  return VR_INDEX;
  46   if (strncasecmp(dorade_field_name, "vr", 2) == 0)  return VR_INDEX;
  47   if (strncasecmp(dorade_field_name, "dm", 2) == 0)  return DM_INDEX;
  48   if (strncasecmp(dorade_field_name, "sw", 2) == 0)  return SW_INDEX;
  49   if (strncasecmp(dorade_field_name, "dbz", 3) == 0) return DZ_INDEX;
  50   if (strncasecmp(dorade_field_name, "zdr", 3) == 0) return ZD_INDEX;
  51   if (strncasecmp(dorade_field_name, "phi", 3) == 0) return PH_INDEX;
  52   if (strncasecmp(dorade_field_name, "rho", 3) == 0) return RH_INDEX;
  53   if (strncasecmp(dorade_field_name, "ldr", 3) == 0) return LR_INDEX;
  54   if (strncasecmp(dorade_field_name, "dx", 2) == 0)  return DX_INDEX;
  55   if (strncasecmp(dorade_field_name, "ch", 2) == 0)  return CH_INDEX;
  56   if (strncasecmp(dorade_field_name, "ah", 2) == 0)  return AH_INDEX;
  57   if (strncasecmp(dorade_field_name, "cv", 2) == 0)  return CV_INDEX;
  58   if (strncasecmp(dorade_field_name, "av", 2) == 0)  return AV_INDEX;
  59   if (strncasecmp(dorade_field_name, "vs", 2) == 0)  return VS_INDEX;
  60   if (strncasecmp(dorade_field_name, "vl", 2) == 0)  return VL_INDEX;
  61   if (strncasecmp(dorade_field_name, "vg", 2) == 0)  return VG_INDEX;
  62   if (strncasecmp(dorade_field_name, "vt", 2) == 0)  return VT_INDEX;
  63   if (strncasecmp(dorade_field_name, "ncp", 2) == 0) return NP_INDEX;
  64
  65   return -1;
  66 }
  67
  68 void prt_skipped_field_msg(char *dorade_field_name)
  69 {
  70   char prtname[9];
  71   int i, already_printed;
  72 #define MAXFIELDS 20
  73   static int nskipped = 0;
  74   static char skipped_list[MAXFIELDS][9];
  75
  76   /* Make sure name is a properly formed string. */
  77   strncpy(prtname, dorade_field_name, 8);
  78   prtname[8] = '\0';
  79
  80   /* We don't want to repeat messages for the same field, so keep a list of
  81    * fields already printed.
  82    */
  83   already_printed = 0;
  84   i = 0;
  85   while (!already_printed && i < nskipped) {
  86     if (strncmp(prtname, skipped_list[i], 2) == 0) already_printed = 1;
  87     i++;
  88   }
  89   if (!already_printed) {
  90     fprintf(stderr, "Unknown DORADE parameter type <%s> -- skipping.\n", prtname);
  91     strcpy(skipped_list[nskipped], prtname);
  92     nskipped++;
  93   }
  94 }
  95
  96 /* For date conversion routines. */
  97 static int daytab[2][13] = {
  98   {0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334, 365},
  99   {0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335, 366}
 100 };
 101
 102 /*************************************************************/
 103 /*                                                           */
 104 /*                          julian                           */
 105 /*                                                           */
 106 /*************************************************************/
 107 static int julian(int year, int mo, int day)
 108 {
 109 /* Converts a calendar date (month, day, year) to a Julian date.
 110          Returns:
 111            Julian day.
 112 */
 113   int leap;
 114
 115   leap = (year%4 == 0 && year%100 != 0) || year%400 == 0;
 116   return(day + daytab[leap][mo-1]);
 117 }
 118
 119 /*************************************************************/
 120 /*                                                           */
 121 /*                           ymd                             */
 122 /*                                                           */
 123 /*************************************************************/
 124 static void ymd(int jday, int year, int *mm, int *dd)
 125 {
 126   /* Input: jday, yyyy */
 127   /* Output: mm, dd */
 128   /* Copied from hdf_to_radar.c, written by Mike Kolander. */
 129
 130   int leap;
 131   int i;
 132
 133   leap = (year%4 == 0 && year%100 != 0) || year%400 == 0;
 134   for (i=0; daytab[leap][i]<jday; i++) continue;
 135   *mm = i;
 136   i--;
 137   *dd = jday - daytab[leap][i];
 138 }
 139
 140 /* Secretly defined in uf_to_radar.c */
 141 Volume *copy_sweeps_into_volume(Volume *new_volume, Volume *old_volume);
 142
 143 /**********************************************************************/
 144 /*                                                                    */
 145 /*                       RSL_dorade_to_radar                          */
 146 /*                                                                    */
 147 /**********************************************************************/
 148 Radar *RSL_dorade_to_radar(char *infile)
 149 {
 150   Radar  *radar;
 151   Volume *new_volume;
 152   Sweep  *sweep;
 153   Ray    *ray;
 154   int iv, iray, iparam;
 155   int range_bin1, gate_size;
 156
 157   int nbins, data_len, word_size;
 158   short *ptr2bytes;
 159   int   *ptr4bytes;
 160   float scale, offset, value;
 161   int datum, missing_data_flag;
 162   float prf;
 163
 164   float (*f)(Range x);
 165   Range (*invf)(float x);
 166
 167   FILE  *fp;
 168   Comment_block   *cb;
 169   Volume_desc     *vd;
 170   Sensor_desc    **sd;
 171   Sweep_record    *sr;
 172   Radar_desc      *rd;
 173   Data_ray        *dray;
 174   Parameter_data  *pd;
 175   Ray_info        *ray_info;
 176   Platform_info   *platform_info;
 177
 178   int nsweep;
 179   int i;
 180   char buf[1024];
 181
 182   int degree, minute;
 183   float second;
 184
 185   int year, month, day, jday, jday_vol;
 186
 187   radar = NULL;
 188   if (infile == NULL) {
 189     int save_fd;
 190     save_fd = dup(0);
 191     fp = fdopen(save_fd, "r");
 192   }  else
 193     if((fp=fopen(infile, "r"))==(FILE *)NULL) {
 194       perror(infile);
 195       return radar;
 196     }
 197
 198   fp = uncompress_pipe(fp); /* Transparently, use gunzip. */
 199
 200   cb = dorade_read_comment_block(fp);
 201
 202   /**********************************************************************/
 203
 204   vd = dorade_read_volume_desc(fp);   /* R E A D */
 205   if (radar_verbose_flag)   dorade_print_volume_desc(vd);  /* P R I N T */
 206
 207   /* R E A D */
 208   sd = (Sensor_desc **) calloc(vd->nsensors, sizeof(Sensor_desc *));
 209   for (i=0; i<vd->nsensors; i++) {
 210     sd[i] = dorade_read_sensor(fp);
 211   }
 212
 213   /* P R I N T */
 214   if (radar_verbose_flag) {
 215     for (i=0; i<vd->nsensors; i++) {
 216       fprintf(stderr, "============ S E N S O R   # %d =====================\n", i);
 217       dorade_print_sensor(sd[i]);
 218     }
 219   }
 220   /* R E A D   sweeps. */
 221   if (vd->nsensors > 1) {
 222     fprintf(stderr, "RSL_dorade_to_radar: Unable to process for more than 1 sensor.\n");
 223     fprintf(stderr, "RSL_dorade_to_radar: Number of sensors is %d\n", vd->nsensors);
 224     return NULL;
 225   }
 226
 227   /* Use sensor 0 for vitals. */
 228   rd = sd[0]->radar_desc;
 229   range_bin1 = sd[0]->cell_range_vector->range_cell[0];
 230   gate_size =  sd[0]->cell_range_vector->range_cell[1] - range_bin1;
 231
 232   radar = RSL_new_radar(MAX_RADAR_VOLUMES);
 233   radar->h.month = vd->month;
 234   radar->h.day   = vd->day;
 235   radar->h.year  = vd->year;
 236   radar->h.hour  = vd->hour;
 237   radar->h.minute = vd->minute;
 238   radar->h.sec    = vd->second;
 239   sprintf(radar->h.radar_type, "dorade");
 240   radar->h.number = 0;
 241   strncpy(radar->h.name, vd->flight_num, sizeof(radar->h.name));
 242   strncpy(radar->h.radar_name, rd->radar_name, sizeof(radar->h.radar_name));
 243   strncpy(radar->h.project, vd->project_name, sizeof(radar->h.project));
 244   sprintf(radar->h.city, "Unknown");
 245   strncpy(radar->h.state, "UKN", 3);
 246   sprintf(radar->h.country, "Unknown");
 247   /* Convert lat to d:m:s */
 248   degree = (int)rd->latitude;
 249   minute = (int)((rd->latitude - degree) * 60);
 250   second = (rd->latitude - degree - minute/60.0) * 3600.0;
 251   radar->h.latd = degree;
 252   radar->h.latm = minute;
 253   radar->h.lats = second;
 254   /* Convert lat to d:m:s */
 255   degree = (int)rd->longitude;
 256   minute = (int)((rd->longitude - degree) * 60);
 257   second = (rd->longitude - degree - minute/60.0) * 3600.0;
 258   radar->h.lond = degree;
 259   radar->h.lonm = minute;
 260   radar->h.lons = second;
 261   radar->h.height = rd->altitude * 1000.0;
 262   radar->h.spulse = 0; /* FIXME */
 263   radar->h.lpulse = 0; /* FIXME */
 264
 265   year = vd->year;
 266   month = vd->month;
 267   day = vd->day;
 268   jday_vol = julian(year, month, day);
 269
 270   /* TODO:
 271      Get any threshold values present in parameter descriptors.
 272
 273      Note: The parameter descriptor contains an 8-character string which may
 274      contain the, "Name of parameter upon which this parameter is thresholded".
 275      According to documentation, if there is no thresholding parameter, the
 276      string will be "NONE".  In practice, this string has occasionally been
 277      seen to contain only spaces for this condition.  The corresponding
 278      missing-data-flag for threshold may also vary, the nominal value being
 279      -999, but has also been -32768.
 280
 281      Pseudo code:
 282      nparam = rd->nparam_desc;  [<--useable code, i.e., not pseudo]
 283      create string array for threshold parameter names, size nparams.
 284      for each param:
 285        strcpy thresh param name from param desc to threshold param array.
 286        if string is all blanks, replace with 'NONE'.
 287      endfor
 288    */
 289
 290   /* Begin volume code. */
 291   /* We don't know how many sweeps per volume exist, until we read
 292    * the file.  So allocate a large number of pointers, hope we don't
 293    * exceed it, and adjust the pointer array at the end.  This is
 294    * efficient because we'll be manipulating pointers to the sweeps and
 295    * not the sweeps themselves.
 296    */
 297
 298   if (radar_verbose_flag)
 299     fprintf(stderr, "Number of parameters: %d\n", rd->nparam_desc);
 300
 301   /* All the parameters are together, however, their order within
 302    * the ray is not guarenteed.  For instance, VE could appear after
 303    * DM.  For this we'll keep a list of parameter names and perform
 304    * a (linear) search.  The result will be an index into the RSL
 305    * volume array (radar->v[i]).  It is likely that the order will be
 306    * consistant within a file, therefore, we'll keep track of the index of
 307    * our previous parameter type and begin the search from there; the next
 308    * index should be a match.
 309    *
 310    * The dorade parameter names and the rsl mapping is:
 311    *
 312    * Dorade: VE, DM, SW, DBZ, ZDR, PHI, RHO, LDR,  DX,  CH,  AH,  CV,  AV
 313    *    RSL: VR, DM, SW, DZ,  ZD,  PH,  RH,  LR,  *DX, *CH, *AH, *CV, *AV.
 314    *
 315    *    * means this is a new RSL name.
 316    */
 317
 318 #define DORADE_MAX_SWEEP 20
 319   nsweep = 0;
 320   while((sr = dorade_read_sweep(fp, sd))) {
 321     for(iray = 0; iray < sr->nrays; iray++) {
 322       dray = sr->data_ray[iray];
 323
 324       /* Now, loop through the parameters and fill the rsl structures. */
 325       for (iparam = 0; iparam < dray->nparam; iparam++) {
 326         pd = dray->parameter_data[iparam];
 327         iv = find_rsl_field_index(pd->name);
 328         if (iv < 0) {
 329           prt_skipped_field_msg(pd->name);
 330           continue;
 331         }
 332         if (radar->v[iv] == NULL) {
 333           radar->v[iv] = RSL_new_volume(DORADE_MAX_SWEEP); /* Expandable */
 334         } else if (nsweep >= radar->v[iv]->h.nsweeps) {
 335           /* Must expand the number of sweeps. */
 336           /* Expand by another DORADE_MAX_SWEEP. */
 337           if (radar_verbose_flag) {
 338             fprintf(stderr, "nsweeps (%d) exceeds radar->v[%d]->h.nsweeps (%d)."
 339               "\n", nsweep, iv, radar->v[iv]->h.nsweeps);
 340             fprintf(stderr, "Increasing it to %d sweeps\n",
 341               radar->v[iv]->h.nsweeps+DORADE_MAX_SWEEP);
 342           }
 343           new_volume = RSL_new_volume(radar->v[iv]->h.nsweeps+DORADE_MAX_SWEEP);
 344           /* Look in uf_to_radar.c for 'copy_sweeps_into_volume' */
 345           new_volume = copy_sweeps_into_volume(new_volume, radar->v[iv]);
 346           radar->v[iv] = new_volume;
 347         }
 348         /* Assign f and invf
 349         switch (iv) . . .
 350          * Or just use RSL_f_list[iv] and RSL_invf_list[iv] as in sweep.h below.
 351          */
 352
 353         /* Allocate the ray and load the parameter data. */
 354         if ((sweep = radar->v[iv]->sweep[nsweep]) == NULL) {
 355           sweep = radar->v[iv]->sweep[nsweep] = RSL_new_sweep(sr->s_info->nrays);
 356           sweep->h.sweep_num    = sr->s_info->sweep_num;
 357           sweep->h.elev         = sr->s_info->fixed_angle;
 358           sweep->h.beam_width   = rd->horizontal_beam_width;
 359           sweep->h.vert_half_bw = rd->vertical_beam_width / 2.0;
 360           /*
 361           sweep->h.vert_half_bw = radar->v[iv]->sweep[nsweep]->h.beam_width / 2.0;
 362           */
 363           sweep->h.horz_half_bw = rd->horizontal_beam_width / 2.0;
 364           sweep->h.f = RSL_f_list[iv];
 365           sweep->h.invf = RSL_invf_list[iv];
 366         }
 367
 368
 369         data_len = dray->data_len[iparam];
 370         word_size = dray->word_size[iparam];
 371         if (word_size != 2 && word_size != 4) {
 372           fprintf(stderr,"RSL_dorade_to_radar: dray->word_size[%d] = %d\n",
 373                   iparam, dray->word_size[iparam]);
 374           fprintf(stderr,"Expected value is 2 or 4.\n");
 375           return NULL;
 376         }
 377         nbins = data_len / word_size;
 378
 379         if ((ray = sweep->ray[iray]) == NULL) {
 380           if (radar_verbose_flag)
 381             fprintf(stderr, "Allocating %d bins for ray %d\n",
 382               dray->data_len[iparam], iray);
 383           ray = sweep->ray[iray] = RSL_new_ray(nbins);
 384         }
 385
 386         /* Load ray header. */
 387
 388         ray_info = dray->ray_info;
 389         jday = ray_info->jday;
 390         if (jday != jday_vol) {
 391           if (jday > 0 && jday < 367) {
 392             /* Recompute year month day */
 393             if (jday < jday_vol) year++;
 394             ymd(jday, year, &month, &day);
 395             jday_vol = jday;
 396           }
 397           else { /* Invalid jday */
 398           }
 399         }
 400         ray->h.year     = year;
 401         ray->h.month    = month;
 402         ray->h.day      = day;
 403         ray->h.hour     = ray_info->hour;
 404         ray->h.minute   = ray_info->minute;
 405         ray->h.sec      = ray_info->second + ray_info->msec / 1000.;
 406         ray->h.azimuth  = ray_info->azimuth;
 407         ray->h.ray_num  = iray + 1;
 408         ray->h.elev     = ray_info->elevation;
 409         ray->h.elev_num = ray_info->sweep_num;
 410         ray->h.unam_rng = rd->unambiguous_range;
 411         ray->h.nyq_vel = rd->unambiguous_velocity;
 412         ray->h.azim_rate = ray_info->scan_rate;
 413         /* TODO
 414         ray->h.fix_angle = ;
 415         */
 416         ray->h.range_bin1 = range_bin1;
 417         ray->h.gate_size = gate_size;
 418         platform_info = dray->platform_info;
 419         ray->h.pitch = platform_info->pitch;
 420         ray->h.roll = platform_info->roll;
 421         ray->h.heading = platform_info->heading;
 422         ray->h.pitch_rate = platform_info->pitch_rate;
 423         ray->h.heading_rate = platform_info->heading_rate;
 424         ray->h.lat = platform_info->latitude;
 425         ray->h.lon = platform_info->longitude;
 426         ray->h.alt = platform_info->altitude;
 427         ray->h.vel_east = platform_info->ew_speed;
 428         ray->h.vel_north = platform_info->ns_speed;
 429         ray->h.vel_up = platform_info->v_speed;
 430         /* TODO: Does DORADE have Doppler velocity res?
 431                  A: No.
 432               ray->h.vel_res = N/A
 433         */
 434         ray->h.beam_width = rd->horizontal_beam_width;
 435         /* TODO
 436         ray->h.frequency = @Radar Descriptor
 437           Get parm_desc.xmit_freq
 438           Find first set-bit, use frequency from corresponding index in rad_desc frequenies.
 439
 440         ray->h.pulse_count =  Is it number samples used? @Param desc.  Probably not.
 441
 442         * The following are DONE *
 443             ray->h.pulse_width = @Parameter Descriptor--not same--it's in meters--
 444                                     we use micro-sec. They're using pulse length.
 445                                     pulse width = pulse length/c
 446           ray->h.wavelength = Can compute using Nyquist velocity and PRF or PRT:
 447               wl = 4*vmax/PRF or wl = 4*vmax*PRT
 448           ray->h.prf = Can compute if nothing else: prf = c/(2*Rmax)
 449         */
 450         /* DORADE is actually using pulse length.  Convert to pulse width. */
 451         ray->h.pulse_width = (sd[0]->p_desc[iparam]->pulse_width/
 452             RSL_SPEED_OF_LIGHT) * 1e6;
 453         prf = RSL_SPEED_OF_LIGHT/(2.*rd->unambiguous_range*1000.);
 454         ray->h.prf = prf;
 455         ray->h.wavelength = (4.*rd->unambiguous_velocity)/prf;
 456         ray->h.nbins = nbins;
 457         ray->h.f = RSL_f_list[iv];
 458         ray->h.invf = RSL_invf_list[iv];
 459
 460         /* Copy the ray data into the RSL ray. */
 461
 462         /* .... fill here .... */
 463
 464         /* Assign pointer to data.
 465          * Get datum using word size and proper cast.
 466          * Convert and store in rsl ray->range.
 467          * Increment pointer to data based on word size.
 468          */
 469
 470         ptr2bytes = (short *) pd->data;
 471         ptr4bytes = (int *)   pd->data;
 472         scale  = sd[0]->p_desc[iparam]->scale_factor;
 473         offset = sd[0]->p_desc[iparam]->offset_factor;
 474         missing_data_flag = sd[0]->p_desc[iparam]->missing_data_flag;
 475
 476         for (i=0; i<nbins; i++) {
 477           if (word_size == 2) datum = *ptr2bytes++;
 478           else datum = *ptr4bytes++;
 479           /*
 480             TODO: If there is a threshold parameter for this parameter then
 481             apply threshold value.  I think threshold works like this: If there's a
 482             threshold parameter, as with VT (threshold parm = NCP), then use
 483             threshold value from that parameter unless it is the missing data value.
 484           */
 485           if (datum != missing_data_flag) {
 486             value = ((float)datum - offset) / scale;
 487           }
 488           else value = BADVAL;
 489
 490           ray->range[i] = ray->h.invf(value);
 491         }
 492
 493         if (iray == 0) {
 494           radar->v[iv]->h.nsweeps = nsweep + 1;
 495           radar->v[iv]->h.f = RSL_f_list[iv];
 496           radar->v[iv]->h.invf = RSL_invf_list[iv];
 497         }
 498       }
 499     }
 500     nsweep++;
 501     if (radar_verbose_flag) fprintf(stderr, "______NEW SWEEP__<%d>____\n", nsweep);
 502     /* Save for loading into volume structure. */
 503     dorade_free_sweep(sr);
 504   }
 505
 506   /* The following avoids a broken pipe message, since a VOLD at the end
 507    * is not read yet.
 508    */
 509   while(fread(buf, sizeof(buf), 1, fp)) continue;  /* Read til EOF */
 510
 511   rsl_pclose(fp);
 512
 513   return radar;
 514 }