3 This is the TRMM Office Radar Software Library.
4 Copyright (C) 1996, 1997
6 Space Applications Corporation
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.
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.
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.
26 #ifdef HAVE_LIBTSDISTK
27 /******************************************************************
28 Reads one volume scan from a HDF file into a RSL radar structure.
30 -----------------------------------------------------------------
31 Libraries required for execution of this code :
32 -ltsdistk : tsdis toolkit
33 -lmfhdf -ldf -ljpeg -lz : HDF
37 -----------------------------------------------------------------
38 *******************************************************************/
47 /* TSDIS toolkit function and structure definitions. */
50 /* RSL function and structure definitions. */
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"
59 /*************************************************************/
61 /* Function Prototypes */
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,
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);
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);
98 static float (*f)(Range x);
99 static Range (*invf)(float x);
101 extern int radar_verbose_flag;
104 static void ymd(int jday, int yy, int *mm, int *dd);
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}
111 static void ymd(int jday, int year, int *mm, int *dd)
113 /* Input: jday, yyyy */
118 leap = (year%4 == 0 && year%100 != 0) || year%400 == 0;
119 for (i=0; daytab[leap][i]<jday; i++) continue;
122 *dd = jday - daytab[leap][i];
125 /*************************************************************/
127 /* RayFillFrom1B51 */
129 /*************************************************************/
130 void RayFillFrom1B51(Ray *ray, int16 *rayData, PARAMETER_DESCRIPTOR *parmDesc)
133 Fill the RSL bin slots of one ray of any volume using the corresponding
134 ray data from a 1B-51 HDF file.
138 for (j=0; j<ray->h.nbins; j++)
140 if (rayData[j] <= AP_VALUE) /* Handle anomalous condition flags. */
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);
147 else /* Valid data value */
149 ray->range[j] = invf( (rayData[j] - parmDesc->offsetFactor) /
150 parmDesc->scaleFactor );
152 } /* end for (j=0... */
155 /*************************************************************/
157 /* RayFillFrom1C51 */
159 /*************************************************************/
160 void RayFillFrom1C51(Ray *ray, int vindex, int16 *rayData, int8 *rayMaskData,
161 PARAMETER_DESCRIPTOR *parmDesc, float calibr)
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.
169 for (j=0; j<ray->h.nbins; j++)
171 if (rayData[j] <= AP_VALUE) /* Handle anomalous condition flags. */
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);
178 else /* Valid data value */
180 if ((vindex == CZ_INDEX) || (vindex == CD_INDEX))
182 if (rayMaskData[j] == 1)
183 ray->range[j] = invf((float)BADVAL);
185 ray->range[j] = invf( (rayData[j] - parmDesc->offsetFactor) /
186 parmDesc->scaleFactor );
187 } /* end if (vindex == CZ_INDEX) */
188 else if (vindex == DZ_INDEX)
190 ray->range[j] = invf(
191 ((rayData[j] - parmDesc->offsetFactor) / parmDesc->scaleFactor) -
192 calibr + X * rayMaskData[j] );
194 } /* end else if DZ_INDEX */
195 else if (vindex == ZD_INDEX)
197 ray->range[j] = invf(
198 ((rayData[j] - parmDesc->offsetFactor) / parmDesc->scaleFactor) +
199 X * rayMaskData[j] );
200 } /* end else if ZD_INDEX */
202 fprintf(stderr, "RayFillFrom1C51(): Illegal volume index..\n");
203 } /* else valid data value */
204 } /* end for (j=0; ... */
207 /*************************************************************/
211 /*************************************************************/
212 void Ray_headerFill(Ray *ray, L1B_1C_GV *gvl1, VosSize *vs,
213 int pindex, int tk_sindex, int rindex)
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);
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*/
232 ray->h.range_bin1 = (int)
233 (gvl1->sensor.parm[pindex]->cellRangeVector.distanceToCell[0] -
234 0.5 * ray->h.gate_size); /* meters */
236 ray->h.vel_res = MISSING_VAL; /* ?? */
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 /
246 ray->h.beam_width = (float)gvl1->sensor.radarDesc.horBeamWidth;
247 /* Carrier freq (GHz) */
248 ray->h.frequency = (float)gvl1->sensor.radarDesc.frequency1;
250 if (ray->h.frequency != 0.0)
251 ray->h.wavelength = (RSL_SPEED_OF_LIGHT / ray->h.frequency) * 1.0e-9;
253 ray->h.wavelength = 0.0;
254 ray->h.nyq_vel = (float) (ray->h.prf * ray->h.wavelength / 4.0);
256 ray->h.unam_rng = (float) RSL_SPEED_OF_LIGHT / (2.0 * ray->h.prf * 1000.0);
258 ray->h.unam_rng = (float) 0.0;
259 ray->h.nbins = vs->tk.ncell[tk_sindex][pindex];
264 /*************************************************************/
268 /*************************************************************/
269 Ray *RayBuild(L1B_1C_GV *gvl1, VosSize *vs, float calibr,
270 int vindex, int pindex, int tk_sindex, int rindex)
274 /* Create a Ray structure. */
275 ray = RSL_new_ray(vs->tk.ncell[tk_sindex][pindex]);
278 perror("RayBuild(): RSL_new_ray failed\n");
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 */
291 gvl1->sensor.parm[pindex]->parmData2byte[tk_sindex][rindex],
292 &gvl1->sensor.parm[pindex]->parmDesc);
296 /*************************************************************/
298 /* Sweep_headerFill */
300 /*************************************************************/
301 void Sweep_headerFill(Sweep *sweep, SENSORS *sensor, int tk_sindex, int nrays)
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;
310 sweep->h.invf = invf;
313 /*************************************************************/
317 /*************************************************************/
318 Sweep *SweepBuild(L1B_1C_GV *gvl1, VosSize *vs, float calibr,
319 int vindex, int pindex, int tk_sindex)
324 /* Create a Sweep structure. */
325 sweep = RSL_new_sweep(vs->rsl.maxNray);
328 perror("SweepBuild(): RSL_new_sweep failed\n");
331 /* Initialize the Sweep_header values. */
332 Sweep_headerFill(sweep, &gvl1->sensor, tk_sindex, vs->rsl.maxNray);
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,
341 /*************************************************************/
343 /* Volume_headerFill */
345 /*************************************************************/
346 void Volume_headerFill(Volume *volume, char *parmDesc, int vindex,
347 int nsweeps, float calibr)
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);
355 volume->h.invf = invf;
356 volume->h.nsweeps = nsweeps;
357 volume->h.calibr_const = calibr;
360 /*************************************************************/
364 /*************************************************************/
365 Volume *VolumeBuild(L1B_1C_GV *gvl1, VosSize *vs, float calibr,
366 int vindex, int pindex)
369 int sindex, tk_sindex;
370 extern int *rsl_qsweep; /* See RSL_read_these_sweeps in volume.c */
371 extern int rsl_qsweep_max;
373 /* Create a Volume structure. */
374 v = RSL_new_volume(vs->tk.nsweep);
377 perror("VolumeBuild(): RSL_new_volume failed\n");
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);
386 /* Build each of the sweeps of this radar volume structure. */
388 for (tk_sindex=0; tk_sindex<vs->tk.nsweep; tk_sindex++)
390 if (rsl_qsweep != NULL) {
391 if (tk_sindex > rsl_qsweep_max) break;
392 if (rsl_qsweep[tk_sindex] == 0) continue;
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)
399 v->sweep[sindex] = SweepBuild(gvl1, vs, calibr, vindex, pindex,
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]);
407 } /* end for (tk_sindex=0;...*/
411 /*************************************************************/
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
420 Upon failure, return -1 .
425 if (strcmp(parmName, "Z") == 0)
427 vindex = DZ_INDEX; invf = DZ_INVF; f = DZ_F;
429 else if (strcmp(parmName, "V") == 0)
431 vindex = VR_INDEX; invf = VR_INVF; f = VR_F;
433 else if (strcmp(parmName, "QCZ") == 0)
435 vindex = CZ_INDEX; invf = CZ_INVF; f = CZ_F;
437 else if (strcmp(parmName, "ZDR") == 0)
439 vindex = ZD_INDEX; invf = ZD_INVF; f = ZD_F;
441 else if (strcmp(parmName, "QCZDR") == 0)
443 vindex = CD_INDEX; invf = CD_INVF; f = CD_F;
445 else if (strcmp(parmName, "QCMZ") == 0)
447 vindex = MZ_INDEX; invf = MZ_INVF; f = MZ_F;
449 else if (strcmp(parmName, "QCMZDR") == 0)
451 vindex = MD_INDEX; invf = MD_INVF; f = MD_F;
461 /*************************************************************/
463 /* Radar_headerFill */
465 /*************************************************************/
466 void Radar_headerFill(Radar *radar, L1B_1C_GV *gvl1)
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);
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)
491 radar->h.latd = -radar->h.latd;
492 radar->h.latm = -radar->h.latm;
493 radar->h.lats = -radar->h.lats;
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)
505 radar->h.lond = -radar->h.lond;
506 radar->h.lonm = -radar->h.lonm;
507 radar->h.lons = -radar->h.lons;
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 */
515 /*************************************************************/
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.
524 If success, returns a pointer to the radar structure.
525 If failure, returns NULL.
529 extern int rsl_qfield[];
532 if (radar_verbose_flag)
534 fprintf(stderr, "\n****** Moving VOS from toolkit L1GV structure -> RSL structure ...\n");
536 /* Create a structure of type Radar */
537 radar = (Radar *)RSL_new_radar(MAX_RADAR_VOLUMES);
540 perror("RadarBuild(): Error creating radar structure.\n");
543 /* Initialize the Radar_header values. */
544 Radar_headerFill(radar, gvl1);
546 /* Build each of the 'nparm' volumes of the radar structure. */
547 for (pindex=0; pindex<vs->tk.nparm; pindex++)
549 /* Identify parameter type, so we know which RSL volume to load the data
551 vindex = parmIdentify(gvl1->sensor.parm[pindex]->parmDesc.parmName);
555 "RadarBuild(): Unexpected parameter type: %s found in HDF file.\n",
556 gvl1->sensor.parm[pindex]->parmDesc.parmName);
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. */
562 if (radar_verbose_flag)
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");
571 } /* end else if (rsl_qfield[vindex] == 0) */
572 else if (vindex == CZ_INDEX) /* Handle CZ and DZ volumes. */
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])
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);
583 } /* end if (vindex == CZ_INDEX) */
584 else if (vindex == CD_INDEX) /* Handle CD and ZD volumes. */
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])
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);
595 } /* end if (vindex == CD_INDEX) */
596 else /* Handle all 1B-51 fields. (DZ, ZD, VR) */
598 radar->v[vindex] = VolumeBuild(gvl1, vs, 0.0, vindex, pindex);
600 } /* end for (pindex=0; ...) */
605 /*************************************************************/
609 /*************************************************************/
610 int commentsRead(VosSize *vs, float *zCal, char *comments, int productID)
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.
616 Returns: OK if success.
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];
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.*/
631 strcpy(parseString, "");
632 strcat(parseString, "%[^*] %*[*\n]");
634 if (productID == TK_L1C_GV)
636 strcat(parseString, "%[^\n]");
639 /* Read all records from the comments field into the record buffers. */
640 if (sscanf(comments, parseString, record[0], record[1]) != nrecords)
644 if (sscanf(record[0], "nSweep=%d", &vs->tk.nsweep) != 1) goto quit;
646 strcpy(parseString, "nRay=%d\n");
647 for (pindex=0; pindex<vs->tk.nparm; pindex++)
648 strcat(parseString, "nCell_parm[%*d]=%d\n");
650 spointer = record[0];
651 for (tk_sindex=0; tk_sindex<vs->tk.nsweep; tk_sindex++)
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;
662 /* If 1C-51 file, read the QC parameters into the qcParm array. */
663 if (productID == TK_L1C_GV)
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)
670 /* Print out the QC parameters we've just read in. */
672 if (radar_verbose_flag)
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]);
683 if (qcParm[ZCAL] <= NOVAL_FLOAT) *zCal = 0.0;
684 else *zCal = qcParm[ZCAL];
685 } /* end if (productID == TK_L1C_GV) */
690 if (radar_verbose_flag)
691 fprintf(stderr, "commentsRead(): Failure reading comments field\n");
695 /*************************************************************/
699 /*************************************************************/
700 static L1B_1C_GV *GVL1Build(IO_HANDLE *granuleHandle, int vosNum,
703 /* Build a toolkit 'L1B_1C_GV' structure sized for the VOS we
704 will later read in from the HDF file.
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)
719 fprintf(stderr, "GVL1Build(): TKgetNparm() failed.\n");
722 vs->tk.nsweep = TKgetNsweep(granuleHandle, vosNum);
723 if (vs->tk.nsweep <= 0)
725 fprintf(stderr, "GVL1Build(): TKgetNsweep() failed.\n");
728 vs->rsl.maxNray = TKgetNray(granuleHandle, vosNum);
729 if (vs->rsl.maxNray <= 0)
731 fprintf(stderr, "GVL1Build(): TKgetNray() failed.\n");
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++)
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)
745 gvl1->sensor.parm[pindex]->parmData2byte =
746 (int16 ***)TKnewParmData2byte(vs->tk.nsweep, vs->rsl.maxNray, ncell);
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);
755 gvl1->sensor.parm[pindex]->parmData2byte =
756 (int16 ***)TKnewParmData2byte(vs->tk.nsweep, vs->rsl.maxNray, ncell);
758 } /* end for (pindex=0; ... */
763 /*************************************************************/
767 /*************************************************************/
768 int metaDataRead(Radar *radar, IO_HANDLE *granuleHandle)
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);
779 /*************************************************************/
783 /*************************************************************/
784 static int hdfFileOpen(char *infile, IO_HANDLE *granuleHandle,
785 char *hdfFileName, int *vosNum)
787 /* Opens, if necessary, an HDF file. Checks that a VOS, not previously
788 retrieved, exists in the HDF file.
795 int productID, nvos, status;
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. */
802 if (strcmp(hdfFileName, infile) != 0)
804 if (TKclose(granuleHandle) != TK_SUCCESS)
806 fprintf(stderr, "hdfFileOpen(): *** TKclose() error\n");
811 } /* end if (*vosNum != 0) */
813 /* If first VOS of HDF file, we need first to open the file. */
816 strncpy(hdfFileName, infile, TK_MAX_FILENAME-1);
817 /* Get the desired product out of the HDF filename. */
818 product = strrchr(hdfFileName, '/');
820 product = hdfFileName;
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;
829 fprintf(stderr, "hdfFileOpen(): Invalid HDF filename.\n");
832 status = TKopen(hdfFileName, productID, TK_READ_ONLY, granuleHandle);
833 if (status != TK_SUCCESS)
835 fprintf(stderr, "hdfFileOpen(): *** TKopen() error\n");
838 } /* end if (*vosNum == 0) */
840 /* Check if the requested VOS exists in the HDF granule. */
841 nvos = (int)TKgetNvos(granuleHandle);
844 if (radar_verbose_flag)
845 fprintf(stderr, "\nEmpty granule.\n");
848 else if (*vosNum+1 > nvos)
850 if (radar_verbose_flag)
851 fprintf(stderr, "\nAll VOSs read from HDF file: %s\n", hdfFileName);
856 fprintf(stderr, "hdfFileOpen():*** TKgetNvos() error\n");
863 /*************************************************************/
865 /* RSL_hdf_to_radar */
867 /*************************************************************/
868 Radar *RSL_hdf_to_radar(char *infile)
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.
877 - pointer to a filled radar structure, if success.
878 - NULL pointer, if failure.
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) */
890 /* Open, if necessary, an HDF file. */
891 status = hdfFileOpen(infile, &granuleHandle, hdfFileName, &vosNum);
892 if (status < 0) goto quit;
894 /* Initialize the 'VosSize' structure. */
895 memset(&vs, '\0', sizeof(VosSize));
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;
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)
908 fprintf(stderr, "RSL_hdf_to_radar(): *** TKreadL1GV() error\n");
912 if (radar_verbose_flag)
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);
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);
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.
932 if (commentsRead(&vs, &zCal, gvl1->comments, granuleHandle.productID) < 0)
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
939 status = metaDataRead(radar, &granuleHandle);
941 /* Free memory allocated to the toolkit GVL1 structure. */
943 if (radar == NULL) goto quit;
950 TKclose(&granuleHandle);
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.
962 Radar *RSL_hdf_to_radar(char *infile)