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