Pileus Git - ~andy/rsl/blob - hdf_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 #ifdef HAVE_CONFIG_H
  24 #include "config.h"
  25 #endif
  26 #ifdef HAVE_LIBTSDISTK
  27 /******************************************************************
  28    Reads one volume scan from a HDF file into a RSL radar structure.
  29
  30   -----------------------------------------------------------------
  31          Libraries required for execution of this code :
  32       -ltsdistk                    : tsdis toolkit
  33       -lmfhdf -ldf -ljpeg -lz      : HDF
  34       -lrsl                        : rsl
  35       -lm                          : C math
  36
  37   -----------------------------------------------------------------
  38 *******************************************************************/
  39
  40
  41 #include <math.h>
  42 #include <time.h>
  43 #include <string.h>
  44 #include <stdio.h>
  45 #include <stdlib.h>
  46 #include <unistd.h>
  47 /* TSDIS toolkit function and structure definitions. */
  48 #include "IO.h"
  49 #include "IO_GV.h"
  50 /* RSL function and structure definitions. */
  51 #include "rsl.h"
  52 /* Parameter definitions for 1B-51 and 1C-51 HDF
  53          file handling applications using the TSDIS toolkit. */
  54 #include "toolkit_1BC-51_appl.h"
  55
  56 #define MISSING_VAL 0
  57
  58
  59 /*************************************************************/
  60 /*                                                           */
  61 /*                Function Prototypes                        */
  62 /*                                                           */
  63 /*************************************************************/
  64 void RayFillFrom1B51(Ray *ray, int16 *rayData, PARAMETER_DESCRIPTOR *parmDesc);
  65 void RayFillFrom1C51(Ray *ray, int vindex, int16 *rayData, int8 *rayMaskData,
  66                                                                                  PARAMETER_DESCRIPTOR *parmDesc, float calibr);
  67 static void Ray_headerFill(Ray *ray, L1B_1C_GV *gvl1, VosSize *vs,
  68                                                                                                          int pindex, int tk_sindex, int rindex);
  69 Ray *RayBuild(L1B_1C_GV *gvl1, VosSize *vs, float calibr,
  70                                                         int vindex, int pindex, int sindex, int rindex);
  71 static void Sweep_headerFill(Sweep *sweep, SENSORS *sensor, int sindex, int nrays);
  72 Sweep *SweepBuild(L1B_1C_GV *gvl1, VosSize *vs, float calibr,
  73                                                                         int vindex, int pindex, int sindex);
  74 static void Volume_headerFill(Volume *volume, char *parmDesc, int vindex,
  75                                                                                          int nsweeps, float calibr);
  76 Volume *VolumeBuild(L1B_1C_GV *gvl1, VosSize *vs, float calibr,
  77                                                                                 int vindex, int pindex);
  78 int parmIdentify(char *parmName);
  79 static void Radar_headerFill(Radar *radar, L1B_1C_GV *gvl1);
  80 Radar *RadarBuild(L1B_1C_GV *gvl1, VosSize *vs, float zCal);
  81 int commentsRead(VosSize *vs, float *zCal, char *comments, int productID);
  82 static L1B_1C_GV *GVL1Build(IO_HANDLE *granuleHandle, int vosNum,
  83                                                                                                                 VosSize *vs);
  84 int metaDataRead(Radar *radar, IO_HANDLE *granuleHandle);
  85 static int hdfFileOpen(char *infile, IO_HANDLE *granuleHandle,
  86                                                                                          char *hdfFileName, int *vosNum);
  87 Radar *RSL_hdf_to_radar(char *infile);
  88
  89
  90 /* Toolkit memory management functions. */
  91 extern void TKfreeGVL1(L1B_1C_GV *gvl1);
  92 extern int8 ***TKnewParmData1byte(int nsweep, int nray, int ncell);
  93 extern int16 ***TKnewParmData2byte(int nsweep, int nray, int ncell);
  94 extern PARAMETER *TKnewGVL1parm(void);
  95 extern L1B_1C_GV *TKnewGVL1(void);
  96
  97
  98 static float (*f)(Range x);
  99 static Range (*invf)(float x);
 100
 101 extern int radar_verbose_flag;
 102
 103
 104 static void ymd(int jday, int yy, int *mm, int *dd);
 105
 106 static int daytab[2][13] = {
 107   {0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334, 365},
 108   {0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335, 366}
 109 };
 110
 111 static void ymd(int jday, int year, int *mm, int *dd)
 112 {
 113   /*  Input: jday, yyyy */
 114   /* Output: mm, dd */
 115   int leap;
 116   int i;
 117
 118   leap = (year%4 == 0 && year%100 != 0) || year%400 == 0;
 119   for (i=0; daytab[leap][i]<jday; i++) continue;
 120   *mm = i;
 121   i--;
 122   *dd = jday - daytab[leap][i];
 123 }
 124
 125 /*************************************************************/
 126 /*                                                           */
 127 /*                       RayFillFrom1B51                     */
 128 /*                                                           */
 129 /*************************************************************/
 130 void RayFillFrom1B51(Ray *ray, int16 *rayData, PARAMETER_DESCRIPTOR *parmDesc)
 131 {
 132 /*
 133          Fill the RSL bin slots of one ray of any volume using the corresponding
 134          ray data from a 1B-51 HDF file.
 135 */
 136         int j;
 137
 138         for (j=0; j<ray->h.nbins; j++)
 139         {
 140                 if (rayData[j] <= AP_VALUE)  /* Handle anomalous condition flags. */
 141                 {
 142                         if (rayData[j] == NO_VALUE) ray->range[j] = invf((float)BADVAL);
 143                         else if (rayData[j] == RNG_AMBIG_VALUE) ray->range[j] = invf((float)RFVAL);
 144                         else if (rayData[j] == NOECHO_VALUE) ray->range[j] = invf((float)NOECHO);
 145                         else ray->range[j] = invf((float)APFLAG);
 146                 }
 147                 else  /* Valid data value */
 148                 {
 149                         ray->range[j] = invf( (rayData[j] - parmDesc->offsetFactor) /
 150                                                      parmDesc->scaleFactor );
 151                 }
 152         } /* end for (j=0... */
 153 }
 154
 155 /*************************************************************/
 156 /*                                                           */
 157 /*                         RayFillFrom1C51                   */
 158 /*                                                           */
 159 /*************************************************************/
 160 void RayFillFrom1C51(Ray *ray, int vindex, int16 *rayData, int8 *rayMaskData,
 161                                                                                  PARAMETER_DESCRIPTOR *parmDesc, float calibr)
 162 {
 163 /*
 164          Fill the RSL bin slots of one ray of a CZ or a DZ volume using the
 165          corresponding ray data and ray_mask data from a 1C-51 HDF file.
 166 */
 167         int j;
 168
 169         for (j=0; j<ray->h.nbins; j++)
 170         {
 171                 if (rayData[j] <= AP_VALUE)  /* Handle anomalous condition flags. */
 172                 {
 173                         if (rayData[j] == NO_VALUE) ray->range[j] = invf((float)BADVAL);
 174                         else if (rayData[j] == RNG_AMBIG_VALUE) ray->range[j] = invf((float)RFVAL);
 175                         else if (rayData[j] == NOECHO_VALUE) ray->range[j] = invf((float)NOECHO);
 176                         else ray->range[j] = invf((float)APFLAG);
 177                 }
 178                 else  /* Valid data value */
 179                 {
 180                         if ((vindex == CZ_INDEX) || (vindex == CD_INDEX))
 181                         {
 182                                 if (rayMaskData[j] == 1)
 183                                   ray->range[j] = invf((float)BADVAL);
 184                                 else
 185                                   ray->range[j] = invf( (rayData[j] - parmDesc->offsetFactor) /
 186                                          parmDesc->scaleFactor );
 187                         }  /* end if (vindex == CZ_INDEX) */
 188                         else if (vindex == DZ_INDEX)
 189                         {
 190                                 ray->range[j] = invf(
 191                                                         ((rayData[j] - parmDesc->offsetFactor) / parmDesc->scaleFactor) -
 192                                                                 calibr + X * rayMaskData[j] );
 193
 194                         }  /* end else if DZ_INDEX */
 195                         else if (vindex == ZD_INDEX)
 196                         {
 197                                 ray->range[j] = invf(
 198                                                 ((rayData[j] - parmDesc->offsetFactor) / parmDesc->scaleFactor) +
 199                                                         X * rayMaskData[j] );
 200                         }  /* end else if ZD_INDEX */
 201                         else
 202                           fprintf(stderr, "RayFillFrom1C51(): Illegal volume index..\n");
 203                 } /* else valid data value */
 204         }  /* end for (j=0; ... */
 205 }
 206
 207 /*************************************************************/
 208 /*                                                           */
 209 /*                      Ray_headerFill                       */
 210 /*                                                           */
 211 /*************************************************************/
 212 void Ray_headerFill(Ray *ray, L1B_1C_GV *gvl1, VosSize *vs,
 213                                                                                 int pindex, int tk_sindex, int rindex)
 214 {
 215         ray->h.year = (int)gvl1->volDes.year;
 216         /* Get calendar date (month, day) from (year, Julian day) */
 217         ymd((int)gvl1->sensor.rayInfoInteger[tk_sindex][rindex][1],
 218                 ray->h.year, &ray->h.month, &ray->h.day);
 219         ray->h.hour = (int)gvl1->sensor.rayInfoInteger[tk_sindex][rindex][2];
 220         ray->h.minute = (int)gvl1->sensor.rayInfoInteger[tk_sindex][rindex][3];
 221         ray->h.sec = (float)(gvl1->sensor.rayInfoInteger[tk_sindex][rindex][4] +
 222                                                                                          gvl1->sensor.rayInfoInteger[tk_sindex][rindex][5]/1000.0);
 223
 224         ray->h.azimuth = gvl1->sensor.rayInfoFloat[tk_sindex][rindex][0];
 225         ray->h.ray_num = rindex + 1;
 226         ray->h.elev = gvl1->sensor.rayInfoFloat[tk_sindex][rindex][1]; /* degrees */
 227         ray->h.elev_num = tk_sindex + 1;
 228         ray->h.gate_size = (int)
 229         (gvl1->sensor.parm[pindex]->cellRangeVector.distanceToCell[2] -
 230          gvl1->sensor.parm[pindex]->cellRangeVector.distanceToCell[1]); /*meters*/
 231
 232         ray->h.range_bin1 = (int)
 233                   (gvl1->sensor.parm[pindex]->cellRangeVector.distanceToCell[0] -
 234                                                 0.5 * ray->h.gate_size); /* meters */
 235
 236         ray->h.vel_res = MISSING_VAL;           /* ?? */
 237         /* Sweeps/min */
 238         ray->h.sweep_rate = (float)(gvl1->sensor.radarDesc.nomScanRate / 6.0);
 239         ray->h.prf = (int)gvl1->sensor.rayInfoFloat[tk_sindex][rindex][3];
 240         ray->h.azim_rate = (float)gvl1->sensor.radarDesc.nomScanRate;
 241         ray->h.fix_angle = (float)gvl1->sensor.sweepInfo[tk_sindex].fixedAngle;
 242         ray->h.pulse_count = (int)gvl1->sensor.rayInfoFloat[tk_sindex][rindex][2];
 243         /* Pulse width (microsec) */
 244         ray->h.pulse_width = (float)(gvl1->sensor.parm[pindex]->parmDesc.pulseWidth /
 245                                                                                                                          300.0);
 246         ray->h.beam_width = (float)gvl1->sensor.radarDesc.horBeamWidth;
 247         /* Carrier freq (GHz) */
 248         ray->h.frequency = (float)gvl1->sensor.radarDesc.frequency1;
 249         /* wavelength (m) */
 250         if (ray->h.frequency != 0.0)
 251           ray->h.wavelength = (RSL_SPEED_OF_LIGHT / ray->h.frequency) * 1.0e-9;
 252         else
 253           ray->h.wavelength = 0.0;
 254         ray->h.nyq_vel = (float) (ray->h.prf * ray->h.wavelength / 4.0);
 255         if (ray->h.prf != 0)
 256           ray->h.unam_rng = (float) RSL_SPEED_OF_LIGHT / (2.0 * ray->h.prf * 1000.0);
 257         else
 258           ray->h.unam_rng = (float) 0.0;
 259         ray->h.nbins = vs->tk.ncell[tk_sindex][pindex];
 260         ray->h.f = f;
 261         ray->h.invf = invf;
 262 }
 263
 264 /*************************************************************/
 265 /*                                                           */
 266 /*                          RayBuild                         */
 267 /*                                                           */
 268 /*************************************************************/
 269 Ray *RayBuild(L1B_1C_GV *gvl1, VosSize *vs, float calibr,
 270                                                         int vindex, int pindex, int tk_sindex, int rindex)
 271 {
 272         Ray *ray;
 273
 274         /* Create a Ray structure. */
 275         ray = RSL_new_ray(vs->tk.ncell[tk_sindex][pindex]);
 276         if (ray == NULL)
 277         {
 278                 perror("RayBuild(): RSL_new_ray failed\n");
 279                 return(NULL);
 280         }
 281         Ray_headerFill(ray, gvl1, vs, pindex, tk_sindex, rindex);
 282         /* Is this a 1C-51 file? */
 283         if ((strcmp(gvl1->sensor.parm[pindex]->parmDesc.parmName, "QCZ") == 0) ||
 284                         (strcmp(gvl1->sensor.parm[pindex]->parmDesc.parmName, "QCZDR") == 0))
 285           RayFillFrom1C51(ray, vindex,
 286                                                 gvl1->sensor.parm[pindex]->parmData2byte[tk_sindex][rindex],
 287                                                 gvl1->sensor.parm[pindex-1]->parmData1byte[tk_sindex][rindex],
 288                                                 &gvl1->sensor.parm[pindex]->parmDesc, calibr);
 289         else /* 1B-51 file */
 290           RayFillFrom1B51(ray,
 291                                                 gvl1->sensor.parm[pindex]->parmData2byte[tk_sindex][rindex],
 292                                                 &gvl1->sensor.parm[pindex]->parmDesc);
 293         return(ray);
 294 }
 295
 296 /*************************************************************/
 297 /*                                                           */
 298 /*                     Sweep_headerFill                      */
 299 /*                                                           */
 300 /*************************************************************/
 301 void Sweep_headerFill(Sweep *sweep, SENSORS *sensor, int tk_sindex, int nrays)
 302 {
 303 /*      sweep->h.sweep_num filled in VolumeBuild() */
 304         sweep->h.elev = sensor->sweepInfo[tk_sindex].fixedAngle;
 305         sweep->h.beam_width = sensor->radarDesc.horBeamWidth;
 306         sweep->h.horz_half_bw = sensor->radarDesc.horBeamWidth / 2.0;
 307         sweep->h.vert_half_bw = sensor->radarDesc.verBeamWidth / 2.0;
 308         sweep->h.nrays = sensor->sweepInfo[tk_sindex].numRays;
 309         sweep->h.f = f;
 310         sweep->h.invf = invf;
 311 }
 312
 313 /*************************************************************/
 314 /*                                                           */
 315 /*                        SweepBuild                         */
 316 /*                                                           */
 317 /*************************************************************/
 318 Sweep *SweepBuild(L1B_1C_GV *gvl1, VosSize *vs, float calibr,
 319                                                                         int vindex, int pindex, int tk_sindex)
 320 {
 321         int rindex;
 322         Sweep *sweep;
 323
 324         /* Create a Sweep structure. */
 325         sweep = RSL_new_sweep(vs->rsl.maxNray);
 326         if (sweep == NULL)
 327         {
 328                 perror("SweepBuild(): RSL_new_sweep failed\n");
 329                 return(NULL);
 330         }
 331         /* Initialize the Sweep_header values. */
 332         Sweep_headerFill(sweep, &gvl1->sensor, tk_sindex, vs->rsl.maxNray);
 333
 334         /* Loop to fill each of the rays of this rsl sweep structure. */
 335         for (rindex=0; rindex<vs->tk.nray[tk_sindex]; rindex++)
 336                 sweep->ray[rindex] = RayBuild(gvl1, vs, calibr, vindex, pindex,
 337                                                                                                                                         tk_sindex, rindex);
 338         return(sweep);
 339 }
 340
 341 /*************************************************************/
 342 /*                                                           */
 343 /*                    Volume_headerFill                      */
 344 /*                                                           */
 345 /*************************************************************/
 346 void Volume_headerFill(Volume *volume, char *parmDesc, int vindex,
 347                                                                                          int nsweeps, float calibr)
 348 {
 349         if (vindex == DZ_INDEX)
 350           volume->h.type_str = strdup("Reflectivity");
 351         else if (vindex == ZD_INDEX)
 352           volume->h.type_str = strdup("Differential Reflectivity");
 353         else volume->h.type_str = strdup(parmDesc);
 354         volume->h.f = f;
 355         volume->h.invf = invf;
 356         volume->h.nsweeps = nsweeps;
 357         volume->h.calibr_const = calibr;
 358 }
 359
 360 /*************************************************************/
 361 /*                                                           */
 362 /*                        VolumeBuild                        */
 363 /*                                                           */
 364 /*************************************************************/
 365 Volume *VolumeBuild(L1B_1C_GV *gvl1, VosSize *vs, float calibr,
 366                                                                                 int vindex, int pindex)
 367 {
 368         Volume *v;
 369         int sindex, tk_sindex;
 370         extern int *rsl_qsweep; /* See RSL_read_these_sweeps in volume.c */
 371         extern int rsl_qsweep_max;
 372
 373         /* Create a Volume structure. */
 374         v = RSL_new_volume(vs->tk.nsweep);
 375         if (v == NULL)
 376         {
 377                 perror("VolumeBuild(): RSL_new_volume failed\n");
 378                 return(NULL);
 379         }
 380         /* Initialize the Volume_header values. */
 381         Volume_headerFill(v, gvl1->sensor.parm[pindex]->parmDesc.parmDesc,
 382                                                                                 vindex, vs->tk.nsweep, calibr);
 383         if (radar_verbose_flag)
 384           fprintf(stderr, "RSL volume type: %s\n", v->h.type_str);
 385
 386         /* Build each of the sweeps of this radar volume structure. */
 387         sindex = -1;
 388         for (tk_sindex=0; tk_sindex<vs->tk.nsweep; tk_sindex++)
 389         {
 390           if (rsl_qsweep != NULL) {
 391                 if (tk_sindex > rsl_qsweep_max) break;
 392                 if (rsl_qsweep[tk_sindex] == 0) continue;
 393           }
 394           /* If data for this parm type exists in this toolkit sweep,
 395                    then move it into a rsl sweep. */
 396           if (vs->tk.ncell[tk_sindex][pindex] > 0)
 397                 {
 398                         sindex++;
 399                   v->sweep[sindex] = SweepBuild(gvl1, vs, calibr, vindex, pindex,
 400                                                                                                                                                 tk_sindex);
 401                         v->sweep[sindex]->h.sweep_num = sindex + 1;
 402                         if (radar_verbose_flag)
 403                           fprintf(stderr, "  rsl_sweep[%02d]  elev=%4.1f  nrays=%d  cells/ray=%d\n",
 404                                                                 v->sweep[sindex]->h.sweep_num-1, v->sweep[sindex]->h.elev,
 405                                                                 vs->tk.nray[tk_sindex], vs->tk.ncell[tk_sindex][pindex]);
 406                 }
 407         } /* end for (tk_sindex=0;...*/
 408         return(v);
 409 }
 410
 411 /*************************************************************/
 412 /*                                                           */
 413 /*                       parmIdentify                        */
 414 /*                                                           */
 415 /*************************************************************/
 416 int parmIdentify(char *parmName)
 417 /* Identify the parameter type stored in the L1B_1C_GV structure.
 418          Upon success, return the corresponding RSL radar volume XX_INDEX
 419          value.
 420          Upon failure, return -1 .
 421 */
 422 {
 423         int vindex;
 424
 425         if (strcmp(parmName, "Z") == 0)
 426         {
 427                 vindex = DZ_INDEX;   invf = DZ_INVF;   f = DZ_F;
 428         }
 429         else if (strcmp(parmName, "V") == 0)
 430         {
 431                 vindex = VR_INDEX;   invf = VR_INVF;   f = VR_F;
 432         }
 433         else if (strcmp(parmName, "QCZ") == 0)
 434         {
 435                 vindex = CZ_INDEX;   invf = CZ_INVF;   f = CZ_F;
 436         }
 437         else if (strcmp(parmName, "ZDR") == 0)
 438         {
 439                 vindex = ZD_INDEX;   invf = ZD_INVF;   f = ZD_F;
 440         }
 441         else if (strcmp(parmName, "QCZDR") == 0)
 442         {
 443                 vindex = CD_INDEX;   invf = CD_INVF;   f = CD_F;
 444         }
 445         else if (strcmp(parmName, "QCMZ") == 0)
 446         {
 447                 vindex = MZ_INDEX;   invf = MZ_INVF;   f = MZ_F;
 448         }
 449         else if (strcmp(parmName, "QCMZDR") == 0)
 450         {
 451                 vindex = MD_INDEX;   invf = MD_INVF;   f = MD_F;
 452         }
 453         else  /* Unknown */
 454         {
 455                 return(-1);
 456         }
 457
 458         return(vindex);
 459 }
 460
 461 /*************************************************************/
 462 /*                                                           */
 463 /*                    Radar_headerFill                       */
 464 /*                                                           */
 465 /*************************************************************/
 466 void Radar_headerFill(Radar *radar, L1B_1C_GV *gvl1)
 467 {
 468         double x;
 469
 470         radar->h.month = (int)gvl1->volDes.month;
 471         radar->h.day = (int)gvl1->volDes.day;
 472         radar->h.year = (int)gvl1->volDes.year;
 473         radar->h.hour = (int)gvl1->volDes.hour;
 474         radar->h.minute = (int)gvl1->volDes.minute;
 475         radar->h.sec = (float)gvl1->volDes.second;
 476         strncpy(radar->h.radar_type, "**", 48);       /*********/
 477         radar->h.nvolumes = MAX_RADAR_VOLUMES;
 478         radar->h.number = MISSING_VAL;
 479         strncpy(radar->h.name, gvl1->sensor.radarDesc.radarName, 7);
 480         strncpy(radar->h.radar_name, gvl1->sensor.sweepInfo[0].radarName, 7);
 481
 482         /* Radar Latitude */
 483         x = fabs(gvl1->sensor.radarDesc.radarLat);
 484         radar->h.latd = (int)floor(x);
 485         x = (x - radar->h.latd) * 60.0;
 486         radar->h.latm = (int)floor(x);
 487         x = (x - radar->h.latm) * 60.0;
 488         radar->h.lats = (int)floor(x + 0.5);  /* round up */
 489         if (gvl1->sensor.radarDesc.radarLat < 0)
 490         {
 491           radar->h.latd = -radar->h.latd;
 492                 radar->h.latm = -radar->h.latm;
 493                 radar->h.lats = -radar->h.lats;
 494         }
 495
 496         /* Radar Longitude */
 497         x = fabs(gvl1->sensor.radarDesc.radarLon);
 498         radar->h.lond = (int)floor(x);
 499         x = (x - radar->h.lond) * 60.0;
 500         radar->h.lonm = (int)floor(x);
 501         x = (x - radar->h.lonm) * 60.0;
 502         radar->h.lons = (int)floor(x + 0.5);  /* round up */
 503         if (gvl1->sensor.radarDesc.radarLon < 0)
 504         {
 505           radar->h.lond = -radar->h.lond;
 506                 radar->h.lonm = -radar->h.lonm;
 507                 radar->h.lons = -radar->h.lons;
 508         }
 509
 510         radar->h.height = (int)(1000.0 * gvl1->sensor.radarDesc.radarAlt + 0.5);
 511         radar->h.spulse = MISSING_VAL;  /* ns */
 512         radar->h.lpulse = MISSING_VAL;  /* ns */
 513 }
 514
 515 /*************************************************************/
 516 /*                                                           */
 517 /*                        RadarBuild                         */
 518 /*                                                           */
 519 /*************************************************************/
 520 Radar *RadarBuild(L1B_1C_GV *gvl1, VosSize *vs, float zCal)
 521 /* Creates and fills a RSL radar structure with data obtained
 522          from the L1B_1C_GV structure.
 523
 524          If success, returns a pointer to the radar structure.
 525          If failure, returns NULL.
 526 */
 527 {
 528         Radar *radar;
 529         extern int rsl_qfield[];
 530         int pindex, vindex;
 531
 532         if (radar_verbose_flag)
 533         {
 534           fprintf(stderr, "\n****** Moving VOS from toolkit L1GV structure -> RSL structure ...\n");
 535         }
 536         /* Create a structure of type Radar */
 537         radar = (Radar *)RSL_new_radar(MAX_RADAR_VOLUMES);
 538         if (radar == NULL)
 539         {
 540                 perror("RadarBuild(): Error creating radar structure.\n");
 541                 return(NULL);
 542         }
 543         /* Initialize the Radar_header values. */
 544         Radar_headerFill(radar, gvl1);
 545
 546         /* Build each of the 'nparm' volumes of the radar structure. */
 547         for (pindex=0; pindex<vs->tk.nparm; pindex++)
 548         {
 549                 /* Identify parameter type, so we know which RSL volume to load the data
 550                          into. */
 551                 vindex = parmIdentify(gvl1->sensor.parm[pindex]->parmDesc.parmName);
 552                 if (vindex < 0)
 553                 {
 554                         fprintf(stderr,
 555                                   "RadarBuild(): Unexpected parameter type: %s found in HDF file.\n",
 556                                                         gvl1->sensor.parm[pindex]->parmDesc.parmName);
 557                 }
 558                 /* Don't build mask volumes. */
 559                 else if ((vindex == MZ_INDEX) || (vindex == MD_INDEX)) continue;
 560                 else if (rsl_qfield[vindex] == 0) /* Don't build unselected volumes. */
 561                 {
 562                   if (radar_verbose_flag)
 563                         {
 564                           fprintf(stderr, "Field %s not selected for retrieval from HDF file.\n",
 565                                                                 gvl1->sensor.parm[pindex]->parmDesc.parmName);
 566                                 if (vindex == CZ_INDEX)
 567                             fprintf(stderr, "Field 'DZ' unselected for retrieval from 1C-51 file.\n");
 568                                 else if (vindex == CD_INDEX)
 569                             fprintf(stderr, "Field 'ZD' unselected for retrieval from 1C-51 file.\n");
 570                         }
 571                 }  /* end else if (rsl_qfield[vindex] == 0) */
 572                 else if (vindex == CZ_INDEX)  /* Handle CZ and DZ volumes. */
 573                 {
 574                         /* Build the RSL CZ volume. */
 575                         radar->v[vindex] = VolumeBuild(gvl1, vs, zCal, vindex, pindex);
 576                         /* If required, build a RSL DZ volume. */
 577                         if (rsl_qfield[DZ_INDEX])
 578                         {
 579                                 if (radar_verbose_flag)
 580                             fprintf(stderr, "Constructing reflectivity volume 'DZ'\n");
 581                                 radar->v[DZ_INDEX] = VolumeBuild(gvl1, vs, zCal, DZ_INDEX, pindex);
 582                         }
 583                 }  /* end if (vindex == CZ_INDEX) */
 584                 else if (vindex == CD_INDEX)  /* Handle CD and ZD volumes. */
 585                 {
 586                         /* Build the RSL CD volume. */
 587                         radar->v[vindex] = VolumeBuild(gvl1, vs, 0.0, vindex, pindex);
 588                         /* If required, build a RSL ZD volume. */
 589                         if (rsl_qfield[ZD_INDEX])
 590                         {
 591                                 if (radar_verbose_flag)
 592                             fprintf(stderr, "Constructing reflectivity volume 'ZD'\n");
 593                                 radar->v[ZD_INDEX] = VolumeBuild(gvl1, vs, 0.0, ZD_INDEX, pindex);
 594                         }
 595                 }  /* end if (vindex == CD_INDEX) */
 596                 else   /* Handle all 1B-51 fields. (DZ, ZD, VR) */
 597                 {
 598                   radar->v[vindex] = VolumeBuild(gvl1, vs, 0.0, vindex, pindex);
 599                 }
 600         }  /* end for (pindex=0; ...) */
 601
 602         return(radar);
 603 }
 604
 605 /*************************************************************/
 606 /*                                                           */
 607 /*                       commentsRead                        */
 608 /*                                                           */
 609 /*************************************************************/
 610 int commentsRead(VosSize *vs, float *zCal, char *comments, int productID)
 611 {
 612 /* Parse the comments field of the 'L1B_1C_GV' structure.
 613          Retrieve the number_of_cells/ray values for each parameter, and
 614          store in the 'VosSize' structure.
 615
 616          Returns: OK if success.
 617                   <0 if failure.
 618 */
 619         char *spointer;
 620         char record[2][2048];  /* 2 records is maximum possible. */
 621         char parseString[1024];
 622         int nrecords, pindex, tk_sindex;
 623         float qcParm[NUMBER_QC_PARAMS];
 624
 625         /* Construct a format string to read the records in the comments
 626                  field. A logical record here is the block of ascii characters
 627                  which details the toolkit dimensions of one VOS.
 628                  For a 1B-51 file, there should be 1 such record.
 629                  For a 1C-51 file, there is one additional record for QC parms.*/
 630
 631         strcpy(parseString, "");
 632         strcat(parseString, "%[^*] %*[*\n]");
 633         nrecords = 1;
 634         if (productID == TK_L1C_GV)
 635         {
 636           strcat(parseString, "%[^\n]");
 637                 nrecords++;
 638         }
 639         /* Read all records from the comments field into the record buffers. */
 640         if (sscanf(comments, parseString, record[0], record[1]) != nrecords)
 641           goto quit;
 642
 643
 644         if (sscanf(record[0], "nSweep=%d", &vs->tk.nsweep) != 1) goto quit;
 645
 646         strcpy(parseString, "nRay=%d\n");
 647         for (pindex=0; pindex<vs->tk.nparm; pindex++)
 648     strcat(parseString, "nCell_parm[%*d]=%d\n");
 649
 650         spointer = record[0];
 651         for (tk_sindex=0; tk_sindex<vs->tk.nsweep; tk_sindex++)
 652         {
 653                 spointer = strstr(spointer, "nRay=");
 654                 if (sscanf(spointer, parseString, &vs->tk.nray[tk_sindex],
 655                                                          &vs->tk.ncell[tk_sindex][0], &vs->tk.ncell[tk_sindex][1],
 656                                                          &vs->tk.ncell[tk_sindex][2], &vs->tk.ncell[tk_sindex][3])
 657                                 != vs->tk.nparm+1) goto quit;
 658                 spointer = spointer + 5;
 659         }
 660
 661
 662         /* If 1C-51 file, read the QC parameters into the qcParm array. */
 663         if (productID == TK_L1C_GV)
 664         {
 665                 if (sscanf(record[1], "-hThresh1 %f -hThresh2 %f -hThresh3 %f -zThresh0 %f -zThresh1 %f -zThresh2 %f -zThresh3 %f -hFreeze %f -dbzNoise %f -zCal %f",
 666                                                          &qcParm[HTHRESH1], &qcParm[HTHRESH2], &qcParm[HTHRESH3],
 667                                                          &qcParm[ZTHRESH0], &qcParm[ZTHRESH1], &qcParm[ZTHRESH2], &qcParm[ZTHRESH3],
 668                                                          &qcParm[HFREEZE], &qcParm[DBZNOISE], &qcParm[ZCAL]) != NUMBER_QC_PARAMS)
 669                   goto quit;
 670                 /* Print out the QC parameters we've just read in. */
 671                 /*
 672                 if (radar_verbose_flag)
 673                 {
 674                         fprintf(stderr, "\n****** Reading VOS QC Parameters from HDF file...\n");
 675                         fprintf(stderr, "hThresh1: %.2f   hThresh2: %.2f   hThresh3: %.2f\n",
 676                                                         qcParm[HTHRESH1], qcParm[HTHRESH2], qcParm[HTHRESH3]);
 677                         fprintf(stderr, "zThresh0: %.2f  zThresh1: %.2f  zThresh2: %.2f  zThresh3: %.2f\n",
 678                                                         qcParm[ZTHRESH0], qcParm[ZTHRESH1], qcParm[ZTHRESH2], qcParm[ZTHRESH3]);
 679                         fprintf(stderr, "hFreeze: %.2f  dbzNoise: %.2f   zCal: %.2f\n\n",
 680                                                         qcParm[HFREEZE], qcParm[DBZNOISE], qcParm[ZCAL]);
 681                 }
 682                 */
 683                 if (qcParm[ZCAL] <= NOVAL_FLOAT) *zCal = 0.0;
 684                 else *zCal = qcParm[ZCAL];
 685         } /* end if (productID == TK_L1C_GV) */
 686
 687         return(OK);
 688
 689  quit:
 690         if (radar_verbose_flag)
 691           fprintf(stderr, "commentsRead(): Failure reading comments field\n");
 692         return(ABORT);
 693 }
 694
 695 /*************************************************************/
 696 /*                                                           */
 697 /*                          GVL1Build                        */
 698 /*                                                           */
 699 /*************************************************************/
 700 static L1B_1C_GV *GVL1Build(IO_HANDLE *granuleHandle, int vosNum,
 701                                                                                                                 VosSize *vs)
 702 {
 703 /* Build a toolkit 'L1B_1C_GV' structure sized for the VOS we
 704          will later read in from the HDF file.
 705          Returns:
 706            gvl1 if success.
 707                  NULL if fails.
 708 */
 709         int ncell, pindex;
 710         L1B_1C_GV *gvl1;
 711
 712         /* Using the toolkit, get the toolkit VOS dimensions from
 713                  the HDF file. Note that the toolkit dimensions are distinct
 714                  from the RSL VOS dimensions. */
 715         vs->tk.nparm = TKgetNparm(granuleHandle, vosNum);
 716 /*  TK_FAIL is now defined as 1?????? */
 717         if (vs->tk.nparm <= 0)
 718         {
 719                 fprintf(stderr, "GVL1Build(): TKgetNparm() failed.\n");
 720                 return(NULL);
 721         }
 722         vs->tk.nsweep = TKgetNsweep(granuleHandle, vosNum);
 723         if (vs->tk.nsweep <= 0)
 724         {
 725                 fprintf(stderr, "GVL1Build(): TKgetNsweep() failed.\n");
 726                 return(NULL);
 727         }
 728         vs->rsl.maxNray = TKgetNray(granuleHandle, vosNum);
 729         if (vs->rsl.maxNray <= 0)
 730         {
 731                 fprintf(stderr, "GVL1Build(): TKgetNray() failed.\n");
 732                 return(NULL);
 733         }
 734
 735         /* Allocate memory for a TSDIS 'L1B_1C_GV' structure. */
 736         gvl1 = (L1B_1C_GV *)TKnewGVL1();
 737         for (pindex=0; pindex<vs->tk.nparm; pindex++)
 738         {
 739                 /* Allocate memory for a new parameter. */
 740                 gvl1->sensor.parm[pindex] = (PARAMETER *)TKnewGVL1parm();
 741                 ncell = TKgetNcell(granuleHandle, vosNum, pindex);
 742                 /* Allocate memory for a 3D array to contain mask values and/or data. */
 743                 if (granuleHandle->productID == TK_L1B_GV)
 744                 {
 745                   gvl1->sensor.parm[pindex]->parmData2byte =
 746                             (int16 ***)TKnewParmData2byte(vs->tk.nsweep, vs->rsl.maxNray, ncell);
 747                 }
 748                 else   /* 1C-51 */
 749                 {
 750                         /* Odd parameters contain data, even parameters contain masks. */
 751                         if ((pindex % 2) == 0)  /* Mask? */
 752                     gvl1->sensor.parm[pindex]->parmData1byte =
 753                                    (int8 ***)TKnewParmData1byte(vs->tk.nsweep, vs->rsl.maxNray, ncell);
 754                         else  /* data */
 755                     gvl1->sensor.parm[pindex]->parmData2byte =
 756                                    (int16 ***)TKnewParmData2byte(vs->tk.nsweep, vs->rsl.maxNray, ncell);
 757                 }
 758         } /* end for (pindex=0; ... */
 759
 760         return(gvl1);
 761 }
 762
 763 /*************************************************************/
 764 /*                                                           */
 765 /*                   metaDataRead                            */
 766 /*                                                           */
 767 /*************************************************************/
 768 int metaDataRead(Radar *radar, IO_HANDLE *granuleHandle)
 769 {
 770         char buf[64];
 771
 772         TKreadMetadataChar(granuleHandle, TK_RADAR_CITY, buf);
 773         strncpy(radar->h.city, buf, 14);
 774         TKreadMetadataChar(granuleHandle, TK_RADAR_STATE, buf);
 775         strncpy(radar->h.state, buf, 2);
 776         return(OK);
 777 }
 778
 779 /*************************************************************/
 780 /*                                                           */
 781 /*                        hdfFileOpen                        */
 782 /*                                                           */
 783 /*************************************************************/
 784 static int hdfFileOpen(char *infile, IO_HANDLE *granuleHandle,
 785                                                                                          char *hdfFileName, int *vosNum)
 786 {
 787 /* Opens, if necessary, an HDF file. Checks that a VOS, not previously
 788          retrieved, exists in the HDF file.
 789
 790    Returns:
 791            OK, if success.
 792                  <0, if failure.
 793 */
 794         char *product;
 795         int productID, nvos, status;
 796
 797         /* If we presently have an open HDF file, check that it is the
 798            file requested; ie, 'infile'. If it's not, we first close
 799                  the open HDF file. */
 800         if (*vosNum != 0)
 801         {
 802                 if (strcmp(hdfFileName, infile) != 0)
 803                 {
 804                         if (TKclose(granuleHandle) != TK_SUCCESS)
 805                         {
 806                           fprintf(stderr, "hdfFileOpen(): *** TKclose() error\n");
 807                                 return(ABORT);
 808                         }
 809                         *vosNum = 0;
 810                 }
 811         } /* end if (*vosNum != 0) */
 812
 813         /* If first VOS of HDF file, we need first to open the file. */
 814         if (*vosNum == 0)
 815         {
 816                 strncpy(hdfFileName, infile, TK_MAX_FILENAME-1);
 817                 /* Get the desired product out of the HDF filename. */
 818                 product = strrchr(hdfFileName, '/');
 819                 if (product == NULL)
 820                   product = hdfFileName;
 821                 else
 822                   product = (char *)(product + 1);
 823                 if (strncmp(product, "1B51", 4) == 0)
 824                   productID = TK_L1B_GV;
 825                 else if (strncmp(product, "1C51", 4) == 0)
 826                   productID = TK_L1C_GV;
 827                 else
 828                 {
 829                         fprintf(stderr, "hdfFileOpen(): Invalid HDF filename.\n");
 830                         return(ABORT);
 831                 }
 832                 status = TKopen(hdfFileName, productID, TK_READ_ONLY, granuleHandle);
 833                 if (status != TK_SUCCESS)
 834                 {
 835                         fprintf(stderr, "hdfFileOpen(): *** TKopen() error\n");
 836                         return(ABORT);
 837                 }
 838         }  /* end if (*vosNum == 0) */
 839
 840         /* Check if the requested VOS exists in the HDF granule. */
 841         nvos = (int)TKgetNvos(granuleHandle);
 842         if (nvos == 0)
 843         {
 844                 if (radar_verbose_flag)
 845                   fprintf(stderr, "\nEmpty granule.\n");
 846                 return(QUIT);
 847         }
 848         else if (*vosNum+1 > nvos)
 849         {
 850                 if (radar_verbose_flag)
 851                   fprintf(stderr, "\nAll VOSs read from HDF file: %s\n", hdfFileName);
 852                 return(QUIT);
 853         }
 854         else if (nvos < 0)
 855         {
 856                 fprintf(stderr, "hdfFileOpen():*** TKgetNvos() error\n");
 857                 return(QUIT);
 858         }
 859
 860         return(OK);
 861 }
 862
 863 /*************************************************************/
 864 /*                                                           */
 865 /*                      RSL_hdf_to_radar                     */
 866 /*                                                           */
 867 /*************************************************************/
 868 Radar *RSL_hdf_to_radar(char *infile)
 869 {
 870         /* Reads one volume scan from a HDF file into a RSL radar structure.
 871                  It is envisioned that all VOSs will normally be retrieved,
 872                  one after the other, from an HDF file. Therefore, in the
 873                  interest of efficiency, this function is designed to keep
 874                  an HDF file open between VOS retrievals.
 875
 876                  Returns:
 877                    - pointer to a filled radar structure, if success.
 878                          - NULL pointer, if failure.
 879         */
 880         int pindex, status;
 881         float zCal=0.0;         /* Z Calibration constant. */
 882         VosSize vs;             /* VOS dimensions storage. */
 883         Radar *radar;           /* RSL structure for VOS storage. */
 884         L1B_1C_GV *gvl1;    /* TSDIS structure for VOS storage. */
 885         /* Following values must remain between invocations of this function. */
 886         static IO_HANDLE granuleHandle;
 887         static char hdfFileName[TK_MAX_FILENAME];
 888         static int vosNum=0;  /* No. of last VOS read from HDF file. (1...N) */
 889
 890         /* Open, if necessary, an HDF file. */
 891         status = hdfFileOpen(infile, &granuleHandle, hdfFileName, &vosNum);
 892         if (status < 0) goto quit;
 893
 894         /* Initialize the 'VosSize' structure. */
 895         memset(&vs, '\0', sizeof(VosSize));
 896
 897         /* Build a toolkit 'L1B_1C_GV' structure correctly sized for the VOS we
 898          will next read in from the HDF file. */
 899         gvl1 = (L1B_1C_GV *)GVL1Build(&granuleHandle, vosNum, &vs);
 900         if (gvl1 == NULL) goto quit;
 901
 902         /* Read VOS from HDF file into the toolkit L1B_1C_GV structure. */
 903         if (radar_verbose_flag)
 904           fprintf(stderr, "\n\n***** Moving VOS from HDF file -> toolkit L1GV structure ...\n");
 905         status = TKreadL1GV(&granuleHandle, gvl1);
 906         if (status != TK_SUCCESS)
 907         {
 908                 fprintf(stderr, "RSL_hdf_to_radar(): *** TKreadL1GV() error\n");
 909                 goto quit;
 910         }
 911
 912         if (radar_verbose_flag)
 913         {
 914                 fprintf(stderr, "Input file: %s\n", hdfFileName);
 915                 fprintf(stderr, "VOS date:   %.2d/%.2d/%d\n", gvl1->volDes.month,
 916                                                 gvl1->volDes.day, gvl1->volDes.year);
 917                 fprintf(stderr, "VOS time:   %.2d:%.2d:%.2d\n", gvl1->volDes.hour,
 918                                                 gvl1->volDes.minute, gvl1->volDes.second);
 919
 920                 fprintf(stderr, "Granule VOS #: %d\n", gvl1->volDes.volNum);
 921                 fprintf(stderr, "VOS Fields:\n");
 922                 for (pindex=0; pindex<vs.tk.nparm ; pindex++)
 923                   fprintf(stderr, "  %d: %s\n", pindex+1,
 924                                                         gvl1->sensor.parm[pindex]->parmDesc.parmDesc);
 925         }
 926
 927         /* Scan thru the comments field of the 'gvl1' structure for
 928                  the number_of_cells/ray/sweep for each parameter, and store
 929                  the values in vs, so we can next build a correctly sized RSL
 930                  structure to contain the VOS.
 931         */
 932         if (commentsRead(&vs, &zCal, gvl1->comments, granuleHandle.productID) < 0)
 933           goto quit;
 934
 935         /* Move VOS from L1B_1C_GV structure to radar structure. */
 936         radar = (Radar *)RadarBuild(gvl1, &vs, zCal);
 937         /* There are a couple metadata items needed for insertion into the
 938                  radar->header. */
 939         status = metaDataRead(radar, &granuleHandle);
 940
 941         /* Free memory allocated to the toolkit GVL1 structure. */
 942         TKfreeGVL1(gvl1);
 943         if (radar == NULL) goto quit;
 944
 945         vosNum++;
 946         return(radar);
 947
 948  quit:
 949         if (status == QUIT)
 950           TKclose(&granuleHandle);
 951         return(NULL);
 952 }
 953
 954 #else
 955 /*
 956  * Just declare and return something when we're told we don't have
 957  * TSDISTK nor HDF.  Do this because RSL_anyformat_to_radar references
 958  * this routine; linking won't fail because of no HDF.
 959  */
 960
 961 #include "rsl.h"
 962 Radar *RSL_hdf_to_radar(char *infile)
 963 {
 964   return NULL;
 965 }
 966 #endif