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.
23 /*******************************************************************
24 * Ingest a Mcgill format file and fill a RSL Radar structure
27 * Object code from this file must be linked with the following libraries:
30 *-----------------------------------------------------------------
31 * NOTE: Adjust Mcgill dbz bias in function RayFill() .
32 * Presently, if (Mcgill_dbz_value != 0)
33 * then Mcgill_dbz_value + 16.5 is stored in RSL Ray.
34 * if (Mcgill dbz value == 0)
35 * then 0 is stored in RSL Ray.
36 * I don't find the Mcgill documentation very clear in this matter.
37 *-----------------------------------------------------------------
38 * Presently sets radar_head->radar_type = DARWIN_TOGA_SIGMET
39 * (Need a radar_type defined in rsl.h for PAFB or MCGILL)
40 *----------------------------------------------------------------
41 * Minor changes will be required in the associated Mcgill library to
42 * accomodate Mcgill data from Sao Paulo. Since I don't have any files
43 * from Sao Paulo to test, I haven't incorporated any such changes.
44 *-------------------------------------------------------------------
45 * Functions defined in this file:
47 * void RayFill(Ray *rsl_ray, mcgRay_t *mcg_ray);
48 * void Ray_headerInit(Ray *ray, mcgHeader_t *head,
49 * mcgRay_t *mcg_ray, int ray_num, int num_bins_rsl);
50 * void Sweep_headerInit(Sweep *sweep, mcgRay_t *mcg_ray,
52 * void Volume_headerInit(Volume *volume, short vol_scan_format);
53 * void Radar_headerInit(Radar *radar, mcgHeader_t *mcg_head);
54 * Radar *RSL_mcgill_to_radar(char *infile);
59 *******************************************************************/
69 #define MCG_DBZ_BIAS 16.5
70 #define MCG_NOISE_BIAS 0.0
72 extern int radar_verbose_flag;
74 /*********************** Function Prototypes ***************************/
75 static void RayFill(Ray *rsl_ray, mcgRay_t *mcg_ray);
76 static void Ray_headerInit(Ray *ray, mcgHeader_t *head,
77 mcgRay_t *mcg_ray, int ray_num, int num_bins_rsl);
78 static void Sweep_headerInit(Sweep *sweep, mcgRay_t *mcg_ray, int nrays);
79 static void Volume_headerInit(Volume *volume, short vol_scan_format);
80 static void Radar_headerInit(Radar *radar, mcgHeader_t *mcg_head);
81 Radar *RSL_mcgill_to_radar(char *infile);
83 static float (*f)(Range x);
84 static Range (*invf)(float x);
86 /* Fixed number_of_RSL_bins for Mcgill sweeps. Note that
87 the number of bins in the RSL structure differs from
88 the number of bins in the Mcgill ray stucture, since
89 the km/bin spacing differs. SEE MCGILL DOCUMENTATION */
90 static int num_bins_normal[24] =
92 /* PAFB volume scan formats 1 and 3 (normal mode) */
93 480, 240, 240, 240, 240, 240, 240, 240, 240, 240,
94 240, 240, 240, 240, 240, 240, 240, 240, 240, 240,
97 static int num_bins_compressed[24] =
99 /* PAFB volume scan formats 2 and 4 (compressed mode) */
100 240, 120, 120, 120, 120, 120, 120, 120, 120, 120,
101 120, 120, 120, 120, 120, 120, 120, 120, 120, 120,
106 /********************************************************************/
107 void RayFill(Ray *rsl_ray, mcgRay_t *mcg_ray)
108 /********************************************************************/
109 /* transfer ray data values from the mcgRay structure
110 to the radar->v[]->sweep[]->ray[] structure */
114 int index=1; /* Used to accomodate varying Mcgill bin spacing */
116 if (rsl_ray == NULL) return;
118 /* Note that the number of bins in the Mcgill ray structure
119 is different from the number of bins in the RSL ray structure,
120 since the Mcgill bin spacing changes with range. RSL ray bin
122 ---------------------------------------------------------
123 MCGILL NORMAL FORMAT:
124 The first 120 mcg_bins have a length of 1 km, and span the
125 the range [0, 120] km.
126 Mcg_bins 120 to 180 have a length of 2 km, and span the
128 Mcg_bins 180 to 240 have a length of 4 km, and span the
130 ---------------------------------------------------------
131 MCGILL COMPRESSED FORMAT:
132 The first 120 mcg_bins have a length of 2 km, and span the
134 Mcg_bins 120 to 180 have a length of 4 km, and span the
137 for (mcg_bin=0; mcg_bin<mcg_ray->num_bins; mcg_bin++)
139 /* I do the case dbz_value of zero separately, so produced
140 color images have a black background. Otherwise, if >>all<<
141 dbz values are biased by 16.5, we get a colored background
142 where there is no precipitation. Looks peculiar. Somebody
143 knowledgable should look into this matter. */
144 if (mcg_ray->data[mcg_bin] == 0)
146 for (j=0; j<index; j++)
148 rsl_ray->range[rsl_bin] = invf(MCG_NOISE_BIAS);
154 for (j=0; j<index; j++)
156 rsl_ray->range[rsl_bin] = invf(mcg_ray->data[mcg_bin] +
161 /* Mcgill bin spacing changes at bins 120 and 180 */
162 if ((mcg_bin == 119) || (mcg_bin == 179))
164 /* In normal format, put range rings at 120km and 240km.
165 In compressed format, put range rings at 240km and 480km */
166 /* rsl_ray->range[rsl_bin-1] = RSL_INVF(60.0); */
167 index = index*2; /* Set up new bin spacing index */
169 } /* end for (mcg_bin=0...*/
173 /*************************************************************************/
174 void Ray_headerInit(Ray *ray, mcgHeader_t *head, mcgRay_t *mcg_ray,
175 int ray_num, int num_bins_rsl)
176 /*************************************************************************/
178 if (ray == NULL) return;
179 ray->h.month = (int)head->month;
180 ray->h.day = (int)head->day;
181 ray->h.year = 1900 + (int)head->year;
182 if (ray->h.year < 1980) ray->h.year += 100; /* Year > 2000 */
183 ray->h.hour = (int)head->hour;
184 ray->h.minute = (int)head->min;
185 ray->h.sec = (float)head->sec;
186 ray->h.unam_rng = MISSING_VAL; /***** ?? ******/
187 ray->h.azimuth = mcg_ray->azm;
188 ray->h.ray_num = ray_num;
189 ray->h.elev = mcg_ray->elev;
190 ray->h.elev_num = mcg_ray->sweep_num;
191 ray->h.range_bin1 = 0;
192 switch (head->vol_scan_format)
195 ray->h.gate_size = 1000; /* 1 km/bin */
196 ray->h.fix_angle = 24; /* 24 sweeps */
199 ray->h.gate_size = 2000; /* 2 km/bin */
200 ray->h.fix_angle = 24; /* 24 sweeps */
203 ray->h.gate_size = 1000; /* 1 km/bin */
204 ray->h.fix_angle = 12; /* 12 sweeps */
207 ray->h.gate_size = 2000; /* 2 km/bin */
208 ray->h.fix_angle = 12; /* 12 sweeps */
213 ray->h.vel_res = MISSING_VAL; /* ?? */
214 ray->h.sweep_rate = MISSING_VAL; /* ?? */
215 ray->h.prf = MISSING_VAL;
216 ray->h.azim_rate = MISSING_VAL; /* ?? */
217 ray->h.pulse_count = MISSING_VAL; /* ?? */
218 ray->h.pulse_width = MISSING_VAL; /* ?? */
219 ray->h.frequency = MISSING_VAL; /* ?? */
220 ray->h.wavelength = MISSING_VAL; /* ?? */
221 ray->h.nyq_vel = MISSING_VAL; /* ?? */
222 ray->h.nbins = num_bins_rsl;
230 /*********************************************************************/
231 void Sweep_headerInit(Sweep *sweep, mcgRay_t *mcg_ray, int nrays)
232 /*********************************************************************/
233 /* Arrive here _after_ sweep ray data has been filled in. */
235 if (radar_verbose_flag)
236 fprintf(stderr,"sweep_num:%02d num_rays:%d\n",mcg_ray->sweep_num, nrays);
237 if (sweep == NULL) return;
238 sweep->h.sweep_num = mcg_ray->sweep_num;
239 sweep->h.elev = mcg_ray->elev;
240 sweep->h.beam_width = 2.0; /* From Dennis' function mcg2uf().
241 I don't know where he got it. */
242 sweep->h.horz_half_bw = 1.0;
243 sweep->h.vert_half_bw = 1.0;
244 sweep->h.nrays = nrays;
246 sweep->h.invf = invf;
251 /************************************************************************/
252 void Volume_headerInit(Volume *volume, short vol_scan_format)
253 /************************************************************************/
255 if (volume == NULL) return;
256 volume->h.type_str = strdup("Reflectivity");
258 volume->h.invf = invf;
260 switch (vol_scan_format)
263 volume->h.nsweeps = 24;
266 volume->h.nsweeps = 12;
274 /***********************************************************************/
275 void Radar_headerInit(Radar *radar, mcgHeader_t *mcg_head)
276 /***********************************************************************/
278 radar->h.month = (int)mcg_head->month;
279 radar->h.day = (int)mcg_head->day;
280 radar->h.year = 1900 + (int)mcg_head->year;
281 radar->h.hour = (int)mcg_head->hour;
282 radar->h.minute = (int)mcg_head->min;
283 radar->h.sec = (float)mcg_head->sec;
284 strcpy(radar->h.radar_type, "mcgill");
285 radar->h.nvolumes = 1; /* Mcgill contains only refl data */
286 radar->h.number = MISSING_VAL; /****************/
287 strcpy(radar->h.name, "PAFB");
288 strcpy(radar->h.radar_name, "MCGILL");
289 strcpy(radar->h.city, "MELB");
290 strcpy(radar->h.state, "FL");
298 radar->h.spulse = MISSING_VAL; /* ns */ /*************/
299 radar->h.lpulse = MISSING_VAL; /* ns */ /*************/
304 /*********************************************************************/
305 Radar *RSL_mcgill_to_radar(char *infile)
306 /*********************************************************************/
308 /* Ingest a Mcgill format radar data file and fill a Radar RSL
309 structure with the data.
315 mcgRay_t *mcg_ray, *mcg_ray_last, *swap;
316 extern int *rsl_qsweep; /* See RSL_read_these_sweeps in volume.c */
317 extern int rsl_qsweep_max;
319 /* If no filename has been passed, there's nothing to do. */
323 /* Default conversion functions. */
327 /* Create a structure of type Radar */
328 radar = (Radar *)RSL_new_radar(MAX_RADAR_VOLUMES);
330 perror("RSL_mcgill_to_radar: RSL_new_radar:");
335 /* Open the Mcgill data file and read Mcgill file header into the
336 mcgFile.head structure */
337 file = (mcgFile_t *)mcgFileOpen(&code, infile);
341 if (radar_verbose_flag)
343 fprintf(stderr,"Input file: %s\n", infile);
344 fprintf(stderr,"Scan_date: %d/%d/%d\n", file->head.month, file->head.day,
346 fprintf(stderr,"Scan_time: %d:%d:%d\n", file->head.hour, file->head.min,
348 fprintf(stderr,"num_recs:%d format:%d\n\n",
349 file->head.num_records, file->head.vol_scan_format);
352 /* Initialize num_bins_rsl to point to the proper bin count array for
354 if ((file->head.vol_scan_format == 1) ||
355 (file->head.vol_scan_format == 3)) /* Normal mode */
357 num_bins_rsl = num_bins_normal;
359 else /* Compressed mode */
360 num_bins_rsl = num_bins_compressed;
362 /* Allocate space for two mcg_rays. mcg_ray will contain the current
363 ray, and mcg_ray_last will contain the previous ray. The previous ray
364 is needed only at sweep increments. */
365 mcg_ray = (mcgRay_t *)malloc(sizeof(mcgRay_t));
366 mcg_ray_last = (mcgRay_t *)malloc(sizeof(mcgRay_t));
368 /* Initialize the values of the Radar_header structure */
369 Radar_headerInit(radar, &file->head);
370 /* Create a structure of type Volume */
371 radar->v[DZ_INDEX] = RSL_new_volume(MAX_SWEEPS);
372 if (radar->v[DZ_INDEX] == NULL) {
373 perror("mcgill_to_radar: RSL_new_volume");
376 /* Initialize the values of the Volume_header structure. */
377 Volume_headerInit(radar->v[DZ_INDEX], file->head.vol_scan_format);
379 /* initialize counters */
381 mcg_ray_last->sweep_num = -1;
383 /* Main loop to read in a ray from the Mcgill data file and store
384 into RSL radar structure */
385 while ((code = mcgRayBuild(mcg_ray, file)) == MCG_OK)
387 /* Discard rays with bogus azimuth values. */
388 if (mcg_ray->azm > 360.0)
390 if (radar_verbose_flag)
391 fprintf(stderr,"**** Bogus azm:%.1f, discarding ray.\n", mcg_ray->azm);
394 /* Check for end_of_sweep. Fill the sweep header _after_ we've
395 read in the sweep, since we don't know the number of
396 rays in the sweep until we find the start of the next sweep*/
397 if (mcg_ray_last->sweep_num < mcg_ray->sweep_num) /* new sweep? */
399 if (mcg_ray_last->sweep_num >= 0)
401 /* fill the sweep header for the previous sweep */
402 Sweep_headerInit(radar->v[DZ_INDEX]->sweep[mcg_ray_last->sweep_num],
403 mcg_ray_last, ray_num+1);
405 ray_num = -1; /* Reset ray counter at start of each sweep */
407 /* Check for too many sweeps. If we've exceeded MAX_SWEEPS,
409 if (MAX_SWEEPS < mcg_ray->sweep_num + 1)
411 perror("RSL_mcgill_to_radar: Exceeded expected no. of sweeps");
413 RSL_free_radar(radar);
416 if (rsl_qsweep != NULL) {
417 if (mcg_ray->sweep_num > rsl_qsweep_max) break;
418 if (rsl_qsweep[mcg_ray->sweep_num] == 0) continue;
421 /* Create new sweep structure. */
422 radar->v[DZ_INDEX]->sweep[mcg_ray->sweep_num] = RSL_new_sweep(MAX_RAYS);
423 if (radar->v[DZ_INDEX]->sweep[mcg_ray->sweep_num] == NULL) {
424 perror("mcgill_to_radar: RSL_new_sweep");
427 } /* end if new sweep */
429 ray_num += 1; /* Increment ray counter */
430 /* Check for too many rays. */
431 if (ray_num > MAX_RAYS)
433 perror("mcgill_to_radar: Exceeded maximal no. of rays");
435 RSL_free_radar(radar);
439 /* Create a new RSL ray containing correct # of bins, and fill it. */
440 radar->v[DZ_INDEX]->sweep[mcg_ray->sweep_num]->ray[ray_num] =
441 RSL_new_ray(num_bins_rsl[mcg_ray->sweep_num] + 8);
442 if (radar->v[DZ_INDEX]->sweep[mcg_ray->sweep_num]->ray[ray_num] == NULL ) {
443 perror("mcgill_to_radar: RSL_new_ray");
447 Ray_headerInit(radar->v[DZ_INDEX]->sweep[mcg_ray->sweep_num]->ray[ray_num],
448 &file->head, mcg_ray, ray_num, num_bins_rsl[mcg_ray->sweep_num]);
449 RayFill(radar->v[DZ_INDEX]->sweep[mcg_ray->sweep_num]->ray[ray_num],
452 /* Swap mcg_ray pointers so that the current mcg_ray, which we've
453 finished loading into an RSL_ray, becomes mcg_ray_last . */
454 swap = (mcgRay_t *)mcg_ray_last;
455 mcg_ray_last = (mcgRay_t *)mcg_ray;
456 mcg_ray = (mcgRay_t *)swap;
459 /* At this point, we're finished reading file records.
460 Fill in the last sweep header. */
461 Sweep_headerInit(radar->v[DZ_INDEX]->sweep[mcg_ray_last->sweep_num],
462 mcg_ray_last, ray_num+1);
464 /* Check which flag the mcgill routines returned, and
465 print appropriate terminating message */
466 quit: if (radar_verbose_flag)
471 fprintf(stderr,"MCG_EOD: Reached end of data file\n");
473 case MCG_OPEN_FILE_ERR:
474 fprintf(stderr,"MCG_OPEN_FILE_ERR: Error opening Mcgill data file\n");
477 fprintf(stderr,"MCG_EOF: End of data file reached prematurely\n");
480 fprintf(stderr,"MCG_READ_ERR: Error occurred while reading from data file\n");
483 fprintf(stderr,"MCG_FORMAT_ERR: Format error encountered in data file\n");
486 fprintf(stderr,"Error reading data file \n");
489 } /* end if (radar_verbose_flag) */
491 if (code == MCG_EOD) /* Successfully read in Mcgill file? */
494 radar = RSL_prune_radar(radar);
497 else /* Fatal error occurred */
499 RSL_free_radar(radar);