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 /******************************************************************
29 Writes one VOS from a RSL radar structure into one 1B-51/1C-51 HDF
32 A 1B-51/1C-51 HDF file contains multiple VOS's recorded by a radar
33 site during a 1-hour time period, and the HDF file is named using
34 the date/hour of the constituent VOS's.
36 Functions defined herein perform tasks preparatory to building a
37 TSDIS toolkit 'L1B_1C_GV' structure, the contents of which are
38 written into the HDF file by the TSDIS toolkit. Construction of the
39 toolkit 'L1B_1C_GV' structure is done via the subroutines defined
40 in RSL file 'radar_to_hdf_2.c'.
42 -----------------------------------------------------------------
43 Libraries required for execution of this code :
44 -ltsdistk : TSDIS toolkit
45 -lmfhdf -ldf -ljpeg -lz : HDF
49 -----------------------------------------------------------------
50 *******************************************************************/
58 #include <sys/types.h>
62 /* TSDIS toolkit function and structure definitions. */
65 /* RSL function and structure definitions. */
67 /* Parameter definitions for 1B-51 and 1C-51 HDF
68 file handling applications using the TSDIS toolkit. */
69 #include "toolkit_1BC-51_appl.h"
71 /*************************************************************/
73 /* Function Prototypes */
75 /*************************************************************/
76 void RSL_set_tkMetaDataString(char *string, int param);
77 void RSL_set_hdf_qc_parameters(float *inParm);
78 void radarVolumesSave(Volume *v[MAX_RADAR_VOLUMES], Radar *radar);
79 void radarVolumesRestore(Radar *radar, Volume *v[MAX_RADAR_VOLUMES]);
80 int nextVolume(Radar *radar, int last_volume);
81 Volume *maskBuild(Volume *cv, Volume *ucv);
82 void metaDataWrite(IO_HANDLE *granuleHandle, Radar *radar,
83 char *hdfFileName, int fileAccessMode);
84 int tkVosDimensions(VosSize *vs, Radar *radar);
85 Ray *first_ray_in_sweep(Sweep *sweep);
86 int rslVosDimensions(VosSize *vs, Radar *radar, float maxRange);
87 int L1GVtemplateInit(VosSize *vs, Radar *radar, char *hdfFileName,
89 static int hdfFileOpen(IO_HANDLE *granuleHandle, VosSize *vs,
90 char *hdfFileName, Radar *radar);
91 int radarPrep1B51(Radar *radar);
92 int radarPrep1C51(Radar *radar);
93 int nullGranuleCreate(char *hdfFileName, IO_HANDLE *granuleHandle,
95 int RSL_radar_to_hdf(Radar *radar, char *hdfFileName);
97 extern L1B_1C_GV *gvl1Build(Radar *radar, float *qcParm, VosSize *vs,
99 extern int radar_verbose_flag;
101 /* The 1st non-NULL ray in each volume is widely used. Hence global. */
102 Ray *first_ray_in_volume[MAX_RADAR_VOLUMES];
106 /*************************************************************/
108 /* RSL_set_tkMetaDataString */
110 /*************************************************************/
113 char GenInputDate[128];
114 char AlgorithmVersion[128];
116 char SoftwareVersion[128];
119 void RSL_set_tkMetaDataString(char *string, int param)
121 /* Allows the application 'level_1' to store metadata strings into
122 the tkMetaDataString buffers for later insertion by RSL function
123 'metaDataWrite()' into the HDF file.
124 Call **before** RSL_radar_to_hdf().
127 #define CP_TKMETA(s, str) \
128 memset(tkMetaDataString.s, '\0', sizeof(tkMetaDataString.s));\
129 strncpy(tkMetaDataString.s, str, sizeof(tkMetaDataString.s)-1)
132 if (param == TK_GEN_DATE_INPUT_FILES) {
133 CP_TKMETA(GenInputDate, string);
135 To decide whether or not to write a VOS into a 1C-51 HDF file
136 (See 'hdf1C51Create()' in file 'level_1.c'), I need to know
137 something about the times of the VOSs already in the file;
138 ie, if no satellite overpass, the file is to contain just one
139 VOS from each half-hour. The only efficient way to do this is
140 by encoding time info into a metadata field.
142 The 'GenInputDate' metaData string contains info indicating
143 which hourly time slots are filled by the VOSs contained in the
144 HDF file. (See function 'metaDataWrite()' in this file for details.)
146 Called by 'hdfFilePeek()' in file 'level_1.c'.
149 else if (param == TK_ALGORITHM_VERSION) {
150 CP_TKMETA(AlgorithmVersion, string);
152 } else if (param == TK_SOFTWARE_VERSION) {
153 CP_TKMETA(SoftwareVersion, string);
155 if (radar_verbose_flag)
156 fprintf(stderr, "TK_SOFTWARE_VERSION = <%s>\n", string);
158 } else if (param == TK_PRODUCT_VERSION) {
159 sscanf(string, "%d", &tkMetaDataString.ProductVersion);
163 fprintf(stderr, "RSL_set_tkMetaDataString: Unknown param==%d!!\n", param);
167 /*************************************************************/
169 /* RSL_set_hdf_qc_parameters */
171 /*************************************************************/
172 static float qcParm[NUMBER_QC_PARAMS] =
174 NOVAL_FLOAT, NOVAL_FLOAT, NOVAL_FLOAT, NOVAL_FLOAT, NOVAL_FLOAT,
175 NOVAL_FLOAT, NOVAL_FLOAT, NOVAL_FLOAT, NOVAL_FLOAT, NOVAL_FLOAT
178 void RSL_set_hdf_qc_parameters(float *inParm)
180 /* Stores 1C-51 QC parameters for later insertion into 1C-51 HDF
181 file. Call **before** RSL_radar_to_hdf().
185 for (j=0; j<NUMBER_QC_PARAMS; j++)
186 qcParm[j] = inParm[j];
189 /*************************************************************/
191 /* radarVolumesSave */
193 /*************************************************************/
194 void radarVolumesSave(Volume *v[MAX_RADAR_VOLUMES], Radar *radar)
196 /* Save array of radar volume pointers. */
199 for (j=0; j<radar->h.nvolumes; j++)
203 /*************************************************************/
205 /* radarVolumesRestore */
207 /*************************************************************/
208 void radarVolumesRestore(Radar *radar, Volume *v[MAX_RADAR_VOLUMES])
210 /* Restore array of radar volume pointers. */
213 for (j=0; j<radar->h.nvolumes; j++)
217 /*************************************************************/
221 /*************************************************************/
222 int nextVolume(Radar *radar, int last_volume)
223 /* Find the index of the next volume in the radar structure, given
224 the index of the last volume found.
226 index of next volume, if success.
232 for (j=last_volume+1; j<radar->h.nvolumes; j++)
233 if (radar->v[j] != NULL)
234 return(j); /* Found volume. Return the index. */
236 return(-1); /* No volume found. Return bogus index. */
239 /*************************************************************/
243 /*************************************************************/
244 Volume *maskBuild(Volume *cv, Volume *ucv)
246 /* Create a mask volume. A mask value, in conjunction with a corresponding
247 corrected reflectivity value, enables the future recovery of a
248 uncorrected value from a 1C-51 HDF file, since the original, uncorrected
249 reflectivity values are not stored in 1C-51 HDF files.
251 This function used for both reflectivity 'DZ' and differential reflectivity
254 cv: corrected volume pointer. (CZ or CD)
255 ucv: uncorrected volume pointer. (DZ or ZD)
256 mv: mask volume pointer. (MZ or MD)
259 int sindex, rindex, bindex;
262 mv = RSL_copy_volume(cv);
263 mv->h.f = MZ_F; /* MZ_F identical to MD_F. Can use either. */
264 mv->h.invf = MZ_INVF; /* MZ_INVF identical to MD_INVF. Can use either. */
265 for (sindex=0; sindex<cv->h.nsweeps; sindex++)
267 if (cv->sweep[sindex] == NULL) continue;
268 mv->sweep[sindex]->h.f = MZ_F;
269 mv->sweep[sindex]->h.invf = MZ_INVF;
270 for (rindex=0; rindex<cv->sweep[sindex]->h.nrays; rindex++)
272 if (cv->sweep[sindex]->ray[rindex] == NULL) continue;
273 mv->sweep[sindex]->ray[rindex]->h.f = MZ_F;
274 mv->sweep[sindex]->ray[rindex]->h.invf = MZ_INVF;
275 for (bindex=0; bindex<cv->sweep[sindex]->ray[rindex]->h.nbins; bindex++)
276 /* Has the Z_value been corrected by the QC algorithm? */
278 if (cv->sweep[sindex]->ray[rindex]->range[bindex] == cv->h.invf(BADVAL))
279 mv->sweep[sindex]->ray[rindex]->range[bindex] = 1;
280 else /* Uncorrected Z_value. */
281 mv->sweep[sindex]->ray[rindex]->range[bindex] = 0;
282 } /* end for (rindex=0; ... */
283 } /* end for (sindex=0; ... */
288 /*************************************************************/
292 /*************************************************************/
293 void metaDataWrite(IO_HANDLE *granuleHandle, Radar *radar,
294 char *hdfFileName, int fileAccessMode)
296 /* Write out some metadata values into toolkit structures. */
298 * YOU REALLY DON'T KNOW, A PRIORI, WHAT DATATYPE THE METADATA ITEM IS.
299 * TRIAL and ERROR is how we figured some of them out.
301 * THEY SHOULD ALL BE STRINGS!!!!!
309 vindex = nextVolume(radar, -1); /* Find 1st non-NULL rsl volume. */
311 /* ----------------- Core MetaData ---------------*/
312 /* Begin/end date of the HDF granule. */
313 tkdate.tkyear = (short) radar->h.year;
314 tkdate.tkmonth = (short) radar->h.month;
315 tkdate.tkday = (short) radar->h.day;
316 TKwriteMetadataInt(granuleHandle, TK_BEGIN_DATE, &tkdate);
317 TKwriteMetadataInt(granuleHandle, TK_END_DATE, &tkdate);
318 /* Begin time of the HDF granule. */
319 tktime.tkhour = (int8) radar->h.hour;
320 tktime.tkminute = (int8) 0;
321 tktime.tksecond = (int8) 0;
322 TKwriteMetadataInt(granuleHandle, TK_BEGIN_TIME, &tktime);
323 /* End time of the HDF granule. */
324 tktime.tkhour = (int8) radar->h.hour;
325 tktime.tkminute = (int8) 59;
326 tktime.tksecond = (int8) 59;
327 TKwriteMetadataInt(granuleHandle, TK_END_TIME, &tktime);
329 /* Longitude & Bounding Coordinates */
330 floatVal = radar->h.lond + radar->h.lonm/60.0 + radar->h.lons/3600.0;
331 TKwriteMetadataFloat(granuleHandle, TK_CENTER_POINT_LON, &floatVal);
332 floatVal = floatVal - 2.0;
333 TKwriteMetadataFloat(granuleHandle, TK_WEST_BOUND_COORD, &floatVal);
334 floatVal = floatVal + 4.0;
335 TKwriteMetadataFloat(granuleHandle, TK_EAST_BOUND_COORD, &floatVal);
336 /* Latitude & Bounding Coordinates */
337 floatVal = radar->h.latd + radar->h.latm/60.0 + radar->h.lats/3600.0;
338 TKwriteMetadataFloat(granuleHandle, TK_CENTER_POINT_LAT, &floatVal);
339 floatVal = floatVal - 2.0;
340 TKwriteMetadataFloat(granuleHandle, TK_SOUTH_BOUND_COORD, &floatVal);
341 floatVal = floatVal + 4.0;
342 TKwriteMetadataFloat(granuleHandle, TK_NORTH_BOUND_COORD, &floatVal);
344 TKwriteMetadataChar(granuleHandle, TK_CONTACT, "Danny Rosenfeld");
348 /* ----------------- PS MetaData ---------------*/
349 TKwriteMetadataChar(granuleHandle, TK_ALGORITHM_VERSION,
350 tkMetaDataString.AlgorithmVersion);
351 TKwriteMetadataInt(granuleHandle, TK_PRODUCT_VERSION,
352 &tkMetaDataString.ProductVersion);
353 TKwriteMetadataChar(granuleHandle, TK_SOFTWARE_VERSION,
354 tkMetaDataString.SoftwareVersion);
355 TKwriteMetadataChar(granuleHandle, TK_MAX_VALID_CHANNEL, "70 dBz");
356 TKwriteMetadataChar(granuleHandle, TK_MIN_VALID_CHANNEL, "-20 dBz");
357 if (radar->h.nvolumes > 0)
359 floatVal = first_ray_in_volume[vindex]->h.wavelength;
361 else /* No radar volumes. Creating empty granule. */
363 TKwriteMetadataChar(granuleHandle, TK_ANOMALY_FLAG, "EMPTY: REASON UNKNOWN");
366 TKwriteMetadataFloat(granuleHandle, TK_RADAR_WAVELENGTH, &floatVal);
368 floatVal = -20.0; /* MIN_REFL_THRESHOLD */
369 TKwriteMetadataFloat(granuleHandle, TK_MIN_REF_THRESHOLD, &floatVal);
370 TKwriteMetadataChar(granuleHandle, TK_RADAR_NAME, radar->h.radar_name);
371 TKwriteMetadataChar(granuleHandle, TK_RADAR_CITY, radar->h.city);
372 TKwriteMetadataChar(granuleHandle, TK_RADAR_STATE, radar->h.state);
373 TKwriteMetadataChar(granuleHandle, TK_RADAR_COUNTRY, radar->h.country);
375 intVal = (int)TKgetNvos(granuleHandle);
376 TKwriteMetadataInt(granuleHandle, TK_NUM_VOS, &intVal);
377 TKwriteMetadataFloat(granuleHandle, TK_GV_DZCAL, &qcParm[ZCAL]);
378 floatVal = X; /* Mask scale factor */
379 TKwriteMetadataFloat(granuleHandle, TK_GV_L1C_SCALE, &floatVal);
380 floatVal = 0.0; /* Correction for gaseous two-way attenuation */
381 TKwriteMetadataFloat(granuleHandle, TK_GV_ALPHA, &floatVal);
383 TKwriteMetadataChar(granuleHandle, TK_INPUT_FILES, "8mm tape files");
384 TKwriteMetadataChar(granuleHandle, TK_DATA_CENTER_SRC, radar->h.radar_name);
386 /* I really, really hate this 1C-51 kludge, but...
387 To decide (within function 'hdf1C51Create') whether or not to write
388 a VOS into a 1C-51 file, I need to know something about the times
389 of the VOSs already in the file; ie, if no satellite overpass, the
390 file is to contain just one VOS from each half-hour.
391 The only efficient way to do this is by encoding info into a metadata
394 So... I use the 'GEN_DATE_INPUT_FILES' metaData field. The intial
395 value of this field is the string 'unKNOWN'. When a VOS from the
396 first half_of_the_hour is written into the file, change the first
397 char 'u' to upper_case. When a VOS from the second half_of_the_hour
398 is written into the file, change the second char 'n' to upper_case.
399 The final string, after both time slots are filled, is entirely
400 upper_case: 'UNKOWN'.
402 if (granuleHandle->productID == TK_L1C_GV) /* 1C-51 file? */
404 /* If new file, then intitialize 'GEN_DATE_INPUT_FILES' string.
405 If file exists, then get the existing string from the file. */
406 if (fileAccessMode == TK_NEW_FILE) strcpy(buf, "unKNOWN");
407 else strcpy(buf, tkMetaDataString.GenInputDate);
408 if (radar->h.nvolumes > 0) /* Don't mess with an empty granule. */
410 if (radar->h.minute < 30) /* 1st half_of_the_hour? */
411 buf[0] = (char) toupper((int)buf[0]); /* Convert char to upper_case */
412 else /* 2nd half_of_the_hour */
413 buf[1] = (char) toupper((int)buf[1]); /* Convert char to upper_case */
416 else /* 1B-51 file */
418 strcpy(buf, "UNKNOWN");
420 TKwriteMetadataChar(granuleHandle, TK_GEN_DATE_INPUT_FILES, buf);
423 /*************************************************************/
425 /* newPhysicalSweep */
427 /*************************************************************/
428 int newPhysicalSweep(Sweep *phys_sweep, Sweep *sweep)
430 /* Checks if the rsl 'sweep' belongs to a different physical sweep
432 Returns: 1, if sweep is a new physical sweep
433 0, if not new physical sweep
438 /* If elevations don't match, new sweep. */
439 if (phys_sweep->h.elev != sweep->h.elev) return(1);
441 /* Check the first 50 rays of both sweeps.
442 * If azimuths don't match, new sweep.
445 if (sweep->h.nrays < nrays) nrays = sweep->h.nrays;
446 if (phys_sweep->h.nrays < nrays) nrays = phys_sweep->h.nrays;
447 for (iray=0; iray<nrays; iray++)
449 if (sweep->ray[iray] == NULL) continue;
450 if (phys_sweep->ray[iray] == NULL) continue;
451 if (phys_sweep->ray[iray]->h.azimuth != sweep->ray[iray]->h.azimuth)
452 return(1); /* New physical sweep */
455 /* No new physical sweep. */
459 /*************************************************************/
461 /* tkVosDimensions */
463 /*************************************************************/
464 int tkVosDimensions(VosSize *vs, Radar *radar)
466 /* Logically reconfigures the VOS in the rsl radar structure into
467 a sequence of physical sweeps as required for the toolkit
470 Records required toolkit gvl1 dimensions in the 'VosSize->tk'
471 structure. When we later load the toolkit gvl1 structure with
472 actual data values, vs->tk.ncell[tk_sindex][pindex] will be zero
473 for those data types not collected by the radar during physical
476 The algorithm used herein to locate physical sweeps is as simple
477 as possible. It works for all VOSs processed to date. However,
478 if more complicated radar scanning regimes turn up, this function
479 will surely require an overhaul and additional complexity.
481 Algorithm: Within the rsl radar structure, walks horizontally
482 across the radar volumes at one sweep level, checking for elev
483 and azim values different from the last physical sweep located.
484 Then increments the sweep level and repeats at the next rsl sweep
487 Returns: OK, if success.
490 int pindex, sindex, tk_sindex, status;
493 /* Note: 'vs->tk' structure contains all zeroes upon entry into this
495 sindex = 0; /* rsl sweep index */
496 tk_sindex = -1; /* toolkit sweep index */
497 while (sindex < vs->rsl.maxNsweep) /* for each rsl sweep... */
499 for (pindex=0; pindex<vs->tk.nparm; pindex++) /* for each rsl volume */
501 if (sindex >= vs->rsl.nsweep[pindex]) continue;
502 sweep = vs->rsl.v[pindex]->sweep[sindex];
503 /* No null sweep ptrs allowed in rsl ptr array. */
504 if (sweep == NULL) return(QUIT);
505 /* Check for a new physical sweep. */
509 vs->tk.nray[tk_sindex] = vs->rsl.nray[pindex][sindex];
510 vs->rsl.sweep[tk_sindex] = sweep;
512 else if ((status=newPhysicalSweep(vs->rsl.sweep[tk_sindex], sweep)))
514 if (status < 0) return(status);
516 if (tk_sindex >= MAX_SWEEP)
518 if (radar_verbose_flag)
519 fprintf(stderr, "tkVosDimensions(): Too many toolkit sweeps.\n");
522 vs->tk.nray[tk_sindex] = vs->rsl.nray[pindex][sindex];
523 vs->rsl.sweep[tk_sindex] = sweep;
524 } /* end if (newPhysicalSweep */
526 vs->tk.ncell[tk_sindex][pindex] = vs->rsl.ncell[pindex][sindex];
527 } /* end for (pindex=0;... */
530 } /* end while (sindex < vs->rsl.maxNsweep) */
531 vs->tk.nsweep = tk_sindex + 1;
535 /*************************************************************/
537 /* first_ray_in_sweep */
539 /*************************************************************/
540 Ray *first_ray_in_sweep(Sweep *sweep)
543 * Return the first non-NULL ray in the sweep.
544 * Returns NULL, if error.
548 if (sweep == NULL) return(NULL);
549 for (iray=0; iray<sweep->h.nrays; iray++)
551 if (sweep->ray[iray] != NULL)
552 return(sweep->ray[iray]);
558 /*************************************************************/
560 /* rslVosDimensions */
562 /*************************************************************/
563 int rslVosDimensions(VosSize *vs, Radar *radar, float maxRange)
565 /* Scopes out all dimensions of the VOS contained in the rsl
566 radar structure. Records dimensions in the 'VosSize->rsl'
569 The ncell values may be less than the actual number of
570 bins in the radar structure, since the 1B-51/1C-51 standards
571 call for range truncation. Truncates range as necessary.
572 For 1B-51: Max range is the lesser of:
573 1: max_range from radar structure, and
576 For 1C-51: Max range is the lesser of:
577 1: max_range from radar structure, and
584 /* Must differentiate between rsl and vosSize array indices,
585 since the rsl arrays may contain NULL elements. */
586 int ivolume, isweep, iray; /* Indices for rsl arrays. */
587 int Vindex, Sindex; /* Indices for non-NULL rsl array elements. */
589 float maxRangeActual = 0.0;
591 Ray *longest_ray_in_sweep;
593 /* Initialize the 'VosSize' structure. */
594 memset(vs, '\0', sizeof(VosSize));
595 vs->rsl.sweep[0] = NULL;
598 if (radar_verbose_flag)
599 fprintf(stderr, "RSL VOS Dimensions...\n");
602 for (ivolume=0; ivolume<radar->h.nvolumes; ivolume++)
604 if (radar->v[ivolume] == NULL) continue;
606 vs->rsl.v[Vindex] = radar->v[ivolume];
608 for (isweep=0; isweep<radar->v[ivolume]->h.nsweeps; isweep++)
610 if (radar->v[ivolume]->sweep[isweep] == NULL) continue;
612 longest_ray_in_sweep = first_ray_in_sweep(radar->v[ivolume]->sweep[isweep]);
613 if (longest_ray_in_sweep == NULL)
615 fprintf(stderr, "rslVosDimensions(): no rays in sweep:%d ???\n\n",
620 if (Sindex == 0) /* First sweep in volume? */
621 first_ray_in_volume[ivolume] = longest_ray_in_sweep;
623 * Go thru the entire sweep to find the number of
624 * rays, and the longest ray.
626 for (iray=0; iray<radar->v[ivolume]->sweep[isweep]->h.nrays; iray++)
628 if (radar->v[ivolume]->sweep[isweep]->ray[iray] == NULL) continue;
629 ray = radar->v[ivolume]->sweep[isweep]->ray[iray];
631 if (ray->h.nbins > longest_ray_in_sweep->h.nbins)
632 longest_ray_in_sweep = ray;
633 } /* end for (iray=0;... */
635 if (Sindex == 0) /* 1st sweep of this volume? */
637 /* Find max range (km) of data in this sweep. */
638 maxRangeActual = (longest_ray_in_sweep->h.range_bin1 +
639 longest_ray_in_sweep->h.nbins *
640 longest_ray_in_sweep->h.gate_size) / 1000.0; /*km*/
641 if (maxRangeActual > maxRange+0.5) /* Truncate range at maxRange km */
642 vs->rsl.ncell[Vindex][Sindex] = (int)
643 ((maxRange*1000.0 - longest_ray_in_sweep->h.range_bin1) /
644 longest_ray_in_sweep->h.gate_size);
646 vs->rsl.ncell[Vindex][Sindex] = (int)
647 ((maxRangeActual*1000.0 - longest_ray_in_sweep->h.range_bin1) /
648 longest_ray_in_sweep->h.gate_size);
649 if (vs->rsl.ncell[Vindex][Sindex] > MAX_CELL)
651 fprintf(stderr, "rslVosDimensions(): ncell[parm%d]=%d > MAX_CELL=%d ",
652 Vindex, vs->rsl.ncell[Vindex][Sindex], MAX_CELL);
653 fprintf(stderr, " gate_size:%d\n", longest_ray_in_sweep->h.gate_size);
656 } /* end if (Sindex == 0) */
657 else /* Not the 1st sweep of volume. */
659 if (longest_ray_in_sweep->h.nbins > vs->rsl.ncell[Vindex][0])
660 vs->rsl.ncell[Vindex][Sindex] = vs->rsl.ncell[Vindex][0];
662 vs->rsl.ncell[Vindex][Sindex] = longest_ray_in_sweep->h.nbins;
665 vs->rsl.nray[Vindex][Sindex] = nrays_in_sweep;
666 if (vs->rsl.nray[Vindex][Sindex] > vs->rsl.maxNray)
668 vs->rsl.maxNray = vs->rsl.nray[Vindex][Sindex];
669 if (vs->rsl.maxNray > MAX_RAY)
671 fprintf(stderr, "rslVosDimensions(): v[%d]->sweep[%d].nray=%d > MAX_RAY=%d\n",
672 ivolume, isweep, vs->rsl.nray[Vindex][Sindex], MAX_RAY);
676 } /* end for (isweep=0;... */
677 vs->rsl.nsweep[Vindex] = Sindex + 1;
678 if (vs->rsl.nsweep[Vindex] > vs->rsl.maxNsweep)
680 vs->rsl.maxNsweep = vs->rsl.nsweep[Vindex];
681 if (vs->rsl.maxNsweep > MAX_SWEEP)
683 fprintf(stderr, "rslVosDimensions(): v[%d].nsweep=%d > MAX_SWEEP=%d\n",
684 ivolume, vs->rsl.nsweep[Vindex], MAX_SWEEP);
689 if (radar_verbose_flag)
692 " vIndex:%2d nsweeps:%d cellSize(m):%4d ncells:%d maxRng(km):%.1f\n",
693 ivolume, vs->rsl.nsweep[Vindex], longest_ray_in_sweep->h.gate_size,
694 longest_ray_in_sweep->h.nbins, maxRangeActual);
697 } /* for (ivolume=0;... */
698 vs->tk.nparm = Vindex + 1;
699 if ((vs->tk.nparm == 0) || (vs->tk.nparm > MAX_PARM))
701 fprintf(stderr, "rslVosDimensions(): Invalid nparm=%d\n", vs->tk.nparm);
708 /*************************************************************/
710 /* L1GVtemplateInit */
712 /*************************************************************/
713 int L1GVtemplateInit(VosSize *vs, Radar *radar, char *hdfFileName,
716 /* Create a toolkit 'Level_1B_1C_GV' template_node for the VOS
717 contained in the radar structure. To do this:
718 1: Based on the data contained in the rsl structure, initialize
719 the following toolkit arrays:
720 'TKnparm' : no. of volumes in radar structure (DZ, VR, ZD, etc).
721 'TKnsweep': max no. of sweeps/volume, over all volumes.
722 'TKnray' : max no. of rays per sweep, over all volumes.
723 'TKncell': max no. of cells (bins) per ray, over all volumes, all sweeps.
724 These values are parameters for the 'TKsetL1GVtemplate()'
725 toolkit function call.
726 2: Call TKsetL1GVtemplate().
733 int32 TKnparm[MAX_VOS], TKnsweep[MAX_VOS], TKnray[MAX_VOS];
734 int32 TKncell[MAX_VOS][MAX_PARM];
736 /* Scope out all dimensions of the VOS contained in the rsl structure. */
737 status = rslVosDimensions(vs, radar, maxRange);
738 if (status < 0) return(status);
739 /* Set the toolkit VOS dimensions. These are different from the rsl
740 dimensions. The toolkit L1GV structure is organized as a
741 sequence of physical sweeps, while RSL is organized as a
742 sequence of logical "volumes".
744 status = tkVosDimensions(vs, radar);
745 if (status < 0) return(status);
747 /* Fill required toolkit array values for TKsetL1GVtemplate() call. */
748 TKnparm[0] = vs->tk.nparm;
749 TKnsweep[0] = vs->tk.nsweep;
750 TKnray[0] = vs->rsl.maxNray;
752 if (radar_verbose_flag)
753 fprintf(stderr, "Toolkit VOS Dimensions...\n");
755 for (pindex=0; pindex<vs->tk.nparm; pindex++)
757 TKncell[0][pindex] = vs->rsl.ncell[pindex][0];
759 if (radar_verbose_flag)
760 fprintf(stderr, " pIndex:%d nsweep:%d nray:%d ncell:%d\n", pindex,
761 (int)TKnsweep[0], (int)TKnray[0], (int)TKncell[0][pindex]);
765 /* Create a toolkit template_node for this VOS. */
766 status = TKsetL1GVtemplate(1, TKnparm, TKncell, TKnray, TKnsweep, hdfFileName);
767 if (status != TK_SUCCESS)
769 fprintf(stderr, "L1GVtemplateInit(): *****TKsetL1GVtemplate() error\n");
773 return(OK); /* Successful template_node creation. */
776 /*************************************************************/
780 /*************************************************************/
781 static int hdfFileOpen(IO_HANDLE *granuleHandle, VosSize *vs,
782 char *hdfFileName, Radar *radar)
784 /* Create a toolkit template node for the new VOS, and then open
785 a HDF file in which to write the VOS.
792 int productID, status;
796 /* Based on the product desired, set the maximum range of radar data
797 values for the HDF file.
798 For 1B-51: Max range = 230.0 km.
799 For 1C-51: Max range = 200.0 km
801 productID = (int)granuleHandle->productID;
802 if (productID == TK_L1B_GV)
803 maxRange = MAX_RANGE_1B51;
804 else if (productID == TK_L1C_GV)
805 maxRange = MAX_RANGE_1C51;
808 fprintf(stderr, "hdfFileOpen(): Invalid product type\n");
812 /* Create a L1BGV template node for the new VOS. */
814 if (radar_verbose_flag)
815 fprintf(stderr, "\n****** Creating toolkit template_node for VOS ...\n");
817 status = L1GVtemplateInit(vs, radar, hdfFileName, maxRange);
818 if (status < 0) return(status);
820 /* Determine if the HDF file to which this VOS belongs already exists.*/
821 if (stat(hdfFileName, &buf) == 0)
823 /* The HDF file already exists. We will append this VOS to it. */
824 fileAccessMode = TK_APPEND;
825 if (radar_verbose_flag)
826 fprintf(stderr, "\n****** Opening HDF file: %s to append VOS ...\n",
829 else /* The HDF file does not exist. We must create it. */
831 fileAccessMode = TK_NEW_FILE;
832 if (radar_verbose_flag)
833 fprintf(stderr, "\n****** Opening new HDF file: %s to write VOS ...\n",
837 /* Finally, open the HDF file. */
838 status = TKopen(hdfFileName, productID, fileAccessMode, granuleHandle);
839 if (status != TK_SUCCESS)
841 if (radar_verbose_flag)
842 fprintf(stderr, "level_1(): ***** TKopen() error\n");
846 /* Find the slot number in the HDF granule for this VOS. */
847 vs->vos_num = TKgetNvos(granuleHandle) - 1;
849 /* Write metadata fields into HDF file. */
850 metaDataWrite(granuleHandle, radar, hdfFileName, fileAccessMode);
855 /*************************************************************/
859 /*************************************************************/
860 int radarPrep1B51(Radar *radar)
863 Prepare the RSL radar structure for 1B-51 processing:
864 Null the pointers to RSL volumes not required for 1B-51.
870 /* For 1B-51, we need only the DZ, VR, and ZD data volumes. NULL the
871 pointers to all other volumes. */
872 for (j=0; j<radar->h.nvolumes; j++)
874 if (radar->v[j] == NULL)
876 if ((j == DZ_INDEX) || (j == VR_INDEX) || (j == ZD_INDEX))
879 } /* end for (j=0; ... */
884 /*************************************************************/
888 /*************************************************************/
889 int radarPrep1C51(Radar *radar)
892 Prepare the RSL radar structure for 1C-51 processing.
893 1. Check that the necessary RSL volumes exist.
894 2. Create mask volumes.
895 3. Null the pointers to RSL volumes not required for 1C-51.
903 if (radar->v[DZ_INDEX] != NULL) /* Is there a DZ volume? */
905 if (radar->v[CZ_INDEX] == NULL) /* DZ exists, hence CZ should exist. */
907 fprintf(stderr, "RSL_radar_to_hdf(): CZ volume expected but not found.\n");
910 else /* Both DZ and CZ volumes exist. */
912 /* Construct mask volume MZ */
913 radar->v[MZ_INDEX] = maskBuild(radar->v[CZ_INDEX], radar->v[DZ_INDEX]);
914 /* We're now finished with the RSL CZ data, so let the CZ volume point
915 to the uncorrected volume DZ. We will later use it,
916 in combination with the mask volume MZ, to create the HDF CZ volume.
917 (The HDF CZ values differ from the RSL CZ values.)
919 radar->v[CZ_INDEX] = radar->v[DZ_INDEX];
921 } /* end if (radar->v[DZ_INDEX] != NULL) */
923 /* For 1C-51, we need only the following volumes:
924 CZ: QC'ed reflectivity, which we will later obtain from DZ and MZ
925 MZ: mask to obtain DZ from CZ
926 NULL the pointers to all other volumes.
928 for (j=0; j<radar->h.nvolumes; j++)
930 if (radar->v[j] == NULL)
932 if ((j == CZ_INDEX) || (j == MZ_INDEX))
935 } /* end for (j=0; ... */
940 /*************************************************************/
942 /* nullGranuleCreate */
944 /*************************************************************/
945 int nullGranuleCreate(char *hdfFileName, IO_HANDLE *granuleHandle,
949 Create an HDF file containing an empty granule.
953 /* Following arrays required by toolkit function 'TKsetL1GVtemplate()'*/
954 int32 TKnparm[MAX_VOS], TKnsweep[MAX_VOS], TKnray[MAX_VOS];
955 int32 TKncell[MAX_VOS][MAX_PARM];
957 /* Check if this HDF file already exists. If it exists, abort. */
958 if (stat(hdfFileName, &buf) == 0)
960 if (radar_verbose_flag)
961 fprintf(stderr, "\nnullGranuleCreate(): File %s already exists.\n",
966 /* Fill toolkit array values for subsequent TKsetL1GVtemplate() call. */
967 TKnparm[0] = TKnsweep[0] = TKnray[0] = TKncell[0][0] = 0;
968 /* Create a L1 GV template_node. */
969 status = TKsetL1GVtemplate(0, TKnparm, TKncell, TKnray, TKnsweep, hdfFileName);
970 if (status != TK_SUCCESS)
972 fprintf(stderr, "nullGranuleCreate(): ***** TKsetL1GVtemplate() error\n");
976 /* Open the HDF file. */
977 if (radar_verbose_flag)
978 fprintf(stderr, "\n\n****** Opening new HDF file: %s for empty granule...\n",
980 status = TKopen(hdfFileName, TK_L1C_GV, TK_NEW_FILE, granuleHandle);
981 if (status != TK_SUCCESS)
983 if (radar_verbose_flag)
984 fprintf(stderr, "nullGranuleCreate(): ***** TKopen() error\n");
987 /* Write metadata fields into HDF file. */
988 metaDataWrite(granuleHandle, radar, hdfFileName, TK_NEW_FILE);
989 /* Close the HDF file */
990 if (radar_verbose_flag)
991 fprintf(stderr, "\n****** Closing HDF file: %s ...\n\n", hdfFileName);
992 status = TKclose(granuleHandle);
993 if (status != TK_SUCCESS)
995 if (radar_verbose_flag)
996 fprintf(stderr, "nullGranuleCreate(): ***** TKclose() error\n");
1003 /*************************************************************/
1005 /* RSL_radar_to_hdf */
1007 /*************************************************************/
1008 int RSL_radar_to_hdf(Radar *radar, char *hdfFileName)
1011 Writes one VOS from a RSL radar structure into one HDF file.
1014 <0 , if failure. (See Error code definitions at top of file.)
1018 L1B_1C_GV *gvl1; /* Toolkit structure for VOS storage. */
1019 IO_HANDLE granuleHandle; /* Toolkit file_descriptor structure. */
1020 VosSize vs; /* Storage of VOS dimensions. */
1021 Volume *v[MAX_RADAR_VOLUMES]; /* Storage of radar volume pointers. */
1023 if (radar == NULL) return(ABORT);
1025 if (radar->h.nvolumes == 0)
1026 /* Create an HDF file to contain an empty granule. */
1028 status = nullGranuleCreate(hdfFileName, &granuleHandle, radar);
1031 /* We will, within functions radarPrep1B51() and radarPrep1C51(),
1032 manipulate the array of radar volume pointers radar->v[].
1033 Hence we here save the array radar->v[] in v[], so we can restore
1034 radar->v[] to its original condition before leaving this function.
1036 radarVolumesSave(v, radar);
1038 /* Get the desired product out of the HDF filename. */
1039 sscanf(strrchr(hdfFileName, '/'), "/%4s", product);
1040 if (strcmp(product, "1B51") == 0)
1042 granuleHandle.productID = TK_L1B_GV;
1043 status = radarPrep1B51(radar);
1045 else if (strcmp(product, "1C51") == 0)
1047 granuleHandle.productID = TK_L1C_GV;
1048 status = radarPrep1C51(radar); /* Creates mask volumes. */
1050 else /* Unknown product. */
1053 if (radar_verbose_flag)
1054 fprintf(stderr, "RSL_radar_to_hdf(): Unknown product requested: %s.\n",
1057 if (status < 0) goto quit;
1059 /* Open the HDF file 'hdfFileName'. */
1060 status = hdfFileOpen(&granuleHandle, &vs, hdfFileName, radar);
1061 if (status < 0) goto quit;
1063 /* Build toolkit 'L1B_1C_GV' structure using data from radar structure.*/
1064 if (radar_verbose_flag)
1065 fprintf(stderr, "\n****** Moving VOS from RSL structure --> toolkit structure ...\n");
1066 gvl1 = gvl1Build(radar, qcParm, &vs, (int)granuleHandle.productID);
1068 /* Write data from toolkit 'L1B_1C_GV' structure to HDF file. */
1069 if (radar_verbose_flag)
1070 fprintf(stderr, "\n****** Writing VOS to HDF file: %s ...\n", hdfFileName);
1071 status = TKwriteL1GV(&granuleHandle, gvl1);
1072 if (status != TK_SUCCESS)
1074 TKclose(&granuleHandle);
1076 if (radar_verbose_flag)
1077 fprintf(stderr, "RSL_radar_to_hdf(): *** TKwriteL1GV() error\n");
1078 goto free_memory_and_quit;
1081 /* Close the HDF file */
1082 if (radar_verbose_flag)
1083 fprintf(stderr, "\n****** Closing HDF file: %s ...\n\n", hdfFileName);
1084 status = TKclose(&granuleHandle);
1085 if (status == TK_SUCCESS)
1089 if (radar_verbose_flag)
1090 fprintf(stderr, "RSL_radar_to_hdf(): *** TKclose() error\n");
1094 free_memory_and_quit:
1095 /* Free memory allocated to the toolkit 'L1B_1C_GV' structure. */
1097 /* If RSL mask volumes for 1C-51 were created above, free them. */
1098 if (radar->v[MZ_INDEX] != NULL)
1099 RSL_free_volume(radar->v[MZ_INDEX]);
1100 if (radar->v[MD_INDEX] != NULL)
1101 RSL_free_volume(radar->v[MD_INDEX]);
1103 /* Restore the array of radar volume pointers in radar->v[]. */
1104 radarVolumesRestore(radar, v);