X-Git-Url: http://pileus.org/git/?p=~andy%2Frsl;a=blobdiff_plain;f=wsr88d_m31.c;h=9e4c5b198603b7172bafd659ed092e1dba08ec44;hp=4373df3220f8033fa034f4056581cd0a9bed7891;hb=71c71381ace86eb987172b22de55fbb837518b26;hpb=d08a7f8a699a044bc4ac5f93917aa7f6c463923b diff --git a/wsr88d_m31.c b/wsr88d_m31.c index 4373df3..9e4c5b1 100644 --- a/wsr88d_m31.c +++ b/wsr88d_m31.c @@ -23,7 +23,9 @@ /* * 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" @@ -64,16 +66,16 @@ typedef struct { 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 { @@ -115,19 +117,19 @@ void wsr88d_swap_m31_hdr(Wsr88d_msg_hdr *msghdr) } -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); } @@ -151,7 +153,7 @@ float wsr88d_get_angle(short bitfield) 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++) { @@ -242,26 +244,25 @@ void wsr88d_get_vcp_data(short *msgtype5) 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); @@ -278,34 +279,24 @@ void get_wsr88d_unamb_and_nyq_vel(Wsr88d_ray_m31 *wsr88d_ray, float *unamb_rng, 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; @@ -315,14 +306,14 @@ int read_wsr88d_ray_m31(Wsr88d_file *wf, int msg_size, } -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; @@ -338,8 +329,8 @@ void wsr88d_load_ray_hdr(Wsr88d_ray_m31 wsr88d_ray, Ray *ray) 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]; @@ -354,7 +345,7 @@ void wsr88d_load_ray_hdr(Wsr88d_ray_m31 wsr88d_ray, Ray *ray) */ 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); @@ -366,76 +357,101 @@ void wsr88d_load_ray_hdr(Wsr88d_ray_m31 wsr88d_ray, Ray *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); @@ -453,11 +469,19 @@ void wsr88d_load_ray(Wsr88d_ray_m31 wsr88d_ray, int data_ptr, 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; @@ -468,24 +492,7 @@ void wsr88d_load_ray(Wsr88d_ray_m31 wsr88d_ray, int data_ptr, 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 */ } @@ -512,23 +519,23 @@ void wsr88d_load_sweep_header(Radar *radar, int isweep) } -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); @@ -549,36 +556,34 @@ Radar *load_wsr88d_m31_into_radar(Wsr88d_file *wf) 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, @@ -590,14 +595,13 @@ Radar *load_wsr88d_m31_into_radar(Wsr88d_file *wf) 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); */ } } @@ -605,7 +609,8 @@ Radar *load_wsr88d_m31_into_radar(Wsr88d_file *wf) 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. */ @@ -617,7 +622,7 @@ Radar *load_wsr88d_m31_into_radar(Wsr88d_file *wf) "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; }