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