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.
31 extern int radar_verbose_flag;
33 * These externals can be found in the wsr88d library; secret code.
35 void print_head(Wsr88d_file_header wsr88d_file_header);
36 void clear_sweep(Wsr88d_sweep *wsr88d_sweep, int x, int n);
37 void free_and_clear_sweep(Wsr88d_sweep *wsr88d_sweep, int x, int n);
39 * Secretly in uf_to_radar.c
41 Volume *copy_sweeps_into_volume(Volume *new_volume, Volume *old_volume);
44 void float_to_range(float *x, Range *c, int n, Range (*function)(float x) )
47 if (*x == WSR88D_BADVAL) *c = function(BADVAL);
48 else if (*x == WSR88D_RFVAL) *c = function(RFVAL);
49 else *c = function(*x);
54 /**********************************************************************/
56 /* done 3/23 wsr88d_load_sweep_into_volume */
58 /* By: John Merritt */
59 /* Space Applications Corporation */
61 /**********************************************************************/
62 int wsr88d_load_sweep_into_volume(Wsr88d_sweep ws,
63 Volume *v, int nsweep, unsigned int vmask)
77 Range (*invf)(float x);
80 /* Allocate memory for MAX_RAYS_IN_SWEEP rays. */
81 v->sweep[nsweep] = RSL_new_sweep(MAX_RAYS_IN_SWEEP);
82 if (v->sweep[nsweep] == NULL) {
83 perror("wsr88d_load_sweep_into_volume: RSL_new_sweep");
87 v->sweep[nsweep]->h.elev = 0;
88 v->sweep[nsweep]->h.nrays = 0;
89 f = (float (*)(Range x))NULL;
90 invf = (Range (*)(float x))NULL;
91 if (vmask & WSR88D_DZ) { invf = DZ_INVF; f = DZ_F; }
92 if (vmask & WSR88D_VR) { invf = VR_INVF; f = VR_F; }
93 if (vmask & WSR88D_SW) { invf = SW_INVF; f = SW_F; }
97 v->sweep[nsweep]->h.invf = invf;
98 v->sweep[nsweep]->h.f = f;
100 for (i=0,iray=0; i<MAX_RAYS_IN_SWEEP; i++) {
101 if (ws.ray[i] != NULL) {
102 wsr88d_ray_to_float(ws.ray[i], vmask, v_data, &n);
103 float_to_range(v_data, c_data, n, invf);
105 wsr88d_get_date(ws.ray[i], &mon, &day, &year);
106 wsr88d_get_time(ws.ray[i], &hh, &mm, &ss, &fsec);
108 fprintf(stderr,"n %d, mon %d, day %d, year %d, hour %d, min %d, sec %d, fsec %f\n",
109 n, mon, day, year, hh, mm, ss, fsec);
112 * Load the sweep/ray headar information.
115 v->sweep[nsweep]->ray[iray] = RSL_new_ray(n);
116 /*(Range *)calloc(n, sizeof(Range)); */
118 ray_ptr = v->sweep[nsweep]->ray[iray]; /* Make code below readable. */
120 ray_ptr->h.invf = invf;
121 ray_ptr->h.month = mon;
122 ray_ptr->h.day = day;
123 ray_ptr->h.year = year + 1900; /* Yes 1900 makes this year 2000 compliant, due to wsr88d using unix time(). */
124 ray_ptr->h.hour = hh;
125 ray_ptr->h.minute = mm;
126 ray_ptr->h.sec = ss + fsec;
127 ray_ptr->h.unam_rng = wsr88d_get_range (ws.ray[i]);
128 ray_ptr->h.azimuth = wsr88d_get_azimuth (ws.ray[i]);
129 /* -180 to +180 is converted to 0 to 360 */
130 if (ray_ptr->h.azimuth < 0) ray_ptr->h.azimuth += 360;
131 ray_ptr->h.ray_num = ws.ray[i]->ray_num;
132 ray_ptr->h.elev = wsr88d_get_elevation_angle(ws.ray[i]);
133 ray_ptr->h.elev_num = ws.ray[i]->elev_num;
134 if (vmask & WSR88D_DZ) {
135 ray_ptr->h.range_bin1 = ws.ray[i]->refl_rng;
136 ray_ptr->h.gate_size = ws.ray[i]->refl_size;
138 ray_ptr->h.range_bin1 = ws.ray[i]->dop_rng;
139 ray_ptr->h.gate_size = ws.ray[i]->dop_size;
141 ray_ptr->h.vel_res = wsr88d_get_velocity_resolution(ws.ray[i]);
142 vol_cpat = wsr88d_get_volume_coverage(ws.ray[i]);
144 case 11: ray_ptr->h.sweep_rate = 16.0/5.0; break;
145 case 12: ray_ptr->h.sweep_rate = 17.0/4.2; break;
146 case 21: ray_ptr->h.sweep_rate = 11.0/6.0; break;
147 case 31: ray_ptr->h.sweep_rate = 8.0/10.0; break;
148 case 32: ray_ptr->h.sweep_rate = 7.0/10.0; break;
149 case 121:ray_ptr->h.sweep_rate = 20.0/5.5; break;
150 default: ray_ptr->h.sweep_rate = 0.0; break;
153 ray_ptr->h.nyq_vel = wsr88d_get_nyquist(ws.ray[i]);
154 ray_ptr->h.azim_rate = wsr88d_get_azimuth_rate(ws.ray[i]);
155 ray_ptr->h.fix_angle = wsr88d_get_fix_angle(ws.ray[i]);
156 ray_ptr->h.pulse_count = wsr88d_get_pulse_count(ws.ray[i]);
157 ray_ptr->h.pulse_width = wsr88d_get_pulse_width(ws.ray[i]);
158 ray_ptr->h.beam_width = .95;
159 ray_ptr->h.prf = wsr88d_get_prf(ws.ray[i]);
160 ray_ptr->h.frequency = wsr88d_get_frequency(ws.ray[i]);
161 ray_ptr->h.wavelength = 0.1071; /* Previously called
162 * wsr88d_get_wavelength(ws.ray[i]).
163 * See wsr88d.c for explanation.
166 /* It is no coincidence that the 'vmask' and wsr88d datatype
167 * values are the same. We expect 'vmask' to be one of
168 * REFL_MASK, VEL_MASK, or SW_MASK. These match WSR88D_DZ,
169 * WSR88D_VR, and WSR88D_SW in the wsr88d library.
171 ray_ptr->h.nbins = n;
172 memcpy(ray_ptr->range, c_data, n*sizeof(Range));
173 v->sweep[nsweep]->h.nrays = iray+1;
174 v->sweep[nsweep]->h.elev += ray_ptr->h.elev;
175 v->sweep[nsweep]->h.sweep_num = ray_ptr->h.elev_num;
180 v->sweep[nsweep]->h.beam_width = .95;
181 v->sweep[nsweep]->h.vert_half_bw = .475;
182 v->sweep[nsweep]->h.horz_half_bw = .475;
183 /* Now calculate the mean elevation angle for this sweep. */
184 if (v->sweep[nsweep]->h.nrays > 0)
185 v->sweep[nsweep]->h.elev /= v->sweep[nsweep]->h.nrays;
187 RSL_free_sweep(v->sweep[nsweep]); /* No rays loaded, free this sweep. */
188 v->sweep[nsweep] = NULL;
194 /**********************************************************************/
196 /* RSL_wsr88d_to_radar */
198 /* By: John Merritt */
199 /* Space Applications Corporation */
201 /**********************************************************************/
203 Radar *RSL_wsr88d_to_radar(char *infile, char *call_or_first_tape_file)
205 * Gets all volumes from the nexrad file. Input file is 'infile'.
206 * Site information is extracted from 'call_or_first_tape_file'; this
207 * is typically a disk file called 'nex.file.1'.
211 * Uses the string in 'call_or_first_tape_file' as the 4 character call sign
212 * for the sight. All UPPERCASE characters. Normally, this call sign
213 * is extracted from the file 'nex.file.1'.
215 * Returns a pointer to a Radar structure; that contains the different
222 Wsr88d_sweep wsr88d_sweep;
223 Wsr88d_file_header wsr88d_file_header;
224 Wsr88d_tape_header wsr88d_tape_header;
230 int volume_mask[] = {WSR88D_DZ, WSR88D_VR, WSR88D_SW};
231 char *field_str[] = {"Reflectivity", "Velocity", "Spectrum width"};
232 Wsr88d_site_info *sitep;
235 int expected_msgtype = 0;
238 extern int rsl_qfield[]; /* See RSL_select_fields in volume.c */
239 extern int *rsl_qsweep; /* See RSL_read_these_sweeps in volume.c */
240 extern int rsl_qsweep_max;
242 Radar *wsr88d_load_m31_into_radar(Wsr88d_file *wf);
245 /* Determine the site quasi automatically. Here is the procedure:
246 * 1. Determine if we have a call sign.
247 * 2. Try reading 'call_or_first_tape_file' from disk. This is done via
248 * wsr88d_read_tape_header.
249 * 3. If no valid site info, abort.
251 if (call_or_first_tape_file == NULL) {
252 fprintf(stderr, "wsr88d_to_radar: No valid site ID info provided.\n");
254 } else if (strlen(call_or_first_tape_file) == 4)
255 sitep = wsr88d_get_site(call_or_first_tape_file);
256 else if (strlen(call_or_first_tape_file) == 0) {
257 fprintf(stderr, "wsr88d_to_radar: No valid site ID info provided.\n");
262 if (wsr88d_read_tape_header(call_or_first_tape_file, &wsr88d_tape_header) > 0) {
263 memcpy(site_id_str, wsr88d_tape_header.site_id, 4);
264 sitep = wsr88d_get_site(site_id_str);
267 fprintf(stderr,"wsr88d_to_radar: No valid site ID info found.\n");
270 if (radar_verbose_flag)
271 fprintf(stderr,"SITE: %c%c%c%c\n", sitep->name[0], sitep->name[1],
272 sitep->name[2], sitep->name[3]);
275 memset(&wsr88d_sweep, 0, sizeof(Wsr88d_sweep)); /* Initialize to 0 a
276 * heavily used variable.
279 /* 1. Open the input wsr88d file. */
280 if (infile == NULL) the_file = "stdin"; /* wsr88d.c understands this to
281 * mean read from stdin.
283 else the_file = infile;
285 if ((wf = wsr88d_open(the_file)) == NULL) {
286 wsr88d_perror(the_file);
292 /* 2. Read wsr88d headers. */
293 /* Return # bytes, 0 or neg. on fail. */
294 n = wsr88d_read_file_header(wf, &wsr88d_file_header);
296 * Get the expected digital radar message type based on version string
297 * from the Archive II header. The message type is 31 for Build 10, and 1
298 * for prior builds. Note that we consider AR2V0001 to be message type 1,
299 * because it has been in the past, but with Build 10 this officially
300 * becomes the version number for Evansville (KVWX), which will use message
301 * type 31. This could be a problem if RSL is used to process KVWX.
304 strncpy(version, wsr88d_file_header.title.filename, 8);
305 if (strncmp(version,"AR2V0007",8) == 0 ||
306 strncmp(version,"AR2V0006",8) == 0 ||
307 strncmp(version,"AR2V0004",8) == 0 ||
308 strncmp(version,"AR2V0003",8) == 0 ||
309 strncmp(version,"AR2V0002",8) == 0) {
310 expected_msgtype = 31;
312 else if (strncmp(version,"ARCHIVE2",8) == 0 ||
313 strncmp(version,"AR2V0001",8) == 0) {
314 expected_msgtype = 1;
318 if (n <= 0 || expected_msgtype == 0) {
319 fprintf(stderr,"RSL_wsr88d_to_radar: ");
321 fprintf(stderr,"wsr88d_read_file_header failed\n");
323 fprintf(stderr,"Archive II header contains unknown version "
324 ": '%s'\n", version);
329 if (radar_verbose_flag)
330 print_head(wsr88d_file_header);
333 if (expected_msgtype == 31) {
335 /* Get radar for message type 31. */
337 radar = wsr88d_load_m31_into_radar(wf);
338 if (radar == NULL) return NULL;
341 /* Get radar for message type 1. */
343 /* Allocate all Volume pointers. */
344 radar = RSL_new_radar(MAX_RADAR_VOLUMES);
345 if (radar == NULL) return NULL;
347 /* Clear the sweep pointers. */
348 clear_sweep(&wsr88d_sweep, 0, MAX_RAYS_IN_SWEEP);
350 /* Allocate a maximum of 30 sweeps for the volume. */
351 /* Order is important. WSR88D_DZ, WSR88D_VR, WSR88D_SW, is
352 * assigned to the indexes DZ_INDEX, VR_INDEX and SW_INDEX respectively.
355 for (iv=0; iv<nvolumes; iv++)
356 if (rsl_qfield[iv]) radar->v[iv] = RSL_new_volume(20);
361 for (;(n = wsr88d_read_sweep(wf, &wsr88d_sweep)) > 0; nsweep++) {
362 if (rsl_qsweep != NULL) {
363 if (nsweep > rsl_qsweep_max) break;
364 if (rsl_qsweep[nsweep] == 0) continue;
366 if (radar_verbose_flag)
367 fprintf(stderr,"Processing for SWEEP # %d\n", nsweep);
369 /* wsr88d_print_sweep_info(&wsr88d_sweep); */
371 for (iv=0; iv<nvolumes; iv++) {
372 if (rsl_qfield[iv]) {
373 /* Exceeded sweep limit.
374 * Allocate more sweeps.
375 * Copy all previous sweeps.
377 if (nsweep >= radar->v[iv]->h.nsweeps) {
378 if (radar_verbose_flag)
379 fprintf(stderr,"Exceeded sweep allocation of %d. "
380 "Adding 20 more.\n", nsweep);
381 new_volume = RSL_new_volume(radar->v[iv]->h.nsweeps+20);
382 new_volume = copy_sweeps_into_volume(new_volume, radar->v[iv]);
383 radar->v[iv] = new_volume;
385 if (wsr88d_load_sweep_into_volume(wsr88d_sweep,
386 radar->v[iv], nsweep, volume_mask[iv]) != 0) {
387 RSL_free_radar(radar);
393 /* Get Volume Coverage Pattern number for radar header. */
395 while (i < MAX_RAYS_IN_SWEEP && wsr88d_sweep.ray[i] == NULL) i++;
396 if (i < MAX_RAYS_IN_SWEEP) radar->h.vcp = wsr88d_get_volume_coverage(
397 wsr88d_sweep.ray[i]);
400 free_and_clear_sweep(&wsr88d_sweep, 0, MAX_RAYS_IN_SWEEP);
403 for (iv=0; iv<nvolumes; iv++) {
404 if (rsl_qfield[iv]) {
405 radar->v[iv]->h.type_str = strdup(field_str[iv]);
406 radar->v[iv]->h.nsweeps = nsweep;
413 * Here we will assign the Radar_header information. Take most of it
414 * from an existing volume's header.
416 radar_load_date_time(radar); /* Magic :-) */
418 radar->h.number = sitep->number;
419 memcpy(&radar->h.name, sitep->name, sizeof(sitep->name));
420 memcpy(&radar->h.radar_name, sitep->name, sizeof(sitep->name)); /* Redundant */
421 memcpy(&radar->h.city, sitep->city, sizeof(sitep->city));
422 memcpy(&radar->h.state, sitep->state, sizeof(sitep->state));
423 strcpy(radar->h.radar_type, "wsr88d");
424 radar->h.latd = sitep->latd;
425 radar->h.latm = sitep->latm;
426 radar->h.lats = sitep->lats;
427 if (radar->h.latd < 0) { /* Degree/min/sec all the same sign */
431 radar->h.lond = sitep->lond;
432 radar->h.lonm = sitep->lonm;
433 radar->h.lons = sitep->lons;
434 if (radar->h.lond < 0) { /* Degree/min/sec all the same sign */
438 radar->h.height = sitep->height;
439 radar->h.spulse = sitep->spulse;
440 radar->h.lpulse = sitep->lpulse;
444 radar = RSL_prune_radar(radar);