#include "rsl.h"
#include "wsr88d.h"
+#include <string.h>
/* Data descriptions in the following data structures are from the "Interface
* Control Document for the RDA/RPG", Build 10.0 Draft, WSR-88D Radar
*/
typedef struct {
- short rpg[6]; /* Unused. Not really part of message header, but is
- * inserted by RPG Communications Manager for each message.
- */
- unsigned short msg_size; /* for this segment, in halfwords */
+ short rpg[6]; /* 12 bytes inserted by RPG Communications Mgr. Ignored. */
+ unsigned short msg_size; /* Message size for this segment, in halfwords */
unsigned char channel; /* RDA Redundant Channel */
- unsigned char msg_type; /* Message Type */
- unsigned short id_seq; /* I.d. Seq = 0 to 7FFF, then roll over to 0 */
- unsigned short msg_date; /* modified Julian date from 1/1/70 */
- unsigned int msg_time; /* packet generation time in ms past midnite */
- unsigned short num_segs;
- unsigned short seg_num;
+ unsigned char msg_type; /* Message type. For example, 31 */
+ unsigned short id_seq; /* Msg seq num = 0 to 7FFF, then roll over to 0 */
+ unsigned short msg_date; /* Modified Julian date from 1/1/70 */
+ unsigned int msg_time; /* Packet generation time in ms past midnight */
+ unsigned short num_segs; /* Number of segments for this message */
+ unsigned short seg_num; /* Number of this segment */
} Wsr88d_msg_hdr;
typedef struct {
float offset;
} Data_moment_hdr;
-typedef struct {
- Data_moment_hdr data_hdr;
- /* TODO: data type will need to be changed to "void" to handle Dual Pole:
- * some data are in 2-byte words.
- */
- unsigned char *data;
-} Data_moment;
+#define MAX_RADIAL_LENGTH 14288
typedef struct {
Ray_header_m31 ray_hdr;
- short unamb_rng;
- short nyq_vel;
- /* Pointers to data */
- /* TODO: Replace data moment pointers with single pointer to radial.
- * Can simply make unsigned char array as in MSG1 version.
- * Maximum radial length is 14288 bytes.
- */
- Data_moment *ref;
- Data_moment *vel;
- Data_moment *sw;
+ float unamb_rng;
+ float nyq_vel;
+ unsigned char data[MAX_RADIAL_LENGTH];
} Wsr88d_ray_m31;
-#define MAXRAYS_M31 800
-#define MAXSWEEPS 20
-enum {START_OF_ELEV, INTERMED_RADIAL, END_OF_ELEV, BEGIN_VOS, END_VOS};
+enum radial_status {START_OF_ELEV, INTERMED_RADIAL, END_OF_ELEV, BEGIN_VOS,
+ END_VOS};
void wsr88d_swap_m31_hdr(Wsr88d_msg_hdr *msghdr)
void wsr88d_swap_m31_ray(Ray_header_m31 *wsr88d_ray)
{
- int *fullword;
+ int *data_ptr;
swap_4_bytes(&wsr88d_ray->ray_time);
swap_2_bytes(&wsr88d_ray->ray_date);
swap_2_bytes(&wsr88d_ray->radial_len);
swap_4_bytes(&wsr88d_ray->elev);
swap_2_bytes(&wsr88d_ray->data_block_count);
- fullword = (int *) &wsr88d_ray->dbptr_vol_const;
- for (; fullword <= (int *) &wsr88d_ray->dbptr_rho; fullword++)
- swap_4_bytes(fullword);
+ data_ptr = (int *) &wsr88d_ray->dbptr_vol_const;
+ for (; data_ptr <= (int *) &wsr88d_ray->dbptr_rho; data_ptr++)
+ swap_4_bytes(data_ptr);
}
-void wsr88d_swap_data_moment(Data_moment_hdr *this_field)
+void wsr88d_swap_data_hdr(Data_moment_hdr *this_field)
{
short *halfword;
halfword = (short *) &this_field->ngates;
}
-void testprt(Ray_header_m31 ray_hdr)
-{
- /* For testing: Print some values from data block header. */
- printf("\nray_time: %d\n",ray_hdr.ray_time);
- printf("ray_date: %d\n",ray_hdr.ray_date);
- printf("azm_num: %d\n",ray_hdr.azm_num);
- printf("azm: %f\n",ray_hdr.azm);
- printf("radial_len: %d\n",ray_hdr.radial_len);
- printf("elev: %f\n",ray_hdr.elev);
- printf("data block count: %d\n",ray_hdr.data_block_count);
- printf("dbptr_vol_const: %d\n",ray_hdr.dbptr_vol_const);
- printf("dbptr_elev_const: %d\n",ray_hdr.dbptr_elev_const);
- printf("dbptr_radial_const: %d\n",ray_hdr.dbptr_radial_const);
- printf("dbptr_ref: %d\n",ray_hdr.dbptr_ref);
- printf("dbptr_vel: %d\n",ray_hdr.dbptr_vel);
- printf("dbptr_sw: %d\n",ray_hdr.dbptr_sw);
-}
-
-
float wsr88d_get_angle(short bitfield)
{
short mask = 1;
}
-Data_moment *read_data_moment(unsigned char *radial, int begin_data_block)
+void get_wsr88d_unamb_and_nyq_vel(Wsr88d_ray_m31 *wsr88d_ray, float *unamb_rng,
+ float *nyq_vel)
{
- Data_moment *this_field;
- unsigned char* data;
- unsigned char* dptr;
- int hdr_size, n;
-
- /*
- printf("read_data_moment: ");
- printf("sizeof(Data_moment) = %d\n", sizeof(Data_moment));
- printf("sizeof(Data_moment_hdr) = %d\n", sizeof(Data_moment_hdr));
- */
- this_field = malloc(sizeof(Data_moment));
- hdr_size = sizeof(Data_moment_hdr);
- /* printf("hdr_size = %d\n", hdr_size); */
-
- dptr = radial + begin_data_block;
- memcpy(&this_field->data_hdr, dptr, hdr_size);
- dptr += hdr_size;
-
- if (little_endian()) wsr88d_swap_data_moment(&this_field->data_hdr);
- data = (unsigned char *) malloc(this_field->data_hdr.ngates);
- memcpy(data, dptr, this_field->data_hdr.ngates);
- this_field->data = data;
- return this_field;
+ int dbptr, found, ray_hdr_len;
+ 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;
+ else {
+ dbptr = wsr88d_ray->ray_hdr.dbptr_elev_const - ray_hdr_len;
+ if (strncmp((char *) &wsr88d_ray->data[dbptr], "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)
+ found = 1;
+ }
+ }
+ if (found) {
+ memcpy(&unamb_rng_sh, &wsr88d_ray->data[dbptr+6], 2);
+ memcpy(&nyq_vel_sh, &wsr88d_ray->data[dbptr+16], 2);
+ if (little_endian()) {
+ swap_2_bytes(&unamb_rng_sh);
+ swap_2_bytes(&nyq_vel_sh);
+ }
+ *unamb_rng = unamb_rng_sh / 10.;
+ *nyq_vel = nyq_vel_sh / 100.;
+ } else {
+ *unamb_rng = 0.;
+ *nyq_vel = 0.;
+ }
}
Wsr88d_ray_m31 *wsr88d_ray)
{
Ray_header_m31 ray_hdr;
- int n, begin_data_block, bytes_to_read, ray_hdr_len;
- unsigned char *radial;
+ int n, bytes_to_read, ray_hdr_len;
+ float nyq_vel, unamb_rng;
/* Read ray data header block */
n = fread(&wsr88d_ray->ray_hdr, sizeof(Ray_header_m31), 1, wf->fptr);
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;
-
- /*
- if (ray_hdr.azm_num == 1) testprt(wsr88d_ray->ray_hdr);
- */
-
ray_hdr_len = sizeof(ray_hdr);
- /* Read in radial
- * Use header_size offset with data pointers
- * Copy data into structures
- */
+ /* Read data portion of radial. */
bytes_to_read = msg_size - ray_hdr_len;
- radial = (unsigned char *) malloc(bytes_to_read);
- n = fread(radial, bytes_to_read, 1, wf->fptr);
+ 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 (ray_hdr.dbptr_radial_const != 0) {
- begin_data_block = ray_hdr.dbptr_radial_const - ray_hdr_len;
- memcpy(&wsr88d_ray->unamb_rng, &radial[begin_data_block+6], 2);
- memcpy(&wsr88d_ray->nyq_vel, &radial[begin_data_block+16], 2);
- if (little_endian()) {
- swap_2_bytes(&wsr88d_ray->unamb_rng);
- swap_2_bytes(&wsr88d_ray->nyq_vel);
- }
- }
- else {
- wsr88d_ray->unamb_rng = 0;
- wsr88d_ray->nyq_vel = 0;
- }
-
- if (ray_hdr.dbptr_ref != 0) {
- begin_data_block = ray_hdr.dbptr_ref - ray_hdr_len;
- wsr88d_ray->ref = read_data_moment(radial, begin_data_block);
- }
- if (ray_hdr.dbptr_vel != 0) {
- begin_data_block = ray_hdr.dbptr_vel - ray_hdr_len;
- wsr88d_ray->vel = read_data_moment(radial, begin_data_block);
- }
- if (ray_hdr.dbptr_sw != 0) {
- begin_data_block = ray_hdr.dbptr_sw - ray_hdr_len;
- wsr88d_ray->sw = read_data_moment(radial, begin_data_block);
- }
-
-
- /* For testing: print reflectivity data. */
- int prtdata = 0;
- { int i;
- if (prtdata) {
- for (i=0; i < wsr88d_ray->ref->data_hdr.ngates; i++) {
- if (i % 10 == 0) printf("\n");
- printf(" %d", wsr88d_ray->ref->data[i]);
- }
- printf("\n");
- }
- }
-
- free(radial);
- return 1;
-}
-
-
-void wsr88d_load_ray_data(Data_moment *data_block, Ray *ray)
-{
- Data_moment_hdr data_hdr;
- int ngates;
- int i;
- float value, scale, offset;
- unsigned char *data;
- Range (*invf)(float x);
- float (*f)(Range x);
-
- data_hdr = data_block->data_hdr;
- data = data_block->data;
-
- /*
- printf("wsr88d_load_ray_data: ");
- printf(" DataName: %s\n", data_hdr.dataname);
- */
-
- ngates = data_hdr.ngates;
- /*
- printf(" ngates = %d\n", ngates);
- printf(" scale = %f\n", data_hdr.scale);
- printf(" offset = %f\n", data_hdr.offset);
- */
- offset = data_hdr.offset;
- scale = data_hdr.scale;
-
- /* Note: data range is 2-255. 0 means signal is below threshold, and 1
- * means range folded.
+ /* We retrieve unambiguous range and Nyquist velocity here so that we don't
+ * have to do it repeatedly for each data moment later.
*/
+ get_wsr88d_unamb_and_nyq_vel(wsr88d_ray, &unamb_rng, &nyq_vel);
+ wsr88d_ray->unamb_rng = unamb_rng;
+ wsr88d_ray->nyq_vel = nyq_vel;
- /* Test print
- for (i=0; i < 10; i++) {
- value = (data[i] - offset) / scale;
- printf(" bytevalue = %d, value = %f\n", data[i], value);
- }
- */
-
- /* Reflectivity */
- if (strncmp(data_hdr.dataname, "DREF", 4) == 0) {
- /* convert data to float
- * convert float to range and put in ray.range
- */
- f = DZ_F;
- invf = DZ_INVF;
- for (i = 0; i < ngates; i++) {
- if (data[i] > 1)
- value = (data[i] - offset) / scale;
- else value = (data[i] == 0) ? BADVAL : RFVAL;
- ray->range[i] = invf(value);
- ray->h.f = f;
- ray->h.invf = invf;
- }
- }
-
- /* Velocity */
- if (strncmp(data_hdr.dataname, "DVEL", 4) == 0) {
- /* convert data to float
- * convert float to range and put in ray.range
- */
- f = VR_F;
- invf = VR_INVF;
- for (i = 0; i < ngates; i++) {
- if (data[i] > 1)
- value = (data[i] - offset) / scale;
- else value = (data[i] == 0) ? BADVAL : RFVAL;
- ray->range[i] = invf(value);
- ray->h.f = f;
- ray->h.invf = invf;
- }
- }
-
- /* Spectrum Width */
- if (strncmp(data_hdr.dataname, "DSW", 3) == 0) {
- /* convert data to float
- * convert float to range and put in ray.range
- */
- f = SW_F;
- invf = SW_INVF;
- for (i = 0; i < ngates; i++) {
- if (data[i] > 1)
- value = (data[i] - offset) / scale;
- else value = (data[i] == 0) ? BADVAL : RFVAL;
- ray->range[i] = invf(value);
- ray->h.f = f;
- ray->h.invf = invf;
- }
- }
- ray->h.range_bin1 = data_hdr.range_first_gate;
- ray->h.gate_size = data_hdr.range_samp_interval;
- ray->h.nbins = ngates;
+ return 1;
}
-
void wsr88d_load_ray_hdr(Wsr88d_ray_m31 wsr88d_ray, Ray *ray)
{
int month, day, year, hour, minute, sec;
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 / 10.;
- ray->h.nyq_vel = wsr88d_ray.nyq_vel / 100.;
+ 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];
ray->h.fix_angle = vcp_data.fixed_angle[elev_index];
ray->h.vel_res = vcp_data.vel_res;
+ if (ray_hdr.azm_res != 1)
+ ray->h.beam_width = 1.0;
+ else ray->h.beam_width = 0.5;
- /* Get some values using message type 1 routines.
- * First load VCP and elevation numbers into msg 1 ray.
+ /* For convenience, use message type 1 routines to get some values.
+ * First load VCP and elevation numbers into a msg 1 ray.
*/
m1_ray.vol_cpat = vcp_data.vcp;
m1_ray.elev_num = ray_hdr.elev_num;
- m1_ray.unam_rng = wsr88d_ray.unamb_rng;
- if (ray_hdr.azm_res != 1)
- ray->h.beam_width = 1.0;
- else ray->h.beam_width = 0.5;
+ 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);
ray->h.pulse_count = wsr88d_get_pulse_count(&m1_ray);
- ray->h.prf = wsr88d_get_prf(&m1_ray);
+ ray->h.prf = (int) wsr88d_get_prf(&m1_ray);
ray->h.wavelength = 0.1071;
}
int wsr88d_get_vol_index(char* dataname)
{
- int vol_index;
+ 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;
- /* TODO: Add the other data moments. */
+ 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;
}
-void wsr88d_load_ray_into_radar(Wsr88d_ray_m31 wsr88d_ray, int isweep, int iray,
- Radar *radar)
+#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)
{
- Ray *ray;
+ /* Load data into ray structure for this field or data moment. */
+
+ Data_moment_hdr data_hdr;
+ int ngates;
+ int i, hdr_size;
+ float value, scale, offset;
+ unsigned char *data;
Range (*invf)(float x);
float (*f)(Range x);
+ Ray *ray;
int vol_index, waveform;
- enum waveforms {surveillance=1, doppler_ambres, doppler_no_ambres, batch};
+ enum waveforms {surveillance=1, doppler_w_amb_res, doppler_no_amb_res,
+ batch};
-/* for each data moment do
- get the name of the data moment
- get new volume if needed
- get new sweep if needed
- get new ray
- put this data moment in ray
- endfor each data moment
-*/
-
- /* Note: The data moment type can only be determined by the name within the
- * individual data moment block. The name of the pointer is not reliable.
- * For example, in legacy resolution mode, dbptr_ref points to the velocity
- * data moment in velocity split cuts. Apparently the location of the first
- * data moment is always stored in the reflectivity pointer (dbptr_ref), and
- * sometimes this is velocity. When this occurs, the velocity data pointer
- * (dbptr_vel) then points to spectrum width.
- * With super resolution, there actually is reflectivity in the velocity
- * split cut. It's there for use with the recombination algorithm.
- */
+ /* 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);
- if (wsr88d_ray.ray_hdr.dbptr_ref > 0) {
- vol_index = wsr88d_get_vol_index(wsr88d_ray.ref->data_hdr.dataname);
- 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;
- default: f = DZ_F; invf = DZ_INVF; break;
- }
- /* If this is reflectivity, check the waveform type to make sure
- * it isn't from a Doppler split cut.
- * We only keep reflectivity if the waveform type is surveillance or
- * batch, or the elevation is above the split cut elevations.
- */
- waveform = vcp_data.waveform[isweep];
- if (vol_index != DZ_INDEX ||
- (waveform == surveillance || waveform == batch ||
- vcp_data.fixed_angle[isweep] >= 6.0))
- {
- 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;
- }
- if (radar->v[vol_index]->sweep[isweep] == NULL) {
- radar->v[vol_index]->sweep[isweep] = RSL_new_sweep(MAXRAYS_M31);
- radar->v[vol_index]->sweep[isweep]->h.f = f;
- radar->v[vol_index]->sweep[isweep]->h.invf = invf;
- }
- ray = RSL_new_ray(wsr88d_ray.ref->data_hdr.ngates);
- wsr88d_load_ray_data(wsr88d_ray.ref, ray);
- wsr88d_load_ray_hdr(wsr88d_ray, ray);
- radar->v[vol_index]->sweep[isweep]->ray[iray] = ray;
- radar->v[vol_index]->sweep[isweep]->h.nrays = iray+1;
- }
+ 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;
}
- if (wsr88d_ray.ray_hdr.dbptr_vel > 0) {
- vol_index = wsr88d_get_vol_index(wsr88d_ray.vel->data_hdr.dataname);
- 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;
- default: f = DZ_F; invf = DZ_INVF; 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)) {
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]->sweep[isweep]->h.f = f;
radar->v[vol_index]->sweep[isweep]->h.invf = invf;
}
- ray = RSL_new_ray(wsr88d_ray.vel->data_hdr.ngates);
- wsr88d_load_ray_data(wsr88d_ray.vel, ray);
- wsr88d_load_ray_hdr(wsr88d_ray, ray);
- radar->v[vol_index]->sweep[isweep]->ray[iray] = ray;
- radar->v[vol_index]->sweep[isweep]->h.nrays = iray+1;
- }
+ ngates = data_hdr.ngates;
+ ray = RSL_new_ray(ngates);
- if (wsr88d_ray.ray_hdr.dbptr_sw > 0) {
- vol_index = wsr88d_get_vol_index(wsr88d_ray.sw->data_hdr.dataname);
- 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;
- default: f = DZ_F; invf = DZ_INVF; break;
- }
- 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;
- }
- if (radar->v[vol_index]->sweep[isweep] == NULL) {
- radar->v[vol_index]->sweep[isweep] = RSL_new_sweep(MAXRAYS_M31);
- radar->v[vol_index]->sweep[isweep]->h.f = f;
- radar->v[vol_index]->sweep[isweep]->h.invf = invf;
+ /* Convert data to float, then use range function to store in ray.
+ * Note: data range is 2-255. 0 means signal is below threshold, and 1
+ * means range folded.
+ */
+
+ offset = data_hdr.offset;
+ scale = data_hdr.scale;
+ if (data_hdr.scale == 0) scale = 1.0;
+ data = &wsr88d_ray.data[data_ptr];
+ for (i = 0; i < ngates; i++) {
+ if (data[i] > 1)
+ value = (data[i] - offset) / scale;
+ else value = (data[i] == 0) ? BADVAL : RFVAL;
+ ray->range[i] = invf(value);
+ ray->h.f = f;
+ ray->h.invf = invf;
}
- ray = RSL_new_ray(wsr88d_ray.sw->data_hdr.ngates);
- wsr88d_load_ray_data(wsr88d_ray.sw, ray);
wsr88d_load_ray_hdr(wsr88d_ray, ray);
+ ray->h.range_bin1 = data_hdr.range_first_gate;
+ ray->h.gate_size = data_hdr.range_samp_interval;
+ 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++;
+ }
+}
+
-void wsr88d_load_sweep_header(Radar *radar, int isweep,
- Wsr88d_ray_m31 wsr88d_ray)
+void wsr88d_load_sweep_header(Radar *radar, int isweep)
{
- int ivolume;
- Ray_header_m31 ray_hdr;
+ int ivolume, nrays;
Sweep *sweep;
- int vcp;
-
- ray_hdr = wsr88d_ray.ray_hdr;
+ Ray *last_ray;
for (ivolume=0; ivolume < MAX_RADAR_VOLUMES; ivolume++) {
- if (radar->v[ivolume] != NULL && radar->v[ivolume]->sweep[isweep] != NULL) {
+ if (radar->v[ivolume] != NULL &&
+ radar->v[ivolume]->sweep[isweep] != NULL) {
sweep = radar->v[ivolume]->sweep[isweep];
- radar->v[ivolume]->sweep[isweep]->h.sweep_num = ray_hdr.elev_num;
- radar->v[ivolume]->sweep[isweep]->h.elev =
- vcp_data.fixed_angle[isweep];
- if (ray_hdr.azm_res != 1)
- sweep->h.beam_width = 1.0;
- else sweep->h.beam_width = 0.5;
+ nrays = sweep->h.nrays;
+ if (nrays == 0) continue;
+ last_ray = sweep->ray[nrays-1];
+ sweep->h.sweep_num = last_ray->h.elev_num;
+ sweep->h.elev = vcp_data.fixed_angle[isweep];
+ sweep->h.beam_width = last_ray->h.beam_width;
sweep->h.vert_half_bw = sweep->h.beam_width / 2.;
sweep->h.horz_half_bw = sweep->h.beam_width / 2.;
}
int msg_hdr_size, msg_size, n;
Radar *radar = NULL;
-/* Message type 31 is a variable length message, whereas all other message
- * types are made up of 2432 byte segments. To handle this, we read the
- * message header and check the message type. If it is not 31, then we read
- * the remainder of the constant size segment. If message type is 31, we read
- * the remainder of the message by the size given in message header.
- * When reading the message header, we must include 12 bytes inserted
- * by RPG, which we ignore, followed by the 8 halfwords (short integers) which
- * make up the actual message header.
- * 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 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.
+ */
n = fread(&msghdr, sizeof(Wsr88d_msg_hdr), 1, wf->fptr);
/* printf("msgtype = %d\n", msghdr.msg_type); */
msg_hdr_size = sizeof(Wsr88d_msg_hdr) - sizeof(msghdr.rpg);
-
radar = RSL_new_radar(MAX_RADAR_VOLUMES);
-
+
while (! end_of_vos) {
if (msghdr.msg_type == 31) {
if (little_endian()) wsr88d_swap_m31_hdr(&msghdr);
- /* Get size in bytes of message following message header.
- * The size given in message header is in halfwords, so double it.
+ /* Get size of the remainder of message. The given size is in
+ * halfwords, but we want it in bytes, so double it.
*/
msg_size = (int) msghdr.msg_size * 2 - msg_hdr_size;
n = read_wsr88d_ray_m31(wf, msg_size, &wsr88d_ray);
+ /* Assume error message was issued from read routine */
if (n <= 0) return NULL;
- /* Load this ray into radar structure ray. */
+ /* 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.
+ */
+ if (wsr88d_ray.ray_hdr.radial_status == START_OF_ELEV &&
+ iray != 0) {
+ 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);
+ wsr88d_load_sweep_header(radar, isweep);
+ isweep++;
+ iray = 0;
+ }
+
+ /* 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;
+ }
}
else { /* msg_type not 31 */
n = fread(&non31_seg_remainder, sizeof(non31_seg_remainder), 1,
fprintf(stderr,"Warning: load_wsr88d_m31_into_radar: ");
if (feof(wf->fptr) != 0)
fprintf(stderr, "Unexpected end of file.\n");
- else fprintf(stderr,"Read failed.\n");
- fprintf(stderr,"Current sweep number: %d\n"
- "Last ray read: %d\n", isweep+1, iray);
- wsr88d_load_sweep_header(radar, isweep, wsr88d_ray);
+ else
+ fprintf(stderr,"Read failed.\n");
+ fprintf(stderr,"Current sweep index: %d\n"
+ "Last ray read: %d\n", isweep, iray);
+ wsr88d_load_sweep_header(radar, isweep);
return radar;
}
if (msghdr.msg_type == 5) {
/* Check for end of sweep */
if (wsr88d_ray.ray_hdr.radial_status == END_OF_ELEV) {
- wsr88d_load_sweep_header(radar, isweep, wsr88d_ray);
+ wsr88d_load_sweep_header(radar, isweep);
isweep++;
iray = 0;
}
+ /* If not at end of volume scan, read next message header. */
if (wsr88d_ray.ray_hdr.radial_status != END_VOS) {
n = fread(&msghdr, sizeof(Wsr88d_msg_hdr), 1, wf->fptr);
if (n < 1) {
if (feof(wf->fptr) != 0) fprintf(stderr,
"Unexpected end of file.\n");
else fprintf(stderr,"Failed reading msghdr.\n");
- fprintf(stderr,"Current sweep number: %d\n"
- "Last ray read: %d\n", isweep+1, iray);
- wsr88d_load_sweep_header(radar, isweep, wsr88d_ray);
+ fprintf(stderr,"Current sweep index: %d\n"
+ "Last ray read: %d\n", isweep, iray);
+ wsr88d_load_sweep_header(radar, isweep);
return radar;
}
}
else {
end_of_vos = 1;
- wsr88d_load_sweep_header(radar, isweep, wsr88d_ray);
+ wsr88d_load_sweep_header(radar, isweep);
}
if (feof(wf->fptr) != 0) end_of_vos = 1;
- }
+ } /* while not end of vos */
return radar;
}
projects / ~andy / rsl / blobdiff