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.
24 * Ingest a Darwin_Toga file and fill a Radar structure with the data.
25 * Darwin Type 1 data contains 4 fields: corrected & uncorrected refl,
26 * velocity, and spectrum width.
27 * Darwin type 19 data contains 1 field: uncorrected reflectivity.
28 * Ignore all other Darwin data types.
29 * Handle only PPI scans; ignore RHI scans.
31 * Object code from this file must be linked with the following
35 * "libtg.a" handles the reading and decompression of Darwin TOGA
36 * data records from the input toga file.
40 * Applied Research Corp.
54 extern int radar_verbose_flag;
55 int tg_open(char *infile, tg_file_str *tg_file);
56 int tg_read_ray(tg_file_str *tg_file);
58 /* Toga radar routines */
59 void fill_ray(Ray *, tg_file_str *, int);
60 void fill_ray_header(Ray *, tg_file_str *, int, int);
61 void fill_sweep_header(Radar *, tg_map_head_str *, int, int);
62 void fill_volume_header(Radar *, tg_map_head_str *);
63 void fill_radar_header(Radar *, tg_map_head_str *);
68 /********************************************************************/
69 void fill_ray(Ray *ray, tg_file_str *tg_file, int datatype)
70 /* transfer ray data values from the toga ray structure
75 if (ray == NULL) return;
76 for (j=0; j<tg_file->ray.num_bins[datatype]; j++)
78 /* check to see if this is valid data */
79 if (tg_file->ray.data[datatype][j] == TG_NO_DATA)
80 ray->range[j] = ray->h.invf(BADVAL); /* bad data, store flag */
81 else /* good data, store into radar ray */
82 ray->range[j] = ray->h.invf(tg_file->ray.data[datatype][j]);
88 void fill_ray_header(Ray *ray, tg_file_str *tg_file, int elev_num,
90 /* Unresolved: Many of the following values */
92 if (ray == NULL) return;
93 ray->h.month = tg_file->ray_head.mon;
94 ray->h.day = tg_file->ray_head.day;
95 ray->h.year = tg_file->ray_head.year;
96 ray->h.hour = tg_file->ray_head.hour;
97 ray->h.minute = tg_file->ray_head.min;
98 ray->h.sec = tg_file->ray_head.hunsec/100.0;
99 ray->h.unam_rng = MISSING_VAL; /* ?? */
100 ray->h.azimuth = tg_file->ray.azm;
101 ray->h.ray_num = tg_file->ray_num;
102 ray->h.elev = tg_file->ray.elev;
103 ray->h.elev_num = elev_num;
104 ray->h.range_bin1 = (int)(tg_file->ray.start_km[datatype]*1000.0);
105 ray->h.gate_size = (int)(tg_file->ray.interval_km[datatype]*1000.0);
107 ray->h.vel_res = MISSING_VAL; /* ?? */
108 ray->h.sweep_rate = MISSING_VAL; /* ?? */
109 ray->h.prf = tg_file->map_head.prf;
110 ray->h.azim_rate = MISSING_VAL; /* ?? */
111 ray->h.fix_angle = MISSING_VAL; /* ?? */
112 ray->h.pulse_count = MISSING_VAL; /* ?? */
113 ray->h.pulse_width = tg_file->map_head.pulsewd/100.0;
114 ray->h.beam_width = 1.65;
115 ray->h.frequency = 5625;
116 ray->h.wavelength = (tg_file->map_head.wavelen/100.0)/100.0; /* m */
117 ray->h.nyq_vel = (ray->h.wavelength * ray->h.prf)/4.0; /* m/s */
118 ray->h.nbins = tg_file->ray.num_bins[datatype];
120 if (datatype == TG_DM_IND)
123 ray->h.invf = DZ_INVF;
125 else if (datatype == TG_DZ_IND)
128 ray->h.invf = DZ_INVF;
130 else if (datatype == TG_VR_IND)
133 ray->h.invf = VR_INVF;
135 else if (datatype == TG_SW_IND)
138 ray->h.invf = SW_INVF;
144 void fill_sweep_header(Radar *radar, tg_map_head_str *map_head,
145 int sweep_num, int nrays)
146 /* Arrive here _after_ sweep ray data has been filled in. */
148 if (radar_verbose_flag)
149 fprintf(stderr,"sweep_num:%02d num_rays:%d\n",sweep_num, nrays);
151 if (map_head->data_set == 1)
153 radar->v[DZ_INDEX]->sweep[sweep_num]->h.f = DZ_F;
154 radar->v[DZ_INDEX]->sweep[sweep_num]->h.invf = DZ_INVF;
155 radar->v[ZT_INDEX]->sweep[sweep_num]->h.f = DZ_F;
156 radar->v[ZT_INDEX]->sweep[sweep_num]->h.invf = DZ_INVF;
157 radar->v[VR_INDEX]->sweep[sweep_num]->h.f = VR_F;
158 radar->v[VR_INDEX]->sweep[sweep_num]->h.invf = VR_INVF;
159 radar->v[SW_INDEX]->sweep[sweep_num]->h.f = SW_F;
160 radar->v[SW_INDEX]->sweep[sweep_num]->h.invf = SW_INVF;
162 radar->v[DZ_INDEX]->sweep[sweep_num]->h.sweep_num = sweep_num;
163 radar->v[ZT_INDEX]->sweep[sweep_num]->h.sweep_num = sweep_num;
164 radar->v[VR_INDEX]->sweep[sweep_num]->h.sweep_num = sweep_num;
165 radar->v[SW_INDEX]->sweep[sweep_num]->h.sweep_num = sweep_num;
167 radar->v[DZ_INDEX]->sweep[sweep_num]->h.elev =
168 map_head->angfix[sweep_num]/10.0;
169 radar->v[ZT_INDEX]->sweep[sweep_num]->h.elev =
170 map_head->angfix[sweep_num]/10.0;
171 radar->v[VR_INDEX]->sweep[sweep_num]->h.elev =
172 map_head->angfix[sweep_num]/10.0;
173 radar->v[SW_INDEX]->sweep[sweep_num]->h.elev =
174 map_head->angfix[sweep_num]/10.0;
176 radar->v[DZ_INDEX]->sweep[sweep_num]->h.beam_width = 1.65;
177 radar->v[ZT_INDEX]->sweep[sweep_num]->h.beam_width = 1.65;
178 radar->v[VR_INDEX]->sweep[sweep_num]->h.beam_width = 1.65;
179 radar->v[SW_INDEX]->sweep[sweep_num]->h.beam_width = 1.65;
181 radar->v[DZ_INDEX]->sweep[sweep_num]->h.horz_half_bw = 0.85;
182 radar->v[ZT_INDEX]->sweep[sweep_num]->h.horz_half_bw = 0.85;
183 radar->v[VR_INDEX]->sweep[sweep_num]->h.horz_half_bw = 0.85;
184 radar->v[SW_INDEX]->sweep[sweep_num]->h.horz_half_bw = 0.85;
186 radar->v[DZ_INDEX]->sweep[sweep_num]->h.vert_half_bw = 0.85;
187 radar->v[ZT_INDEX]->sweep[sweep_num]->h.vert_half_bw = 0.85;
188 radar->v[VR_INDEX]->sweep[sweep_num]->h.vert_half_bw = 0.85;
189 radar->v[SW_INDEX]->sweep[sweep_num]->h.vert_half_bw = 0.85;
191 radar->v[DZ_INDEX]->sweep[sweep_num]->h.nrays = nrays;
192 radar->v[ZT_INDEX]->sweep[sweep_num]->h.nrays = nrays;
193 radar->v[VR_INDEX]->sweep[sweep_num]->h.nrays = nrays;
194 radar->v[SW_INDEX]->sweep[sweep_num]->h.nrays = nrays;
196 else if (map_head->data_set == 19)
198 radar->v[ZT_INDEX]->sweep[sweep_num]->h.sweep_num = sweep_num;
199 radar->v[ZT_INDEX]->sweep[sweep_num]->h.elev =
200 map_head->angfix[sweep_num]/10.0;
201 radar->v[ZT_INDEX]->sweep[sweep_num]->h.beam_width = 1.65;
202 radar->v[ZT_INDEX]->sweep[sweep_num]->h.horz_half_bw = 0.825;
203 radar->v[ZT_INDEX]->sweep[sweep_num]->h.vert_half_bw = 0.825;
204 radar->v[ZT_INDEX]->sweep[sweep_num]->h.nrays = nrays;
205 radar->v[ZT_INDEX]->sweep[sweep_num]->h.f = DZ_F;
206 radar->v[ZT_INDEX]->sweep[sweep_num]->h.invf = DZ_INVF;
208 else /* Another data type? */
215 void fill_volume_header(Radar *radar, tg_map_head_str *map_head)
217 if (map_head->data_set == 1)
219 radar->v[DZ_INDEX] = RSL_new_volume(MAX_SWEEPS);
220 radar->v[ZT_INDEX] = RSL_new_volume(MAX_SWEEPS);
221 radar->v[VR_INDEX] = RSL_new_volume(MAX_SWEEPS);
222 radar->v[SW_INDEX] = RSL_new_volume(MAX_SWEEPS);
223 radar->v[DZ_INDEX]->h.type_str = strdup("Reflectivity");
224 radar->v[ZT_INDEX]->h.type_str = strdup("Reflectivity");
225 radar->v[VR_INDEX]->h.type_str = strdup("Velocity");
226 radar->v[SW_INDEX]->h.type_str = strdup("Spectrum width");
227 radar->v[DZ_INDEX]->h.nsweeps = map_head->numfix_ang;
228 radar->v[ZT_INDEX]->h.nsweeps = map_head->numfix_ang;
229 radar->v[VR_INDEX]->h.nsweeps = map_head->numfix_ang;
230 radar->v[SW_INDEX]->h.nsweeps = map_head->numfix_ang;
231 radar->v[DZ_INDEX]->h.f = DZ_F;
232 radar->v[DZ_INDEX]->h.invf = DZ_INVF;
233 radar->v[ZT_INDEX]->h.f = DZ_F;
234 radar->v[ZT_INDEX]->h.invf = DZ_INVF;
235 radar->v[VR_INDEX]->h.f = VR_F;
236 radar->v[VR_INDEX]->h.invf = VR_INVF;
237 radar->v[SW_INDEX]->h.f = SW_F;
238 radar->v[SW_INDEX]->h.invf = SW_INVF;
240 else if (map_head->data_set == 19)
242 radar->v[ZT_INDEX] = RSL_new_volume(MAX_SWEEPS);
243 radar->v[ZT_INDEX]->h.type_str = strdup("Reflectivity");
244 radar->v[ZT_INDEX]->h.nsweeps = map_head->numfix_ang;
245 radar->v[ZT_INDEX]->h.f = DZ_F;
246 radar->v[ZT_INDEX]->h.invf = DZ_INVF;
248 else /* Another data type? */
254 void fill_radar_header(Radar *radar, tg_map_head_str *map_head)
256 radar->h.month = map_head->scan_mon;
257 radar->h.day = map_head->scan_day;
258 radar->h.year = map_head->scan_year;
259 radar->h.hour = map_head->scan_hour;
260 radar->h.minute = map_head->scan_min;
261 radar->h.sec = map_head->scan_sec/100.0;
262 strcpy(radar->h.radar_type, "toga");
264 strcpy(radar->h.name, "");
265 strcpy(radar->h.radar_name, "TOGA");
266 strcpy(radar->h.city, "Darwin");
267 strcpy(radar->h.state, "");
268 /* The following lat, lon coordinates obtained from Dave Wolff */
275 radar->h.height = 37;
276 radar->h.spulse = 500; /* ns */
277 radar->h.lpulse = 2000; /* ns */
283 Radar *RSL_toga_to_radar(char *infile)
284 /* Ingest a Darwin_Toga file and fill a Radar structure with the data.
285 Darwin Type 1 data contains 4 fields: corrected & uncorrected refl,
286 velocity, and spectrum width.
287 Darwin type 19 data contains 1 field: uncorrected reflectivity.
288 Ignore all other Darwin data types.
289 Handle only PPI scans; ignore RHI scans. */
291 int m, swp_num, num_bins;
295 extern int *rsl_qsweep; /* See RSL_read_these_sweeps in volume.c */
296 extern int rsl_qsweep_max;
298 /* open the toga data file and read toga file header into the
299 tg_file map_head_structure */
300 if (tg_open(infile,&tg_file) < 0) /* Understands NULL 'infile' as stdin. */
302 if (radar_verbose_flag)
303 fprintf(stderr,"Error opening/reading data file\n");
307 if (radar_verbose_flag)
309 fprintf(stderr,"Input file: %s\n",infile);
310 fprintf(stderr,"Scan_date: %.2d/%.2d/%d\n",tg_file.map_head.scan_mon,
311 tg_file.map_head.scan_day,tg_file.map_head.scan_year);
312 fprintf(stderr,"Scan_time: %.2d:%.2d\n",tg_file.map_head.scan_hour,
313 tg_file.map_head.scan_min);
315 /* If this toga file does not contain a PPI (constant elevation)
316 scan, do not bother with it. Just quit. */
317 if (tg_file.map_head.scanmod != 1)
319 if (radar_verbose_flag)
320 fprintf(stderr,"Darwtoga file %s does not contain a PPI scan.\n",infile);
324 /* Can handle only datatypes 1 and 19 */
325 if ((tg_file.map_head.data_set != 1) &&
326 (tg_file.map_head.data_set != 19))
328 if (radar_verbose_flag)
329 fprintf(stderr,"File %s does not contain Doppler or refl data\n",infile);
333 if (radar_verbose_flag)
335 if (tg_file.map_head.data_set == 1) fprintf(stderr,"Type 1 data\n");
336 else if (tg_file.map_head.data_set == 19) fprintf(stderr,"Type 19 data\n");
337 else fprintf(stderr,"Unknown data type\n");
340 radar = RSL_new_radar(MAX_RADAR_VOLUMES);
341 fill_radar_header(radar, &tg_file.map_head);
342 fill_volume_header(radar, &tg_file.map_head);
344 /* initialize counters */
345 tg_file.ray_num = -1;
348 /* Main loop to read in a toga ray from the toga
349 data file and store into radar structure */
350 while ((m=tg_read_ray(&tg_file)) > 0)
352 /* check for end_of_sweep. Fill the sweep header _after_ we've
354 if (swp_num < tg_file.ray_head.tilt - 1) /* new sweep? */
358 /* fill the sweep header for the previous sweep */
359 fill_sweep_header(radar, &tg_file.map_head, swp_num,
362 tg_file.ray_num = -1; /* Reset ray_num. */
363 swp_num += 1; /* increment sweep count */
364 if (rsl_qsweep != NULL) {
365 if (swp_num > rsl_qsweep_max) break;
366 if (rsl_qsweep[swp_num] == 0) continue;
368 /* Check for too many sweeps. */
369 if ((tg_file.map_head.numfix_ang < swp_num + 1) ||
370 (MAX_SWEEPS < swp_num + 1))
372 perror("toga_to_radar: Exceeded expected no. of sweeps");
374 RSL_free_radar(radar);
377 /* Create new sweep structures. */
378 if (tg_file.map_head.data_set == 1)
380 radar->v[DZ_INDEX]->sweep[swp_num] = RSL_new_sweep(MAX_RAYS);
381 radar->v[ZT_INDEX]->sweep[swp_num] = RSL_new_sweep(MAX_RAYS);
382 radar->v[VR_INDEX]->sweep[swp_num] = RSL_new_sweep(MAX_RAYS);
383 radar->v[SW_INDEX]->sweep[swp_num] = RSL_new_sweep(MAX_RAYS);
385 else if (tg_file.map_head.data_set == 19)
387 radar->v[ZT_INDEX]->sweep[swp_num] = RSL_new_sweep(MAX_RAYS);
389 } /* end if new sweep */
391 num_bins = tg_file.ray.num_bins[TG_DM_IND] + 4; /* types 1,19 data*/
392 tg_file.ray_num += 1;
393 /* Check for too many rays. */
394 if (tg_file.ray_num > MAX_RAYS)
396 perror("toga_to_radar: Exceeded maximal no. of rays");
398 RSL_free_radar(radar);
401 /* check for uncorrected reflectivity data */
402 if (tg_file.ray.da_inv[TG_DM_IND] == TRUE)
404 new_ray = RSL_new_ray(num_bins);
405 radar->v[ZT_INDEX]->sweep[swp_num]->ray[tg_file.ray_num] = new_ray;
406 fill_ray_header(new_ray, &tg_file, swp_num, TG_DM_IND);
407 fill_ray(new_ray, &tg_file, TG_DM_IND);
409 /* check for corrected reflectivity data */
410 if (tg_file.ray.da_inv[TG_DZ_IND] == TRUE)
412 new_ray = RSL_new_ray(num_bins);
413 radar->v[DZ_INDEX]->sweep[swp_num]->ray[tg_file.ray_num] = new_ray;
414 fill_ray_header(new_ray, &tg_file, swp_num, TG_DZ_IND);
415 fill_ray(new_ray, &tg_file, TG_DZ_IND);
417 /* check for velocity data */
418 if (tg_file.ray.da_inv[TG_VR_IND] == TRUE)
420 new_ray = RSL_new_ray(num_bins);
421 radar->v[VR_INDEX]->sweep[swp_num]->ray[tg_file.ray_num] = new_ray;
422 fill_ray_header(new_ray, &tg_file, swp_num, TG_VR_IND);
423 fill_ray(new_ray, &tg_file, TG_VR_IND);
425 /* check for spectrum width data */
426 if (tg_file.ray.da_inv[TG_SW_IND] == TRUE)
428 new_ray = RSL_new_ray(num_bins);
429 radar->v[SW_INDEX]->sweep[swp_num]->ray[tg_file.ray_num] = new_ray;
430 fill_ray_header(new_ray, &tg_file, swp_num, TG_SW_IND);
431 fill_ray(new_ray, &tg_file, TG_SW_IND);
435 /* At this point, we are finished reading toga data records.
436 Fill in the last sweep header. */
437 fill_sweep_header(radar, &tg_file.map_head, swp_num,
440 /* Check which flag the TOGA read_record routines returned, and
441 print appropriate terminating message */
442 if (radar_verbose_flag)
447 fprintf(stderr,"Reached end of data file\n");
449 case TG_RAY_READ_ERR:
450 fprintf(stderr,"Error reading toga ray\n");
453 fprintf(stderr,"Found unrecognized toga ray type\n");
454 fprintf(stderr,"ray_type: %d\n",tg_file.ray_head.type);
457 fprintf(stderr,"found out-of-sequence toga record \n");
460 fprintf(stderr,"Error reading toga ray_header \n");
463 fprintf(stderr,"Error reading toga file \n");
466 } /* end if (radar_verbose_flag) */
469 if (m == TG_END_DATA) {
470 radar = RSL_prune_radar(radar);
474 RSL_free_radar(radar);