3 This is the TRMM Office Radar Software Library.
9 This program is free software: you can redistribute it and/or modify
10 it under the terms of the GNU General Public License as published by
11 the Free Software Foundation, either version 3 of the License, or
12 (at your option) any later version.
14 This program 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
17 GNU General Public License for more details.
19 You should have received a copy of the GNU General Public License
20 along with this program. If not, see <http://www.gnu.org/licenses/>.
25 * This file contains routines for processing Message Type 31, the digital
26 * radar message type introduced in WSR-88D Level II Build 10.
33 /* Data descriptions in the following data structures are from the "Interface
34 * Control Document for the RDA/RPG", Build 10.0 Draft, WSR-88D Radar
39 short rpg[6]; /* 12 bytes inserted by RPG Communications Mgr. Ignored. */
40 unsigned short msg_size; /* Message size for this segment, in halfwords */
41 unsigned char channel; /* RDA Redundant Channel */
42 unsigned char msg_type; /* Message type. For example, 31 */
43 unsigned short id_seq; /* Msg seq num = 0 to 7FFF, then roll over to 0 */
44 unsigned short msg_date; /* Modified Julian date from 1/1/70 */
45 unsigned int msg_time; /* Packet generation time in ms past midnight */
46 unsigned short num_segs; /* Number of segments for this message */
47 unsigned short seg_num; /* Number of this segment */
52 unsigned int ray_time; /* Data collection time in milliseconds past midnight GMT */
53 unsigned short ray_date; /* Julian date - 2440586.5 (1/01/1970) */
54 unsigned short azm_num ; /* Radial number within elevation scan */
55 float azm; /* Azimuth angle in degrees (0 to 359.956055) */
56 unsigned char compression_code; /* 0 = uncompressed, 1 = BZIP2, 2 = zlib */
57 unsigned char spare; /* for word alignment */
58 unsigned short radial_len; /* radial length in bytes, including data header block */
59 unsigned char azm_res;
60 unsigned char radial_status;
61 unsigned char elev_num;
62 unsigned char sector_cut_num;
63 float elev; /* Elevation angle in degrees (-7.0 to 70.0) */
64 unsigned char radial_spot_blanking;
65 unsigned char azm_indexing_mode;
66 unsigned short data_block_count;
67 /* Data Block Pointers */
68 unsigned int dbptr_vol_const;
69 unsigned int dbptr_elev_const;
70 unsigned int dbptr_radial_const;
71 unsigned int dbptr_ref;
72 unsigned int dbptr_vel;
73 unsigned int dbptr_sw;
74 unsigned int dbptr_zdr;
75 unsigned int dbptr_phi;
76 unsigned int dbptr_rho;
77 } Ray_header_m31; /* Called Data Header Block in RDA/RPG document. */
81 unsigned int reserved;
82 unsigned short ngates;
83 short range_first_gate;
84 short range_samp_interval;
85 short thresh_not_overlayed;
87 unsigned char controlflag;
88 unsigned char datasize_bits;
93 #define MAX_RADIAL_LENGTH 14288
96 Ray_header_m31 ray_hdr;
99 unsigned char data[MAX_RADIAL_LENGTH];
103 enum radial_status {START_OF_ELEV, INTERMED_RADIAL, END_OF_ELEV, BEGIN_VOS,
107 void wsr88d_swap_m31_hdr(Wsr88d_msg_hdr *msghdr)
109 swap_2_bytes(&msghdr->msg_size);
110 swap_2_bytes(&msghdr->id_seq);
111 swap_2_bytes(&msghdr->msg_date);
112 swap_4_bytes(&msghdr->msg_time);
113 swap_2_bytes(&msghdr->num_segs);
114 swap_2_bytes(&msghdr->seg_num);
118 void wsr88d_swap_m31_ray(Ray_header_m31 *wsr88d_ray)
122 swap_4_bytes(&wsr88d_ray->ray_time);
123 swap_2_bytes(&wsr88d_ray->ray_date);
124 swap_2_bytes(&wsr88d_ray->azm_num);
125 swap_4_bytes(&wsr88d_ray->azm);
126 swap_2_bytes(&wsr88d_ray->radial_len);
127 swap_4_bytes(&wsr88d_ray->elev);
128 swap_2_bytes(&wsr88d_ray->data_block_count);
129 data_ptr = (int *) &wsr88d_ray->dbptr_vol_const;
130 for (; data_ptr <= (int *) &wsr88d_ray->dbptr_rho; data_ptr++)
131 swap_4_bytes(data_ptr);
135 void wsr88d_swap_data_hdr(Data_moment_hdr *this_field)
138 halfword = (short *) &this_field->ngates;
139 for (; halfword < (short *) &this_field->controlflag; halfword++)
140 swap_2_bytes(halfword);
141 swap_4_bytes(&this_field->scale);
142 swap_4_bytes(&this_field->offset);
146 float wsr88d_get_angle(short bitfield)
151 float value[13] = {0.043945, 0.08789, 0.17578, 0.35156, .70313, 1.40625,
152 2.8125, 5.625, 11.25, 22.5, 45., 90., 180.};
154 /* find which bits are set and sum corresponding values to get angle. */
156 bitfield = bitfield >> 3; /* 3 least significant bits aren't used. */
157 for (i = 0; i < 13; i++) {
158 if (bitfield & mask) angle += value[i];
159 bitfield = bitfield >> 1;
165 float wsr88d_get_azim_rate(short bitfield)
170 float value[12] = {0.0109863, 0.021972656, 0.043945, 0.08789, 0.17578,
171 0.35156, .70313, 1.40625, 2.8125, 5.625, 11.25, 22.5};
173 /* Find which bits are set and sum corresponding values to get rate. */
175 bitfield = bitfield >> 3; /* 3 least significant bits aren't used. */
176 for (i = 0; i < 12; i++) {
177 if (bitfield & mask) rate += value[i];
178 bitfield = bitfield >> 1;
180 if (bitfield >> 15) rate = -rate;
184 #define WSR88D_MAX_SWEEPS 20
190 float fixed_angle[WSR88D_MAX_SWEEPS];
191 float azim_rate[WSR88D_MAX_SWEEPS];
192 int waveform[WSR88D_MAX_SWEEPS];
193 int super_res_ctrl[WSR88D_MAX_SWEEPS];
194 int surveil_prf_num[WSR88D_MAX_SWEEPS];
195 int doppler_prf_num[WSR88D_MAX_SWEEPS];
198 static VCP_data vcp_data;
200 void wsr88d_get_vcp_data(short *msgtype5)
202 short azim_rate, fixed_angle, vel_res;
203 short sres_and_survprf; /* super res ctrl and surveil prf, one byte each */
204 short chconf_and_waveform;
207 vcp_data.vcp = (unsigned short) msgtype5[2];
208 vcp_data.num_cuts = msgtype5[3];
209 if (little_endian()) {
210 swap_2_bytes(&vcp_data.vcp);
211 swap_2_bytes(&vcp_data.num_cuts);
213 vel_res = msgtype5[5];
214 if (little_endian()) swap_2_bytes(&vel_res);
215 vel_res = vel_res >> 8;
216 if (vel_res == 2) vcp_data.vel_res = 0.5;
217 else if (vel_res == 4) vcp_data.vel_res = 1.0;
218 else vcp_data.vel_res = 0.0;
219 /* Get elevation related information for each sweep. */
220 for (i=0; i < vcp_data.num_cuts; i++) {
221 fixed_angle = msgtype5[11 + i*23];
222 azim_rate = msgtype5[15 + i*23];
223 chconf_and_waveform = msgtype5[12 + i*23];
224 sres_and_survprf = msgtype5[13 + i*23];
225 vcp_data.doppler_prf_num[i] = msgtype5[23 + i*23];
226 if (little_endian()) {
227 swap_2_bytes(&fixed_angle);
228 swap_2_bytes(&azim_rate);
229 swap_2_bytes(&chconf_and_waveform);
230 swap_2_bytes(&sres_and_survprf);
231 swap_2_bytes(&vcp_data.doppler_prf_num[i]);
233 vcp_data.fixed_angle[i] = wsr88d_get_angle(fixed_angle);
234 vcp_data.azim_rate[i] = wsr88d_get_azim_rate(azim_rate);
235 vcp_data.waveform[i] = chconf_and_waveform & 0xff;
236 vcp_data.super_res_ctrl[i] = sres_and_survprf >> 8;
237 vcp_data.surveil_prf_num[i] = sres_and_survprf & 0xff;
242 void get_wsr88d_unamb_and_nyq_vel(Wsr88d_ray_m31 *wsr88d_ray, float *unamb_rng,
245 int dbptr, found, ray_hdr_len;
246 short nyq_vel_sh, unamb_rng_sh;
249 ray_hdr_len = sizeof(wsr88d_ray->ray_hdr);
250 dbptr = wsr88d_ray->ray_hdr.dbptr_radial_const - ray_hdr_len;
251 if (strncmp((char *) &wsr88d_ray->data[dbptr], "RRAD", 4) == 0) found = 1;
253 dbptr = wsr88d_ray->ray_hdr.dbptr_elev_const - ray_hdr_len;
254 if (strncmp((char *) &wsr88d_ray->data[dbptr], "RRAD", 4) == 0)
257 dbptr = wsr88d_ray->ray_hdr.dbptr_vol_const - ray_hdr_len;
258 if (strncmp((char *) &wsr88d_ray->data[dbptr], "RRAD", 4) == 0)
263 memcpy(&unamb_rng_sh, &wsr88d_ray->data[dbptr+6], 2);
264 memcpy(&nyq_vel_sh, &wsr88d_ray->data[dbptr+16], 2);
265 if (little_endian()) {
266 swap_2_bytes(&unamb_rng_sh);
267 swap_2_bytes(&nyq_vel_sh);
269 *unamb_rng = unamb_rng_sh / 10.;
270 *nyq_vel = nyq_vel_sh / 100.;
278 int read_wsr88d_ray_m31(Wsr88d_file *wf, int msg_size,
279 Wsr88d_ray_m31 *wsr88d_ray)
281 Ray_header_m31 ray_hdr;
282 int n, bytes_to_read, ray_hdr_len;
283 float nyq_vel, unamb_rng;
285 /* Read ray data header block */
286 n = fread(&wsr88d_ray->ray_hdr, sizeof(Ray_header_m31), 1, wf->fptr);
288 fprintf(stderr,"read_wsr88d_ray_m31: Read failed.\n");
292 /* Byte swap header if needed. */
293 if (little_endian()) wsr88d_swap_m31_ray(&wsr88d_ray->ray_hdr);
295 ray_hdr = wsr88d_ray->ray_hdr;
296 ray_hdr_len = sizeof(ray_hdr);
298 /* Read data portion of radial. */
300 bytes_to_read = msg_size - ray_hdr_len;
301 n = fread(wsr88d_ray->data, bytes_to_read, 1, wf->fptr);
303 fprintf(stderr,"read_wsr88d_ray_m31: Read failed.\n");
307 /* We retrieve unambiguous range and Nyquist velocity here so that we don't
308 * have to do it repeatedly for each data moment later.
310 get_wsr88d_unamb_and_nyq_vel(wsr88d_ray, &unamb_rng, &nyq_vel);
311 wsr88d_ray->unamb_rng = unamb_rng;
312 wsr88d_ray->nyq_vel = nyq_vel;
318 void wsr88d_load_ray_hdr(Wsr88d_ray_m31 wsr88d_ray, Ray *ray)
320 int month, day, year, hour, minute, sec;
323 Ray_header_m31 ray_hdr;
325 ray_hdr = wsr88d_ray.ray_hdr;
326 m1_ray.ray_date = ray_hdr.ray_date;
327 m1_ray.ray_time = ray_hdr.ray_time;
329 wsr88d_get_date(&m1_ray, &month, &day, &year);
330 wsr88d_get_time(&m1_ray, &hour, &minute, &sec, &fsec);
331 ray->h.year = year + 1900;
332 ray->h.month = month;
335 ray->h.minute = minute;
336 ray->h.sec = sec + fsec;
337 ray->h.azimuth = ray_hdr.azm;
338 ray->h.ray_num = ray_hdr.azm_num;
339 ray->h.elev = ray_hdr.elev;
340 ray->h.elev_num = ray_hdr.elev_num;
341 ray->h.unam_rng = wsr88d_ray.unamb_rng;
342 ray->h.nyq_vel = wsr88d_ray.nyq_vel;
344 elev_index = ray_hdr.elev_num - 1;
345 ray->h.azim_rate = vcp_data.azim_rate[elev_index];
346 ray->h.fix_angle = vcp_data.fixed_angle[elev_index];
347 ray->h.vel_res = vcp_data.vel_res;
348 if (ray_hdr.azm_res != 1)
349 ray->h.beam_width = 1.0;
350 else ray->h.beam_width = 0.5;
352 /* For convenience, use message type 1 routines to get some values.
353 * First load VCP and elevation numbers into a msg 1 ray.
355 m1_ray.vol_cpat = vcp_data.vcp;
356 m1_ray.elev_num = ray_hdr.elev_num;
357 m1_ray.unam_rng = (short) (wsr88d_ray.unamb_rng * 10.);
358 /* Get values from message type 1 routines. */
359 ray->h.frequency = wsr88d_get_frequency(&m1_ray);
360 ray->h.pulse_width = wsr88d_get_pulse_width(&m1_ray);
361 ray->h.pulse_count = wsr88d_get_pulse_count(&m1_ray);
362 ray->h.prf = (int) wsr88d_get_prf(&m1_ray);
363 ray->h.wavelength = 0.1071;
367 int wsr88d_get_vol_index(char* dataname)
371 if (strncmp(dataname, "DREF", 4) == 0) vol_index = DZ_INDEX;
372 if (strncmp(dataname, "DVEL", 4) == 0) vol_index = VR_INDEX;
373 if (strncmp(dataname, "DSW", 3) == 0) vol_index = SW_INDEX;
374 if (strncmp(dataname, "DZDR", 3) == 0) vol_index = DR_INDEX;
375 if (strncmp(dataname, "DPHI", 3) == 0) vol_index = PH_INDEX;
376 if (strncmp(dataname, "DRHO", 3) == 0) vol_index = RH_INDEX;
382 #define MAXRAYS_M31 800
385 void wsr88d_load_ray(Wsr88d_ray_m31 wsr88d_ray, int data_ptr,
386 int isweep, int iray, Radar *radar)
388 /* Load data into ray structure for this field or data moment. */
390 Data_moment_hdr data_hdr;
393 float value, scale, offset;
395 Range (*invf)(float x);
398 int vol_index, waveform;
400 enum waveforms {surveillance=1, doppler_w_amb_res, doppler_no_amb_res,
403 /* Get data moment header. */
404 hdr_size = sizeof(data_hdr);
405 memcpy(&data_hdr, &wsr88d_ray.data[data_ptr], hdr_size);
406 data_ptr += hdr_size;
407 if (little_endian()) wsr88d_swap_data_hdr(&data_hdr);
409 vol_index = wsr88d_get_vol_index(data_hdr.dataname);
411 fprintf(stderr,"wsr88d_load_ray: Unknown dataname %s. isweep = %d, "
412 "iray = %d.\n", data_hdr.dataname, isweep, iray);
416 case DZ_INDEX: f = DZ_F; invf = DZ_INVF; break;
417 case VR_INDEX: f = VR_F; invf = VR_INVF; break;
418 case SW_INDEX: f = SW_F; invf = SW_INVF; break;
419 case DR_INDEX: f = DR_F; invf = DR_INVF; break;
420 case PH_INDEX: f = PH_F; invf = PH_INVF; break;
421 case RH_INDEX: f = RH_F; invf = RH_INVF; break;
424 waveform = vcp_data.waveform[isweep];
425 /* The following somewhat complicated if-statement says we'll do the
426 * body of the statement if:
427 * a) the data moment is not reflectivity, or
428 * b) the data moment is reflectivity and one of the following is true:
429 * - waveform is surveillance
430 * - waveform is batch
431 * - elevation is greater than highest split cut (i.e., 6 deg.)
433 if (vol_index != DZ_INDEX || (waveform == surveillance ||
434 waveform == batch || vcp_data.fixed_angle[isweep] >= 6.0)) {
435 if (radar->v[vol_index] == NULL) {
436 radar->v[vol_index] = RSL_new_volume(MAXSWEEPS);
437 radar->v[vol_index]->h.f = f;
438 radar->v[vol_index]->h.invf = invf;
440 if (radar->v[vol_index]->sweep[isweep] == NULL) {
441 radar->v[vol_index]->sweep[isweep] = RSL_new_sweep(MAXRAYS_M31);
442 radar->v[vol_index]->sweep[isweep]->h.f = f;
443 radar->v[vol_index]->sweep[isweep]->h.invf = invf;
445 ngates = data_hdr.ngates;
446 ray = RSL_new_ray(ngates);
448 /* Convert data to float, then use range function to store in ray.
449 * Note: data range is 2-255. 0 means signal is below threshold, and 1
450 * means range folded.
453 offset = data_hdr.offset;
454 scale = data_hdr.scale;
455 if (data_hdr.scale == 0) scale = 1.0;
456 data = &wsr88d_ray.data[data_ptr];
457 for (i = 0; i < ngates; i++) {
459 value = (data[i] - offset) / scale;
460 else value = (data[i] == 0) ? BADVAL : RFVAL;
461 ray->range[i] = invf(value);
465 wsr88d_load_ray_hdr(wsr88d_ray, ray);
466 ray->h.range_bin1 = data_hdr.range_first_gate;
467 ray->h.gate_size = data_hdr.range_samp_interval;
468 ray->h.nbins = ngates;
469 radar->v[vol_index]->sweep[isweep]->ray[iray] = ray;
470 radar->v[vol_index]->sweep[isweep]->h.nrays = iray+1;
475 void wsr88d_load_ray_into_radar(Wsr88d_ray_m31 wsr88d_ray, int isweep, int iray,
478 int *data_ptr, hdr_size;
479 int i, ndatablocks, nconstblocks = 3;
481 hdr_size = sizeof(wsr88d_ray.ray_hdr);
483 ndatablocks = wsr88d_ray.ray_hdr.data_block_count;
484 data_ptr = (int *) &wsr88d_ray.ray_hdr.dbptr_ref;
485 for (i=0; i < ndatablocks-nconstblocks; i++) {
486 wsr88d_load_ray(wsr88d_ray, *data_ptr-hdr_size, isweep, iray, radar);
492 void wsr88d_load_sweep_header(Radar *radar, int isweep)
498 for (ivolume=0; ivolume < MAX_RADAR_VOLUMES; ivolume++) {
499 if (radar->v[ivolume] != NULL &&
500 radar->v[ivolume]->sweep[isweep] != NULL) {
501 sweep = radar->v[ivolume]->sweep[isweep];
502 nrays = sweep->h.nrays;
503 if (nrays == 0) continue;
504 last_ray = sweep->ray[nrays-1];
505 sweep->h.sweep_num = last_ray->h.elev_num;
506 sweep->h.elev = vcp_data.fixed_angle[isweep];
507 sweep->h.beam_width = last_ray->h.beam_width;
508 sweep->h.vert_half_bw = sweep->h.beam_width / 2.;
509 sweep->h.horz_half_bw = sweep->h.beam_width / 2.;
515 Radar *load_wsr88d_m31_into_radar(Wsr88d_file *wf)
517 Wsr88d_msg_hdr msghdr;
518 Wsr88d_ray_m31 wsr88d_ray;
519 short non31_seg_remainder[1202]; /* Remainder after message header */
520 int end_of_vos = 0, isweep = 0, iray = 0;
521 int msg_hdr_size, msg_size, n;
524 /* Message type 31 is a variable length message. All other message types
525 * are made up of 2432-byte segments. To handle all types, we read the
526 * message header and check the message type. If it is not 31, we simply
527 * read the remainder of the 2432-byte segment. If message type is 31, we
528 * use the size given in message header to determine how many bytes to
529 * read. For more information, see the "Interface Control Document for the
530 * RDA/RPG" at the WSR-88D Radar Operations Center web site.
533 n = fread(&msghdr, sizeof(Wsr88d_msg_hdr), 1, wf->fptr);
535 /* printf("msgtype = %d\n", msghdr.msg_type); */
536 msg_hdr_size = sizeof(Wsr88d_msg_hdr) - sizeof(msghdr.rpg);
538 radar = RSL_new_radar(MAX_RADAR_VOLUMES);
540 while (! end_of_vos) {
541 if (msghdr.msg_type == 31) {
542 if (little_endian()) wsr88d_swap_m31_hdr(&msghdr);
544 /* Get size of the remainder of message. The given size is in
545 * halfwords, but we want it in bytes, so double it.
547 msg_size = (int) msghdr.msg_size * 2 - msg_hdr_size;
549 n = read_wsr88d_ray_m31(wf, msg_size, &wsr88d_ray);
550 /* Assume error message was issued from read routine */
551 if (n <= 0) return NULL;
553 /* We need to check radial status for start of new elevation.
554 * Sometimes this occurs without an end-of-elevation flag for the
555 * previous sweep, which is the trigger for loading the sweep
556 * header. "iray" is set to zero after end-of-elev is received,
557 * so that's why we check it. We issue a warning because when this
558 * occurs, the number of rays in the previous sweep is usually less
561 if (wsr88d_ray.ray_hdr.radial_status == START_OF_ELEV &&
563 fprintf(stderr,"Warning: Radial status is Start-of-Elevation, "
564 "but End-of-Elevation was not\n"
565 "issued for elevation number %d. Number of rays = %d"
566 "\n", isweep+1, iray);
567 wsr88d_load_sweep_header(radar, isweep);
572 /* Load ray into radar structure. */
573 wsr88d_load_ray_into_radar(wsr88d_ray, isweep, iray, radar);
575 if (iray >= MAXRAYS_M31) {
576 fprintf(stderr,"Error: iray = %d, equals or exceeds MAXRAYS_M31"
577 " (%d)\n", iray, MAXRAYS_M31);
578 fprintf(stderr,"isweep = %d\n", isweep);
579 RSL_free_radar(radar);
583 else { /* msg_type not 31 */
584 n = fread(&non31_seg_remainder, sizeof(non31_seg_remainder), 1,
587 fprintf(stderr,"Warning: load_wsr88d_m31_into_radar: ");
588 if (feof(wf->fptr) != 0)
589 fprintf(stderr, "Unexpected end of file.\n");
591 fprintf(stderr,"Read failed.\n");
592 fprintf(stderr,"Current sweep index: %d\n"
593 "Last ray read: %d\n", isweep, iray);
594 wsr88d_load_sweep_header(radar, isweep);
597 if (msghdr.msg_type == 5) {
598 wsr88d_get_vcp_data(non31_seg_remainder);
599 radar->h.vcp = vcp_data.vcp;
600 /* printf("VCP = %d\n", vcp_data.vcp); */
604 /* Check for end of sweep */
605 if (wsr88d_ray.ray_hdr.radial_status == END_OF_ELEV) {
606 wsr88d_load_sweep_header(radar, isweep);
611 /* If not at end of volume scan, read next message header. */
612 if (wsr88d_ray.ray_hdr.radial_status != END_VOS) {
613 n = fread(&msghdr, sizeof(Wsr88d_msg_hdr), 1, wf->fptr);
615 fprintf(stderr,"Warning: load_wsr88d_m31_into_radar: ");
616 if (feof(wf->fptr) != 0) fprintf(stderr,
617 "Unexpected end of file.\n");
618 else fprintf(stderr,"Failed reading msghdr.\n");
619 fprintf(stderr,"Current sweep index: %d\n"
620 "Last ray read: %d\n", isweep, iray);
621 wsr88d_load_sweep_header(radar, isweep);
627 wsr88d_load_sweep_header(radar, isweep);
629 if (feof(wf->fptr) != 0) end_of_vos = 1;
630 } /* while not end of vos */