/*
* This file contains routines for processing Message Type 31, the digital
- * radar message type introduced in WSR-88D Level II Build 10.
+ * radar message type introduced in WSR-88D Level II Build 10. For more
+ * information, see the "Interface Control Document for the RDA/RPG" at the
+ * WSR-88D Radar Operations Center web site.
*/
#include "rsl.h"
unsigned char radial_spot_blanking;
unsigned char azm_indexing_mode;
unsigned short data_block_count;
- /* Data Block Pointers */
- unsigned int dbptr_vol_const;
- unsigned int dbptr_elev_const;
- unsigned int dbptr_radial_const;
- unsigned int dbptr_ref;
- unsigned int dbptr_vel;
- unsigned int dbptr_sw;
- unsigned int dbptr_zdr;
- unsigned int dbptr_phi;
- unsigned int dbptr_rho;
+ /* Data Block Indexes */
+ unsigned int vol_const;
+ unsigned int elev_const;
+ unsigned int radial_const;
+ unsigned int field1;
+ unsigned int field2;
+ unsigned int field3;
+ unsigned int field4;
+ unsigned int field5;
+ unsigned int field6;
} Ray_header_m31; /* Called Data Header Block in RDA/RPG document. */
typedef struct {
}
-void wsr88d_swap_m31_ray(Ray_header_m31 *wsr88d_ray)
+void wsr88d_swap_m31_ray_hdr(Ray_header_m31 *ray_hdr)
{
int *data_ptr;
- swap_4_bytes(&wsr88d_ray->ray_time);
- swap_2_bytes(&wsr88d_ray->ray_date);
- swap_2_bytes(&wsr88d_ray->azm_num);
- swap_4_bytes(&wsr88d_ray->azm);
- swap_2_bytes(&wsr88d_ray->radial_len);
- swap_4_bytes(&wsr88d_ray->elev);
- swap_2_bytes(&wsr88d_ray->data_block_count);
- data_ptr = (int *) &wsr88d_ray->dbptr_vol_const;
- for (; data_ptr <= (int *) &wsr88d_ray->dbptr_rho; data_ptr++)
+ swap_4_bytes(&ray_hdr->ray_time);
+ swap_2_bytes(&ray_hdr->ray_date);
+ swap_2_bytes(&ray_hdr->azm_num);
+ swap_4_bytes(&ray_hdr->azm);
+ swap_2_bytes(&ray_hdr->radial_len);
+ swap_4_bytes(&ray_hdr->elev);
+ swap_2_bytes(&ray_hdr->data_block_count);
+ data_ptr = (int *) &ray_hdr->vol_const;
+ for (; data_ptr <= (int *) &ray_hdr->field6; data_ptr++)
swap_4_bytes(data_ptr);
}
float value[13] = {0.043945, 0.08789, 0.17578, 0.35156, .70313, 1.40625,
2.8125, 5.625, 11.25, 22.5, 45., 90., 180.};
- /* find which bits are set and sum corresponding values to get angle. */
+ /* Find which bits are set and sum corresponding values to get angle. */
bitfield = bitfield >> 3; /* 3 least significant bits aren't used. */
for (i = 0; i < 13; i++) {
void get_wsr88d_unamb_and_nyq_vel(Wsr88d_ray_m31 *wsr88d_ray, float *unamb_rng,
float *nyq_vel)
{
- int dbptr, found, ray_hdr_len;
+ int dindex, found;
short nyq_vel_sh, unamb_rng_sh;
found = 0;
- ray_hdr_len = sizeof(wsr88d_ray->ray_hdr);
- dbptr = wsr88d_ray->ray_hdr.dbptr_radial_const - ray_hdr_len;
- if (strncmp((char *) &wsr88d_ray->data[dbptr], "RRAD", 4) == 0) found = 1;
+ dindex = wsr88d_ray->ray_hdr.radial_const;
+ if (strncmp((char *) &wsr88d_ray->data[dindex], "RRAD", 4) == 0) found = 1;
else {
- dbptr = wsr88d_ray->ray_hdr.dbptr_elev_const - ray_hdr_len;
- if (strncmp((char *) &wsr88d_ray->data[dbptr], "RRAD", 4) == 0)
+ dindex = wsr88d_ray->ray_hdr.elev_const;
+ if (strncmp((char *) &wsr88d_ray->data[dindex], "RRAD", 4) == 0)
found = 1;
else {
- dbptr = wsr88d_ray->ray_hdr.dbptr_vol_const - ray_hdr_len;
- if (strncmp((char *) &wsr88d_ray->data[dbptr], "RRAD", 4) == 0)
+ dindex = wsr88d_ray->ray_hdr.vol_const;
+ if (strncmp((char *) &wsr88d_ray->data[dindex], "RRAD", 4) == 0)
found = 1;
}
}
if (found) {
- memcpy(&unamb_rng_sh, &wsr88d_ray->data[dbptr+6], 2);
- memcpy(&nyq_vel_sh, &wsr88d_ray->data[dbptr+16], 2);
+ memcpy(&unamb_rng_sh, &wsr88d_ray->data[dindex+6], 2);
+ memcpy(&nyq_vel_sh, &wsr88d_ray->data[dindex+16], 2);
if (little_endian()) {
swap_2_bytes(&unamb_rng_sh);
swap_2_bytes(&nyq_vel_sh);
int read_wsr88d_ray_m31(Wsr88d_file *wf, int msg_size,
Wsr88d_ray_m31 *wsr88d_ray)
{
- Ray_header_m31 ray_hdr;
- int n, bytes_to_read, ray_hdr_len;
+ int n;
float nyq_vel, unamb_rng;
- /* Read ray data header block */
- n = fread(&wsr88d_ray->ray_hdr, sizeof(Ray_header_m31), 1, wf->fptr);
+ /* Read wsr88d ray. */
+
+ n = fread(wsr88d_ray->data, msg_size, 1, wf->fptr);
if (n < 1) {
fprintf(stderr,"read_wsr88d_ray_m31: Read failed.\n");
return 0;
}
- /* Byte swap header if needed. */
- if (little_endian()) wsr88d_swap_m31_ray(&wsr88d_ray->ray_hdr);
-
- ray_hdr = wsr88d_ray->ray_hdr;
- ray_hdr_len = sizeof(ray_hdr);
-
- /* Read data portion of radial. */
+ /* Copy data header block to ray header structure. */
+ memcpy(&wsr88d_ray->ray_hdr, &wsr88d_ray->data, sizeof(Ray_header_m31));
- bytes_to_read = msg_size - ray_hdr_len;
- n = fread(wsr88d_ray->data, bytes_to_read, 1, wf->fptr);
- if (n < 1) {
- fprintf(stderr,"read_wsr88d_ray_m31: Read failed.\n");
- return 0;
- }
+ if (little_endian()) wsr88d_swap_m31_ray_hdr(&wsr88d_ray->ray_hdr);
- /* We retrieve unambiguous range and Nyquist velocity here so that we don't
- * have to do it repeatedly for each data moment later.
+ /* Retrieve unambiguous range and Nyquist velocity here so that we don't
+ * have to do it for each data moment later.
*/
get_wsr88d_unamb_and_nyq_vel(wsr88d_ray, &unamb_rng, &nyq_vel);
wsr88d_ray->unamb_rng = unamb_rng;
}
-void wsr88d_load_ray_hdr(Wsr88d_ray_m31 wsr88d_ray, Ray *ray)
+void wsr88d_load_ray_hdr(Wsr88d_ray_m31 *wsr88d_ray, Ray *ray)
{
int month, day, year, hour, minute, sec;
float fsec;
Wsr88d_ray m1_ray;
Ray_header_m31 ray_hdr;
- ray_hdr = wsr88d_ray.ray_hdr;
+ ray_hdr = wsr88d_ray->ray_hdr;
m1_ray.ray_date = ray_hdr.ray_date;
m1_ray.ray_time = ray_hdr.ray_time;
ray->h.ray_num = ray_hdr.azm_num;
ray->h.elev = ray_hdr.elev;
ray->h.elev_num = ray_hdr.elev_num;
- ray->h.unam_rng = wsr88d_ray.unamb_rng;
- ray->h.nyq_vel = wsr88d_ray.nyq_vel;
+ ray->h.unam_rng = wsr88d_ray->unamb_rng;
+ ray->h.nyq_vel = wsr88d_ray->nyq_vel;
int elev_index;
elev_index = ray_hdr.elev_num - 1;
ray->h.azim_rate = vcp_data.azim_rate[elev_index];
*/
m1_ray.vol_cpat = vcp_data.vcp;
m1_ray.elev_num = ray_hdr.elev_num;
- m1_ray.unam_rng = (short) (wsr88d_ray.unamb_rng * 10.);
+ m1_ray.unam_rng = (short) (wsr88d_ray->unamb_rng * 10.);
/* Get values from message type 1 routines. */
ray->h.frequency = wsr88d_get_frequency(&m1_ray);
ray->h.pulse_width = wsr88d_get_pulse_width(&m1_ray);
int wsr88d_get_vol_index(char* dataname)
{
- int vol_index = -1;
-
- if (strncmp(dataname, "DREF", 4) == 0) vol_index = DZ_INDEX;
- if (strncmp(dataname, "DVEL", 4) == 0) vol_index = VR_INDEX;
- if (strncmp(dataname, "DSW", 3) == 0) vol_index = SW_INDEX;
- if (strncmp(dataname, "DZDR", 3) == 0) vol_index = DR_INDEX;
- if (strncmp(dataname, "DPHI", 3) == 0) vol_index = PH_INDEX;
- if (strncmp(dataname, "DRHO", 3) == 0) vol_index = RH_INDEX;
-
- return vol_index;
+ if (strncmp(dataname, "DREF", 4) == 0) return DZ_INDEX;
+ if (strncmp(dataname, "DVEL", 4) == 0) return VR_INDEX;
+ if (strncmp(dataname, "DSW", 3) == 0) return SW_INDEX;
+ if (strncmp(dataname, "DZDR", 4) == 0) return DR_INDEX;
+ if (strncmp(dataname, "DPHI", 4) == 0) return PH_INDEX;
+ if (strncmp(dataname, "DRHO", 4) == 0) return RH_INDEX;
+
+ return -1;
}
#define MAXRAYS_M31 800
#define MAXSWEEPS 20
-void wsr88d_load_ray(Wsr88d_ray_m31 wsr88d_ray, int data_ptr,
- int isweep, int iray, Radar *radar)
+void wsr88d_load_ray_into_radar(Wsr88d_ray_m31 *wsr88d_ray, int isweep,
+ Radar *radar)
{
- /* Load data into ray structure for this field or data moment. */
+ /* Load data into ray structure for each data field. */
+
+ int data_index;
+ int *field_offset;
+ int ifield, nfields;
+ int iray;
+ const nconstblocks = 3;
Data_moment_hdr data_hdr;
- int ngates;
+ int ngates, do_swap;
int i, hdr_size;
+ unsigned short item;
float value, scale, offset;
unsigned char *data;
Range (*invf)(float x);
float (*f)(Range x);
Ray *ray;
int vol_index, waveform;
+ char *type_str;
+
+ int keep_hi_prf_dz = 0; /* TODO: make this an argument. */
enum waveforms {surveillance=1, doppler_w_amb_res, doppler_no_amb_res,
batch};
- /* Get data moment header. */
- hdr_size = sizeof(data_hdr);
- memcpy(&data_hdr, &wsr88d_ray.data[data_ptr], hdr_size);
- data_ptr += hdr_size;
- if (little_endian()) wsr88d_swap_data_hdr(&data_hdr);
-
- vol_index = wsr88d_get_vol_index(data_hdr.dataname);
- if (vol_index < 0) {
- fprintf(stderr,"wsr88d_load_ray: Unknown dataname %s. isweep = %d, "
- "iray = %d.\n", data_hdr.dataname, isweep, iray);
- return;
- }
- switch (vol_index) {
- case DZ_INDEX: f = DZ_F; invf = DZ_INVF; break;
- case VR_INDEX: f = VR_F; invf = VR_INVF; break;
- case SW_INDEX: f = SW_F; invf = SW_INVF; break;
- case DR_INDEX: f = DR_F; invf = DR_INVF; break;
- case PH_INDEX: f = PH_F; invf = PH_INVF; break;
- case RH_INDEX: f = RH_F; invf = RH_INVF; break;
- }
+ nfields = wsr88d_ray->ray_hdr.data_block_count - nconstblocks;
+ field_offset = (int *) &wsr88d_ray->ray_hdr.radial_const;
+ do_swap = little_endian();
+ iray = wsr88d_ray->ray_hdr.azm_num - 1;
+
+ for (ifield=0; ifield < nfields; ifield++) {
+ field_offset++;
+ data_index = *field_offset;
+ /* Get data moment header. */
+ hdr_size = sizeof(data_hdr);
+ memcpy(&data_hdr, &wsr88d_ray->data[data_index], hdr_size);
+ if (do_swap) wsr88d_swap_data_hdr(&data_hdr);
+ data_index += hdr_size;
+
+ vol_index = wsr88d_get_vol_index(data_hdr.dataname);
+ if (vol_index < 0) {
+ fprintf(stderr,"wsr88d_load_ray_into_radar: Unknown dataname %s. "
+ "isweep = %d, iray = %d.\n", data_hdr.dataname, isweep,
+ iray);
+ return;
+ }
+ switch (vol_index) {
+ case DZ_INDEX: f = DZ_F; invf = DZ_INVF;
+ type_str = "Reflectivity"; break;
+ case VR_INDEX: f = VR_F; invf = VR_INVF;
+ type_str = "Velocity"; break;
+ case SW_INDEX: f = SW_F; invf = SW_INVF;
+ type_str = "Spectrum width"; break;
+ case DR_INDEX: f = DR_F; invf = DR_INVF;
+ type_str = "Diff. Reflectivity"; break;
+ case PH_INDEX: f = PH_F; invf = PH_INVF;
+ type_str = "Diff. Phase"; break;
+ case RH_INDEX: f = RH_F; invf = RH_INVF;
+ type_str = "Correlation Coef (Rho)"; break;
+ }
- waveform = vcp_data.waveform[isweep];
- /* The following somewhat complicated if-statement says we'll do the
- * body of the statement if:
- * a) the data moment is not reflectivity, or
- * b) the data moment is reflectivity and one of the following is true:
- * - waveform is surveillance
- * - waveform is batch
- * - elevation is greater than highest split cut (i.e., 6 deg.)
- */
- if (vol_index != DZ_INDEX || (waveform == surveillance ||
- waveform == batch || vcp_data.fixed_angle[isweep] >= 6.0)) {
+ waveform = vcp_data.waveform[isweep];
+
+ /* Ignore short-range reflectivity from velocity split cuts unless
+ * merging of split cuts is suppressed. The indicators for this type of
+ * reflectivity are surveillance mode is 0 and elevation angle is
+ * below 6 degrees.
+ */
+ if (vol_index == DZ_INDEX && (vcp_data.surveil_prf_num[isweep] == 0 &&
+ vcp_data.fixed_angle[isweep] < 6.0 && !keep_hi_prf_dz))
+ continue;
+
+ /* Load the data for this field. */
if (radar->v[vol_index] == NULL) {
radar->v[vol_index] = RSL_new_volume(MAXSWEEPS);
radar->v[vol_index]->h.f = f;
radar->v[vol_index]->h.invf = invf;
+ radar->v[vol_index]->h.type_str = type_str;
}
if (radar->v[vol_index]->sweep[isweep] == NULL) {
radar->v[vol_index]->sweep[isweep] = RSL_new_sweep(MAXRAYS_M31);
offset = data_hdr.offset;
scale = data_hdr.scale;
if (data_hdr.scale == 0) scale = 1.0;
- data = &wsr88d_ray.data[data_ptr];
+ data = &wsr88d_ray->data[data_index];
for (i = 0; i < ngates; i++) {
- if (data[i] > 1)
- value = (data[i] - offset) / scale;
- else value = (data[i] == 0) ? BADVAL : RFVAL;
+ if (data_hdr.datasize_bits != 16) {
+ item = *data;
+ data++;
+ } else {
+ item = *(unsigned short *)data;
+ if (do_swap) swap_2_bytes(&item);
+ data += 2;
+ }
+ if (item > 1)
+ value = (item - offset) / scale;
+ else value = (item == 0) ? BADVAL : RFVAL;
ray->range[i] = invf(value);
ray->h.f = f;
ray->h.invf = invf;
ray->h.nbins = ngates;
radar->v[vol_index]->sweep[isweep]->ray[iray] = ray;
radar->v[vol_index]->sweep[isweep]->h.nrays = iray+1;
- }
-}
-
-
-void wsr88d_load_ray_into_radar(Wsr88d_ray_m31 wsr88d_ray, int isweep, int iray,
- Radar *radar)
-{
- int *data_ptr, hdr_size;
- int i, ndatablocks, nconstblocks = 3;
-
- hdr_size = sizeof(wsr88d_ray.ray_hdr);
-
- ndatablocks = wsr88d_ray.ray_hdr.data_block_count;
- data_ptr = (int *) &wsr88d_ray.ray_hdr.dbptr_ref;
- for (i=0; i < ndatablocks-nconstblocks; i++) {
- wsr88d_load_ray(wsr88d_ray, *data_ptr-hdr_size, isweep, iray, radar);
- data_ptr++;
- }
+ } /* for each data field */
}
}
-Radar *load_wsr88d_m31_into_radar(Wsr88d_file *wf)
+Radar *wsr88d_load_m31_into_radar(Wsr88d_file *wf)
{
Wsr88d_msg_hdr msghdr;
Wsr88d_ray_m31 wsr88d_ray;
short non31_seg_remainder[1202]; /* Remainder after message header */
- int end_of_vos = 0, isweep = 0, iray = 0;
+ int end_of_vos = 0, isweep = 0;
int msg_hdr_size, msg_size, n;
+ int sweep_hdrs_written = 0, prev_elev_num = 1, prev_raynum = 0, raynum = 0;
Radar *radar = NULL;
- /* Message type 31 is a variable length message. All other message types
- * are made up of 2432-byte segments. To handle all types, we read the
- * message header and check the message type. If it is not 31, we simply
- * read the remainder of the 2432-byte segment. If message type is 31, we
- * use the size given in message header to determine how many bytes to
- * read. For more information, see the "Interface Control Document for the
- * RDA/RPG" at the WSR-88D Radar Operations Center web site.
- */
+ /* Message type 31 is a variable length message. All other types are made
+ * up of 1 or more segments, where each segment is 2432 bytes in length.
+ * To handle these differences, read the message header and check its type.
+ * If it is 31, use the size given in the message header to determine the
+ * number of bytes to read. If not, simply read the remainder of the
+ * 2432-byte segment.
+ */
n = fread(&msghdr, sizeof(Wsr88d_msg_hdr), 1, wf->fptr);
n = read_wsr88d_ray_m31(wf, msg_size, &wsr88d_ray);
/* Assume error message was issued from read routine */
if (n <= 0) return NULL;
+ raynum = wsr88d_ray.ray_hdr.azm_num;
+ if (raynum > MAXRAYS_M31) {
+ fprintf(stderr,"Error: raynum = %d, exceeds MAXRAYS_M31"
+ " (%d)\n", raynum, MAXRAYS_M31);
+ fprintf(stderr,"isweep = %d\n", isweep);
+ RSL_free_radar(radar);
+ return NULL;
+ }
- /* We need to check radial status for start of new elevation.
- * Sometimes this occurs without an end-of-elevation flag for the
- * previous sweep, which is the trigger for loading the sweep
- * header. "iray" is set to zero after end-of-elev is received,
- * so that's why we check it. We issue a warning because when this
- * occurs, the number of rays in the previous sweep is usually less
- * than expected.
+ /* Check for an unexpected start of new elevation, and issue a
+ * warning if this has occurred. This usually means less rays than
+ * expected. It happens, but rarely.
*/
if (wsr88d_ray.ray_hdr.radial_status == START_OF_ELEV &&
- iray != 0) {
+ sweep_hdrs_written != prev_elev_num) {
fprintf(stderr,"Warning: Radial status is Start-of-Elevation, "
"but End-of-Elevation was not\n"
"issued for elevation number %d. Number of rays = %d"
- "\n", isweep+1, iray);
+ "\n", prev_elev_num, prev_raynum);
wsr88d_load_sweep_header(radar, isweep);
isweep++;
- iray = 0;
+ sweep_hdrs_written++;
+ prev_elev_num = wsr88d_ray.ray_hdr.elev_num - 1;
}
/* Load ray into radar structure. */
- wsr88d_load_ray_into_radar(wsr88d_ray, isweep, iray, radar);
- iray++;
- if (iray >= MAXRAYS_M31) {
- fprintf(stderr,"Error: iray = %d, equals or exceeds MAXRAYS_M31"
- " (%d)\n", iray, MAXRAYS_M31);
- fprintf(stderr,"isweep = %d\n", isweep);
- RSL_free_radar(radar);
- return NULL;
- }
+ wsr88d_load_ray_into_radar(&wsr88d_ray, isweep, radar);
+ prev_raynum = raynum;
}
else { /* msg_type not 31 */
n = fread(&non31_seg_remainder, sizeof(non31_seg_remainder), 1,
else
fprintf(stderr,"Read failed.\n");
fprintf(stderr,"Current sweep index: %d\n"
- "Last ray read: %d\n", isweep, iray);
+ "Last ray read: %d\n", isweep, prev_raynum);
wsr88d_load_sweep_header(radar, isweep);
return radar;
}
if (msghdr.msg_type == 5) {
wsr88d_get_vcp_data(non31_seg_remainder);
radar->h.vcp = vcp_data.vcp;
- /* printf("VCP = %d\n", vcp_data.vcp); */
}
}
if (wsr88d_ray.ray_hdr.radial_status == END_OF_ELEV) {
wsr88d_load_sweep_header(radar, isweep);
isweep++;
- iray = 0;
+ sweep_hdrs_written++;
+ prev_elev_num = wsr88d_ray.ray_hdr.elev_num;
}
/* If not at end of volume scan, read next message header. */
"Unexpected end of file.\n");
else fprintf(stderr,"Failed reading msghdr.\n");
fprintf(stderr,"Current sweep index: %d\n"
- "Last ray read: %d\n", isweep, iray);
+ "Last ray read: %d\n", isweep, prev_raynum);
wsr88d_load_sweep_header(radar, isweep);
return radar;
}
projects / ~andy / rsl / blobdiff