X-Git-Url: http://pileus.org/git/?p=~andy%2Frsl;a=blobdiff_plain;f=wsr88d_m31.c;h=4373df3220f8033fa034f4056581cd0a9bed7891;hp=3c260a711bef0c444b46a4afdb11ec79358d4694;hb=d08a7f8a699a044bc4ac5f93917aa7f6c463923b;hpb=0a59ff8412de3e114b7b5fd081a5fd39102542c1 diff --git a/wsr88d_m31.c b/wsr88d_m31.c index 3c260a7..4373df3 100644 --- a/wsr88d_m31.c +++ b/wsr88d_m31.c @@ -28,6 +28,7 @@ #include "rsl.h" #include "wsr88d.h" +#include /* Data descriptions in the following data structures are from the "Interface * Control Document for the RDA/RPG", Build 10.0 Draft, WSR-88D Radar @@ -35,17 +36,15 @@ */ 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 { @@ -91,32 +90,18 @@ 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) @@ -132,7 +117,7 @@ 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); @@ -141,13 +126,13 @@ void wsr88d_swap_m31_ray(Ray_header_m31 *wsr88d_ray) 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; @@ -158,25 +143,6 @@ void wsr88d_swap_data_moment(Data_moment_hdr *this_field) } -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; @@ -273,31 +239,39 @@ void wsr88d_get_vcp_data(short *msgtype5) } -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.; + } } @@ -305,8 +279,8 @@ int read_wsr88d_ray_m31(Wsr88d_file *wf, int msg_size, 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); @@ -314,170 +288,33 @@ int read_wsr88d_ray_m31(Wsr88d_file *wf, int msg_size, 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; @@ -501,118 +338,100 @@ 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 / 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; @@ -623,61 +442,69 @@ void wsr88d_load_ray_into_radar(Wsr88d_ray_m31 wsr88d_ray, int isweep, int iray, 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.; } @@ -694,41 +521,64 @@ Radar *load_wsr88d_m31_into_radar(Wsr88d_file *wf) 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, @@ -737,10 +587,11 @@ Radar *load_wsr88d_m31_into_radar(Wsr88d_file *wf) 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) { @@ -752,11 +603,12 @@ Radar *load_wsr88d_m31_into_radar(Wsr88d_file *wf) /* 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) { @@ -764,18 +616,18 @@ Radar *load_wsr88d_m31_into_radar(Wsr88d_file *wf) 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; }