Pileus Git - ~andy/rsl/blob - src/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   FILE  *fp;
 165   Comment_block   *cb;
 166   Volume_desc     *vd;
 167   Sensor_desc    **sd;
 168   Sweep_record    *sr;
 169   Radar_desc      *rd;
 170   Data_ray        *dray;
 171   Parameter_data  *pd;
 172   Ray_info        *ray_info;
 173   Platform_info   *platform_info;
 174
 175   int nsweep;
 176   int i;
 177   char buf[1024];
 178
 179   int degree, minute;
 180   float second;
 181
 182   int year, month, day, jday, jday_vol;
 183
 184   radar = NULL;
 185   if (infile == NULL) {
 186     int save_fd;
 187     save_fd = dup(0);
 188     fp = fdopen(save_fd, "r");
 189   }  else
 190     if((fp=fopen(infile, "r"))==(FILE *)NULL) {
 191       perror(infile);
 192       return radar;
 193     }
 194
 195   fp = uncompress_pipe(fp); /* Transparently, use gunzip. */
 196
 197   cb = dorade_read_comment_block(fp);
 198
 199   /**********************************************************************/
 200
 201   vd = dorade_read_volume_desc(fp);   /* R E A D */
 202   if (radar_verbose_flag)   dorade_print_volume_desc(vd);  /* P R I N T */
 203
 204   /* R E A D */
 205   sd = (Sensor_desc **) calloc(vd->nsensors, sizeof(Sensor_desc *));
 206   for (i=0; i<vd->nsensors; i++) {
 207     sd[i] = dorade_read_sensor(fp);
 208   }
 209
 210   /* P R I N T */
 211   if (radar_verbose_flag) {
 212     for (i=0; i<vd->nsensors; i++) {
 213       fprintf(stderr, "============ S E N S O R   # %d =====================\n", i);
 214       dorade_print_sensor(sd[i]);
 215     }
 216   }
 217   /* R E A D   sweeps. */
 218   if (vd->nsensors > 1) {
 219     fprintf(stderr, "RSL_dorade_to_radar: Unable to process for more than 1 sensor.\n");
 220     fprintf(stderr, "RSL_dorade_to_radar: Number of sensors is %d\n", vd->nsensors);
 221     return NULL;
 222   }
 223
 224   /* Use sensor 0 for vitals. */
 225   rd = sd[0]->radar_desc;
 226   range_bin1 = sd[0]->cell_range_vector->range_cell[0];
 227   gate_size =  sd[0]->cell_range_vector->range_cell[1] - range_bin1;
 228
 229   radar = RSL_new_radar(MAX_RADAR_VOLUMES);
 230   radar->h.month = vd->month;
 231   radar->h.day   = vd->day;
 232   radar->h.year  = vd->year;
 233   radar->h.hour  = vd->hour;
 234   radar->h.minute = vd->minute;
 235   radar->h.sec    = vd->second;
 236   sprintf(radar->h.radar_type, "dorade");
 237   radar->h.number = 0;
 238   strncpy(radar->h.name, vd->flight_num, sizeof(radar->h.name));
 239   strncpy(radar->h.radar_name, rd->radar_name, sizeof(radar->h.radar_name));
 240   strncpy(radar->h.project, vd->project_name, sizeof(radar->h.project));
 241   sprintf(radar->h.city, "Unknown");
 242   strncpy(radar->h.state, "UKN", 3);
 243   sprintf(radar->h.country, "Unknown");
 244   /* Convert lat to d:m:s */
 245   degree = (int)rd->latitude;
 246   minute = (int)((rd->latitude - degree) * 60);
 247   second = (rd->latitude - degree - minute/60.0) * 3600.0;
 248   radar->h.latd = degree;
 249   radar->h.latm = minute;
 250   radar->h.lats = second;
 251   /* Convert lat to d:m:s */
 252   degree = (int)rd->longitude;
 253   minute = (int)((rd->longitude - degree) * 60);
 254   second = (rd->longitude - degree - minute/60.0) * 3600.0;
 255   radar->h.lond = degree;
 256   radar->h.lonm = minute;
 257   radar->h.lons = second;
 258   radar->h.height = rd->altitude * 1000.0;
 259   radar->h.spulse = 0; /* FIXME */
 260   radar->h.lpulse = 0; /* FIXME */
 261
 262   year = vd->year;
 263   month = vd->month;
 264   day = vd->day;
 265   jday_vol = julian(year, month, day);
 266
 267   /* TODO:
 268      Get any threshold values present in parameter descriptors.
 269
 270      Note: The parameter descriptor contains an 8-character string which may
 271      contain the, "Name of parameter upon which this parameter is thresholded".
 272      According to documentation, if there is no thresholding parameter, the
 273      string will be "NONE".  In practice, this string has occasionally been
 274      seen to contain only spaces for this condition.  The corresponding
 275      missing-data-flag for threshold may also vary, the nominal value being
 276      -999, but has also been -32768.
 277
 278      Pseudo code:
 279      nparam = rd->nparam_desc;  [<--useable code, i.e., not pseudo]
 280      create string array for threshold parameter names, size nparams.
 281      for each param:
 282        strcpy thresh param name from param desc to threshold param array.
 283        if string is all blanks, replace with 'NONE'.
 284      endfor
 285    */
 286
 287   /* Begin volume code. */
 288   /* We don't know how many sweeps per volume exist, until we read
 289    * the file.  So allocate a large number of pointers, hope we don't
 290    * exceed it, and adjust the pointer array at the end.  This is
 291    * efficient because we'll be manipulating pointers to the sweeps and
 292    * not the sweeps themselves.
 293    */
 294
 295   if (radar_verbose_flag)
 296     fprintf(stderr, "Number of parameters: %d\n", rd->nparam_desc);
 297
 298   /* All the parameters are together, however, their order within
 299    * the ray is not guarenteed.  For instance, VE could appear after
 300    * DM.  For this we'll keep a list of parameter names and perform
 301    * a (linear) search.  The result will be an index into the RSL
 302    * volume array (radar->v[i]).  It is likely that the order will be
 303    * consistant within a file, therefore, we'll keep track of the index of
 304    * our previous parameter type and begin the search from there; the next
 305    * index should be a match.
 306    *
 307    * The dorade parameter names and the rsl mapping is:
 308    *
 309    * Dorade: VE, DM, SW, DBZ, ZDR, PHI, RHO, LDR,  DX,  CH,  AH,  CV,  AV
 310    *    RSL: VR, DM, SW, DZ,  ZD,  PH,  RH,  LR,  *DX, *CH, *AH, *CV, *AV.
 311    *
 312    *    * means this is a new RSL name.
 313    */
 314
 315 #define DORADE_MAX_SWEEP 20
 316   nsweep = 0;
 317   while((sr = dorade_read_sweep(fp, sd))) {
 318     for(iray = 0; iray < sr->nrays; iray++) {
 319       dray = sr->data_ray[iray];
 320
 321       /* Now, loop through the parameters and fill the rsl structures. */
 322       for (iparam = 0; iparam < dray->nparam; iparam++) {
 323         pd = dray->parameter_data[iparam];
 324         iv = find_rsl_field_index(pd->name);
 325         if (iv < 0) {
 326           prt_skipped_field_msg(pd->name);
 327           continue;
 328         }
 329         if (radar->v[iv] == NULL) {
 330           radar->v[iv] = RSL_new_volume(DORADE_MAX_SWEEP); /* Expandable */
 331         } else if (nsweep >= radar->v[iv]->h.nsweeps) {
 332           /* Must expand the number of sweeps. */
 333           /* Expand by another DORADE_MAX_SWEEP. */
 334           if (radar_verbose_flag) {
 335             fprintf(stderr, "nsweeps (%d) exceeds radar->v[%d]->h.nsweeps (%d)."
 336               "\n", nsweep, iv, radar->v[iv]->h.nsweeps);
 337             fprintf(stderr, "Increasing it to %d sweeps\n",
 338               radar->v[iv]->h.nsweeps+DORADE_MAX_SWEEP);
 339           }
 340           new_volume = RSL_new_volume(radar->v[iv]->h.nsweeps+DORADE_MAX_SWEEP);
 341           /* Look in uf_to_radar.c for 'copy_sweeps_into_volume' */
 342           new_volume = copy_sweeps_into_volume(new_volume, radar->v[iv]);
 343           radar->v[iv] = new_volume;
 344         }
 345         /* Assign f and invf
 346         switch (iv) . . .
 347          * Or just use RSL_f_list[iv] and RSL_invf_list[iv] as in sweep.h below.
 348          */
 349
 350         /* Allocate the ray and load the parameter data. */
 351         if ((sweep = radar->v[iv]->sweep[nsweep]) == NULL) {
 352           sweep = radar->v[iv]->sweep[nsweep] = RSL_new_sweep(sr->s_info->nrays);
 353           sweep->h.sweep_num    = sr->s_info->sweep_num;
 354           sweep->h.elev         = sr->s_info->fixed_angle;
 355           sweep->h.beam_width   = rd->horizontal_beam_width;
 356           sweep->h.vert_half_bw = rd->vertical_beam_width / 2.0;
 357           /*
 358           sweep->h.vert_half_bw = radar->v[iv]->sweep[nsweep]->h.beam_width / 2.0;
 359           */
 360           sweep->h.horz_half_bw = rd->horizontal_beam_width / 2.0;
 361           sweep->h.f = RSL_f_list[iv];
 362           sweep->h.invf = RSL_invf_list[iv];
 363         }
 364
 365
 366         data_len = dray->data_len[iparam];
 367         word_size = dray->word_size[iparam];
 368         if (word_size != 2 && word_size != 4) {
 369           fprintf(stderr,"RSL_dorade_to_radar: dray->word_size[%d] = %d\n",
 370                   iparam, dray->word_size[iparam]);
 371           fprintf(stderr,"Expected value is 2 or 4.\n");
 372           return NULL;
 373         }
 374         nbins = data_len / word_size;
 375
 376         if ((ray = sweep->ray[iray]) == NULL) {
 377           if (radar_verbose_flag)
 378             fprintf(stderr, "Allocating %d bins for ray %d\n",
 379               dray->data_len[iparam], iray);
 380           ray = sweep->ray[iray] = RSL_new_ray(nbins);
 381         }
 382
 383         /* Load ray header. */
 384
 385         ray_info = dray->ray_info;
 386         jday = ray_info->jday;
 387         if (jday != jday_vol) {
 388           if (jday > 0 && jday < 367) {
 389             /* Recompute year month day */
 390             if (jday < jday_vol) year++;
 391             ymd(jday, year, &month, &day);
 392             jday_vol = jday;
 393           }
 394           else { /* Invalid jday */
 395           }
 396         }
 397         ray->h.year     = year;
 398         ray->h.month    = month;
 399         ray->h.day      = day;
 400         ray->h.hour     = ray_info->hour;
 401         ray->h.minute   = ray_info->minute;
 402         ray->h.sec      = ray_info->second + ray_info->msec / 1000.;
 403         ray->h.azimuth  = ray_info->azimuth;
 404         ray->h.ray_num  = iray + 1;
 405         ray->h.elev     = ray_info->elevation;
 406         ray->h.elev_num = ray_info->sweep_num;
 407         ray->h.unam_rng = rd->unambiguous_range;
 408         ray->h.nyq_vel = rd->unambiguous_velocity;
 409         ray->h.azim_rate = ray_info->scan_rate;
 410         /* TODO
 411         ray->h.fix_angle = ;
 412         */
 413         ray->h.range_bin1 = range_bin1;
 414         ray->h.gate_size = gate_size;
 415         platform_info = dray->platform_info;
 416         ray->h.pitch = platform_info->pitch;
 417         ray->h.roll = platform_info->roll;
 418         ray->h.heading = platform_info->heading;
 419         ray->h.pitch_rate = platform_info->pitch_rate;
 420         ray->h.heading_rate = platform_info->heading_rate;
 421         ray->h.lat = platform_info->latitude;
 422         ray->h.lon = platform_info->longitude;
 423         ray->h.alt = platform_info->altitude;
 424         ray->h.vel_east = platform_info->ew_speed;
 425         ray->h.vel_north = platform_info->ns_speed;
 426         ray->h.vel_up = platform_info->v_speed;
 427         /* TODO: Does DORADE have Doppler velocity res?
 428                  A: No.
 429               ray->h.vel_res = N/A
 430         */
 431         ray->h.beam_width = rd->horizontal_beam_width;
 432         /* TODO
 433         ray->h.frequency = @Radar Descriptor
 434           Get parm_desc.xmit_freq
 435           Find first set-bit, use frequency from corresponding index in rad_desc frequenies.
 436
 437         ray->h.pulse_count =  Is it number samples used? @Param desc.  Probably not.
 438
 439         * The following are DONE *
 440             ray->h.pulse_width = @Parameter Descriptor--not same--it's in meters--
 441                                     we use micro-sec. They're using pulse length.
 442                                     pulse width = pulse length/c
 443           ray->h.wavelength = Can compute using Nyquist velocity and PRF or PRT:
 444               wl = 4*vmax/PRF or wl = 4*vmax*PRT
 445           ray->h.prf = Can compute if nothing else: prf = c/(2*Rmax)
 446         */
 447         /* DORADE is actually using pulse length.  Convert to pulse width. */
 448         ray->h.pulse_width = (sd[0]->p_desc[iparam]->pulse_width/
 449             RSL_SPEED_OF_LIGHT) * 1e6;
 450         prf = RSL_SPEED_OF_LIGHT/(2.*rd->unambiguous_range*1000.);
 451         ray->h.prf = prf;
 452         ray->h.wavelength = (4.*rd->unambiguous_velocity)/prf;
 453         ray->h.nbins = nbins;
 454         ray->h.f = RSL_f_list[iv];
 455         ray->h.invf = RSL_invf_list[iv];
 456
 457         /* Copy the ray data into the RSL ray. */
 458
 459         /* .... fill here .... */
 460
 461         /* Assign pointer to data.
 462          * Get datum using word size and proper cast.
 463          * Convert and store in rsl ray->range.
 464          * Increment pointer to data based on word size.
 465          */
 466
 467         ptr2bytes = (short *) pd->data;
 468         ptr4bytes = (int *)   pd->data;
 469         scale  = sd[0]->p_desc[iparam]->scale_factor;
 470         offset = sd[0]->p_desc[iparam]->offset_factor;
 471         missing_data_flag = sd[0]->p_desc[iparam]->missing_data_flag;
 472
 473         for (i=0; i<nbins; i++) {
 474           if (word_size == 2) datum = *ptr2bytes++;
 475           else datum = *ptr4bytes++;
 476           /*
 477             TODO: If there is a threshold parameter for this parameter then
 478             apply threshold value.  I think threshold works like this: If there's a
 479             threshold parameter, as with VT (threshold parm = NCP), then use
 480             threshold value from that parameter unless it is the missing data value.
 481           */
 482           if (datum != missing_data_flag) {
 483             value = ((float)datum - offset) / scale;
 484           }
 485           else value = BADVAL;
 486
 487           ray->range[i] = ray->h.invf(value);
 488         }
 489
 490         if (iray == 0) {
 491           radar->v[iv]->h.nsweeps = nsweep + 1;
 492           radar->v[iv]->h.f = RSL_f_list[iv];
 493           radar->v[iv]->h.invf = RSL_invf_list[iv];
 494         }
 495       }
 496     }
 497     nsweep++;
 498     if (radar_verbose_flag) fprintf(stderr, "______NEW SWEEP__<%d>____\n", nsweep);
 499     /* Save for loading into volume structure. */
 500     dorade_free_sweep(sr);
 501   }
 502
 503   /* The following avoids a broken pipe message, since a VOLD at the end
 504    * is not read yet.
 505    */
 506   while(fread(buf, sizeof(buf), 1, fp)) continue;  /* Read til EOF */
 507
 508   rsl_pclose(fp);
 509
 510   return radar;
 511 }