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.
29 extern int radar_verbose_flag;
30 /* Missing data flag : -32768 when a signed short. */
31 #define UF_NO_DATA 0X8000
34 /* Field names. Any convensions may be observed. */
36 * DZ = Reflectivity (dBZ).
37 * VR = Radial Velocity.
38 * SW = Spectrum Width.
39 * CZ = Corrected Reflectivity. (Quality controlled: AP removed, etc.)
40 * ZT = Total Reflectivity (dB(mW)). Becomes UZ in UF files.
41 * DR = Differential Reflectivity.
42 * LR = Another DR (LDR).
43 * ZD = Tina Johnson use this one.
44 * DM = Received power.
45 * RH = Rho coefficient.
46 * PH = Phi (MCTEX parameter).
47 * XZ = X-band reflectivity.
49 * MZ = DZ mask for 1C-51 HDF.
50 * MD = ZD mask for 1C-51 HDF.
51 * ZE = Edited reflectivity.
52 * VE = Edited velocity.
53 * KD = KDP wavelength*deg/km
54 * TI = TIME (units unknown).
55 * These fields may appear in any order in the UF file.
56 * There are more fields than appear here. See rsl.h.
60 /* Changed old buffer size (16384) for larger dualpol files. BLK 5/20/2011 */
61 typedef short UF_buffer[20000]; /* Bigger than documented 4096. */
63 void swap_uf_buffer(UF_buffer uf);
64 void swap2(short *buf, int n);
66 /**********************************************************************/
68 /* RSL_radar_to_uf_fp */
70 /* By: John Merritt */
71 /* Space Applications Corporation */
73 /**********************************************************************/
74 void RSL_radar_to_uf_fp(Radar *r, FILE *fp)
78 * 1. Fill the UF buffers with data from the Radar structure.
79 * 2. Write to a stream. Assume open and leave it so.
85 /* These are pointers to various locations within the UF buffer 'uf'.
86 * They are used to index the different components of the UF structure in
87 * a manor consistant with the UF documentation. For instance, uf_ma[1]
88 * will be equivalenced to the second word (2 bytes/each) of the UF
91 short *uf_ma; /* Mandatory header block. */
92 short *uf_op; /* Optional header block. */
93 short *uf_lu; /* Local Use header block. */
94 short *uf_dh; /* Data header. */
95 short *uf_fh; /* Field header. */
96 short *uf_data; /* Data. */
98 /* The length of each header. */
99 int len_ma, len_op, len_lu, len_dh, len_fh, len_data;
101 /* Booleans to flag inclusion of headers. */
102 int q_op, q_lu, q_dh, q_fh;
104 int current_fh_index;
106 int rec_len, save_rec_len;
118 int uf_sweep_mode = 1; /* default PPI */
120 /* Here are the arrays for each field type. Each dimension is the number
121 * of fields in the radar structure. I do this because the radar organization
122 * is by volumes (field types) and the UF demands that each ray contain
123 * all the field types.
129 int nvolumes, maxsweeps, nrays;
131 int sweep_num, ray_num, rec_num;
135 fprintf(stderr, "radar_to_uf_fp: radar pointer NULL\n");
139 /* Do all the headers first time around. Then, prune OP and LU. */
140 q_op = q_lu = q_dh = q_fh = 1;
142 memset(&uf, 0, sizeof(uf)); /* Init to 0 or NULL for pointers. */
144 sweep_num = ray_num = rec_num = 0;
145 true_nvolumes = nvolumes = maxsweeps = nrays = 0;
147 /* PPI and RHI are enum constants defined in rsl.h */
148 if (r->h.scan_mode == PPI) uf_sweep_mode = 1;
149 else if (r->h.scan_mode == RHI) uf_sweep_mode = 3;
152 * The organization of the Radar structure is by volumes, then sweeps, then
153 * rays, then gates. This is different from the UF file organization.
154 * The UF format wants sweeps, rays, then gates for all field types (volumes).
155 * So, we have to do a back flip, here. This is achieved by maintaining
156 * an array of volume pointers and sweep pointers, each dimensioned by
157 * 'nvolumes', which contains the data for the different field types; this
158 * is our innermost loop. The variables are 'volume[i]' and 'sweep[i]' where
159 * 'i' is the volume index.
161 * In other words, we are getting all the field types together, when we
162 * are looping on the number of rays in a sweep, so we can load the UF_buffer
166 nvolumes = r->h.nvolumes;
167 volume = (Volume **) calloc(nvolumes, sizeof(Volume *));
168 sweep = (Sweep **) calloc(nvolumes, sizeof(Sweep *));
169 nsweeps = (int *) calloc(nvolumes, sizeof(int));
171 /* Get the the number of sweeps in the radar structure. This will be
172 * the main controlling loop variable.
174 for (i=0; i<nvolumes; i++) {
177 nsweeps[i] = volume[i]->h.nsweeps;
178 if (nsweeps[i] > maxsweeps) maxsweeps = nsweeps[i];
183 if (radar_verbose_flag) {
184 fprintf(stderr,"True number of volumes for UF is %d\n", true_nvolumes);
185 fprintf(stderr,"Maximum # of volumes for UF is %d\n", nvolumes);
188 max_field_names = MAX_RADAR_VOLUMES;
191 * LOOP for all sweeps (typically 11 or 16 for wsr88d data.
194 for (i=0; i<maxsweeps; i++) {
195 /* Get the array of volume and sweep pointers; one for each field type. */
197 for (k=0; k<nvolumes; k++) {
198 if (volume[k]) sweep[k] = volume[k]->sweep[i];
200 /* Check if we really can access this sweep. Paul discovered that
201 * if the actual number of sweeps is less than the maximum that we
202 * could be chasing a bad pointer (a NON-NULL garbage pointer).
204 if (i >= nsweeps[k]) sweep[k] = NULL;
206 if (sweep[k]) if (sweep[k]->h.nrays > nrays) nrays = sweep[k]->h.nrays;
209 sweep_num++; /* I guess it will be ok to count NULL sweeps. */
211 if (radar_verbose_flag)
212 fprintf(stderr,"Processing sweep %d for %d rays.", i, nrays);
213 if (radar_verbose_flag) {
214 if (little_endian()) fprintf(stderr," ... On Little endian.\n");
215 else fprintf(stderr,"\n");
219 /* Now LOOP for all rays within this particular sweep (i).
220 * Get all the field types together for the ray, see ray[k], and
221 * fill the UF data buffer appropriately.
223 for (j=0; j<nrays; j++) {
224 memset(uf, 0, sizeof(uf));
226 ray_num++; /* And counting, possibly, NULL rays. */
227 current_fh_index = 0;
230 /* Find any ray for header information. It does not matter which
231 * ray, since the information for the MANDITORY, OPTIONAL, and LOCAL
232 * USE headers is common to any field type ray.
235 for (k=0; k<nvolumes; k++) {
237 if (j < sweep[k]->h.nrays)
239 if ((ray = sweep[k]->ray[j])) break;
242 /* If there is no such ray, then continue on to the next ray. */
245 fprintf(stderr,"Ray: %.4d, Time: %2.2d:%2.2d:%f %.2d/%.2d/%.4d\n", ray_num, ray->h.hour, ray->h.minute, ray->h.sec, ray->h.month, ray->h.day, ray->h.year);
249 * ---- Begining of MANDITORY HEADER BLOCK.
252 memcpy(&uf_ma[0], "UF", 2);
253 if (little_endian()) memcpy(&uf_ma[0], "FU", 2);
254 uf_ma[1] = 0; /* Not known yet. */
255 uf_ma[2] = 0; /* Not known yet. Really, I do. */
256 uf_ma[3] = 0; /* Not known yet. */
257 uf_ma[4] = 0; /* Not known yet. */
262 uf_ma[9 ] = sweep_num;
263 memcpy(&uf_ma[10], r->h.radar_name, 8);
264 if (little_endian()) swap2(&uf_ma[10], 8/2);
265 memcpy(&uf_ma[14], r->h.name, 8);
266 if (little_endian()) swap2(&uf_ma[14], 8/2);
267 /* Convert decimal lat/lon to d:m:s */
269 if (ray->h.lat != 0.0) {
270 degree = (int)ray->h.lat;
271 minute = (int)((ray->h.lat - degree) * 60);
272 second = (ray->h.lat - degree - minute/60.0) * 3600.0;
280 if (second > 0.0) uf_ma[20] = second*64 + 0.5;
281 else uf_ma[20] = second*64 - 0.5;
283 if (ray->h.lon != 0.0) {
284 degree = (int)ray->h.lon;
285 minute = (int)((ray->h.lon - degree) * 60);
286 second = (ray->h.lon - degree - minute/60.0) * 3600.0;
294 if (second > 0.0) uf_ma[23] = second*64 + 0.5;
295 else uf_ma[23] = second*64 - 0.5;
297 uf_ma[24] = ray->h.alt;
299 uf_ma[24] = r->h.height;
301 uf_ma[25] = ray->h.year % 100; /* By definition: not year 2000 compliant. */
302 uf_ma[26] = ray->h.month;
303 uf_ma[27] = ray->h.day;
304 uf_ma[28] = ray->h.hour;
305 uf_ma[29] = ray->h.minute;
306 uf_ma[30] = ray->h.sec;
307 memcpy(&uf_ma[31], "UT", 2);
308 if (little_endian()) memcpy(&uf_ma[31], "TU", 2);
309 if (ray->h.azimuth > 0) uf_ma[32] = ray->h.azimuth*64 + 0.5;
310 else uf_ma[32] = ray->h.azimuth*64 - 0.5;
311 uf_ma[33] = ray->h.elev*64 + 0.5;
312 uf_ma[34] = uf_sweep_mode;
313 if (ray->h.fix_angle != 0.)
314 uf_ma[35] = ray->h.fix_angle*64.0 + 0.5;
315 else uf_ma[35] = sweep[k]->h.elev*64.0 + 0.5;
316 uf_ma[36] = ray->h.sweep_rate*(360.0/60.0)*64.0 + 0.5;
318 the_time = time(NULL);
319 tm = gmtime(&the_time);
321 uf_ma[37] = tm->tm_year % 100; /* Same format as data year */
322 uf_ma[38] = tm->tm_mon+1;
323 uf_ma[39] = tm->tm_mday;
324 memcpy(&uf_ma[40], "RSL" VERSION, 8);
325 if (little_endian()) swap2(&uf_ma[40], 8/2);
326 uf_ma[44] = (signed short)UF_NO_DATA;
330 * ---- End of MANDITORY HEADER BLOCK.
333 /* ---- Begining of OPTIONAL HEADER BLOCK. */
336 q_op = 0; /* Only once. */
338 memcpy(&uf_op[0], "TRMMGVUF", 8);
339 if (little_endian()) swap2(&uf_op[0], 8/2);
340 uf_op[4] = (signed short)UF_NO_DATA;
341 uf_op[5] = (signed short)UF_NO_DATA;
342 uf_op[6] = ray->h.hour;
343 uf_op[7] = ray->h.minute;
344 uf_op[8] = ray->h.sec;
345 memcpy(&uf_op[9], "RADAR_UF", 8);
346 if (little_endian()) swap2(&uf_op[9], 8/2);
350 /* ---- End of OPTIONAL HEADER BLOCK. */
352 /* ---- Begining of LOCAL USE HEADER BLOCK. */
355 /* TODO: Code within "#ifdef LUHDR_VR_AZ" below should be removed
356 * once testing of merge_split_cuts is completed.
358 /* 5/18/2010 Temporarily define LUHDR_VR_AZ until merge_split_cuts is
362 /* If DZ and VR azimuths are different, store VR azimuth in Local Use
363 * Header. This is done for WSR-88D split cuts.
365 if (sweep[DZ_INDEX] && sweep[VR_INDEX]) {
366 if (sweep[DZ_INDEX]->ray[j] && sweep[VR_INDEX]->ray[j]) {
367 vr_az = sweep[VR_INDEX]->ray[j]->h.azimuth;
368 if (sweep[DZ_INDEX]->ray[j]->h.azimuth != vr_az)
369 q_lu = 1; /* Set to use Local Use Header block. */
375 /* Store azimuth for WSR-88D VR ray in Local Use Header. */
376 uf_lu = uf+len_ma+len_op;
377 memcpy(&uf_lu[0], "AZ", 2);
378 if (little_endian()) memcpy(&uf_lu[0], "ZA", 2);
379 if (vr_az > 0) uf_lu[1] = vr_az*64 + 0.5;
380 else uf_lu[1] = vr_az*64 - 0.5;
383 /* ---- End of LOCAL USE HEADER BLOCK. */
386 /* Here is where we loop on each field type. We need to keep
387 * track of how many FIELD HEADER and FIELD DATA sections, one
388 * for each field type, we fill. The variable that tracks this
389 * index into 'uf' is 'current_fh_index'. It is bumped by
390 * the length of the FIELD HEADER and FIELD DATA for each field
391 * type encountered. Field types expected are: Reflectivity,
392 * Velocity, and Spectrum width; this is a typicial list but it
393 * is not restricted to it.
396 for (k=0; k<nvolumes; k++) {
398 if (j < sweep[k]->h.nrays && sweep[k]->ray[j])
399 ray = sweep[k]->ray[j];
405 /* ---- Begining of DATA HEADER. */
408 len_dh = 2*true_nvolumes + 3;
409 uf_dh = uf+len_ma+len_op+len_lu;
413 /* 'nfield' indexes the field number.
414 * 'k' indexes the particular field from the volume.
415 * RSL_ftype contains field names and is defined in rsl.h.
417 if (k > max_field_names-1) {
419 "RSL_uf_to_radar: No field name for volume index %d\n", k);
420 fprintf(stderr,"RSL_ftype must be updated in rsl.h for new field.\n");
421 fprintf(stderr,"Quitting now.\n");
424 memcpy(&uf_dh[3+2*(nfield-1)], RSL_ftype[k], 2);
425 if (little_endian()) swap2(&uf_dh[3+2*(nfield-1)], 2/2);
426 if (current_fh_index == 0) current_fh_index = len_ma+len_op+len_lu+len_dh;
427 uf_dh[4+2*(nfield-1)] = current_fh_index + 1;
429 /* ---- End of DATA HEADER. */
431 /* ---- Begining of FIELD HEADER. */
433 uf_fh = uf+current_fh_index;
434 if (k != PH_INDEX) scale_factor = 100;
435 else scale_factor = 10;
436 uf_fh[1] = scale_factor;
437 uf_fh[2] = ray->h.range_bin1/1000.0;
438 uf_fh[3] = ray->h.range_bin1 - (1000*uf_fh[2]);
439 uf_fh[4] = ray->h.gate_size;
440 uf_fh[5] = ray->h.nbins;
441 uf_fh[6] = ray->h.pulse_width*(RSL_SPEED_OF_LIGHT/1.0e6);
442 uf_fh[7] = sweep[k]->h.beam_width*64.0 + 0.5;
443 uf_fh[8] = sweep[k]->h.beam_width*64.0 + 0.5;
444 uf_fh[9] = ray->h.frequency*64.0 + 0.5; /* Bandwidth (mHz). */
445 uf_fh[10] = 0; /* Horizontal polarization. */
446 uf_fh[11] = ray->h.wavelength*64.0*100.0; /* m to cm. */
447 uf_fh[12] = ray->h.pulse_count;
448 memcpy(&uf_fh[13], " ", 2);
449 uf_fh[14] = (signed short)UF_NO_DATA;
450 uf_fh[15] = (signed short)UF_NO_DATA;
451 if (k == DZ_INDEX || k == ZT_INDEX) {
452 uf_fh[16] = volume[k]->h.calibr_const*100.0 + 0.5;
455 memcpy(&uf_fh[16], " ", 2);
458 uf_fh[17] = 1.0/ray->h.prf*1000000.0; /* Pulse repetition time(msec) = 1/prf */
460 uf_fh[17] = (signed short)UF_NO_DATA; /* Pulse repetition time = 1/prf */
462 if (VR_INDEX == k || VE_INDEX == k) {
463 uf_fh[19] = scale_factor*ray->h.nyq_vel;
470 uf_fh[0] = current_fh_index + len_fh + 1;
471 /* ---- End of FIELD HEADER. */
473 /* ---- Begining of FIELD DATA. */
474 uf_data = uf+len_fh+current_fh_index;
475 len_data = ray->h.nbins;
476 for (m=0; m<len_data; m++) {
477 x = ray->h.f(ray->range[m]);
478 if (x == BADVAL || x == RFVAL || x == APFLAG || x == NOECHO)
479 uf_data[m] = (signed short)UF_NO_DATA;
481 uf_data[m] = scale_factor * x;
484 current_fh_index += (len_fh+len_data);
487 /* ---- End of FIELD DATA. */
489 /* Fill in some infomation we didn't know. Like, buffer length,
490 * record number, etc.
493 uf_ma[1] = current_fh_index;
494 uf_ma[3] = len_ma + len_op + 1;
495 uf_ma[4] = len_ma + len_op + len_lu + 1;
498 /* WRITE the UF buffer. */
499 rec_len =(int)uf_ma[1]*2;
500 save_rec_len = rec_len; /* We destroy 'rec_len' when making it
501 big endian on a little endian machine. */
503 swap_4_bytes(&rec_len);
504 if (fwrite(&rec_len, sizeof(int), 1, fp) != 1)
505 perror("RSL_radar_to_uf_fp: short write");
508 if (fwrite(uf, sizeof(char), save_rec_len, fp) != 1)
509 perror("RSL_radar_to_uf_fp: short write");
510 if (fwrite(&rec_len, sizeof(int), 1, fp) != 1)
511 perror("RSL_radar_to_uf_fp: short write");
517 /**********************************************************************/
519 /* RSL_radar_to_uf */
521 /**********************************************************************/
522 void RSL_radar_to_uf(Radar *r, char *outfile)
526 fprintf(stderr, "radar_to_uf: radar pointer NULL\n");
530 if ((fp = fopen(outfile, "wb")) == NULL) {
535 RSL_radar_to_uf_fp(r, fp);
539 /**********************************************************************/
541 /* RSL_radar_to_uf_gzip */
543 /**********************************************************************/
544 void RSL_radar_to_uf_gzip(Radar *r, char *outfile)
548 fprintf(stderr, "radar_to_uf_gzip: radar pointer NULL\n");
552 if ((fp = fopen(outfile, "wb")) == NULL) {
557 fp = compress_pipe(fp);
558 RSL_radar_to_uf_fp(r, fp);