From: Andy Spencer Date: Wed, 28 Oct 2009 01:55:27 +0000 (+0000) Subject: Changes from Bart (2009-10-28) X-Git-Tag: v1.41~4 X-Git-Url: http://pileus.org/git/?p=~andy%2Frsl;a=commitdiff_plain;h=e953e14aeb8c5b83435ee0103f3b069cc1875c96 Changes from Bart (2009-10-28) --- diff --git a/CHANGES b/CHANGES index af29ca9..9266df3 100644 --- a/CHANGES +++ b/CHANGES @@ -1,6 +1,20 @@ /* Changes for RSL * *--------------------------------------------------------------------- + * v1.41 In progress + * + * 1. wsr88d_m31.c: Simplified the WSR-88D ray structure and supporting code. + * Removed function read_data_moment. + * Added function get_wsr88d_unamb_and_nyq_vel. + * 2. Added support for Sigmet 2-byte data types, as well as HydroClass 1 and 2 byte + * types. + * Files involved: nsig_to_radar.c, nsig.c, nsig.h, volume.c. + * [TODO: dual PRF mods] + * 3. Completed DORADE ingest. Files involved: dorade_to_radar.c, dorade.c, volume.c. + * 4. prune.c: Added global variable prune_radar and two functions to set or unset: + * RSL_prune_radar_on() + * RSL_prune_radar_off() + *--------------------------------------------------------------------- * v1.40 Released 10/10/2008 * * 1. wsr88d_m31.c (load_wsr88d_m31_into_radar): Corrected a potential problem diff --git a/anyformat_to_radar.c b/anyformat_to_radar.c index 601a6f3..70be908 100644 --- a/anyformat_to_radar.c +++ b/anyformat_to_radar.c @@ -107,6 +107,7 @@ enum File_type RSL_filetype(char *infile) /* Byte swapped ? */ if (strncmp(" P A B", magic, 6) == 0) return MCGILL_FILE; if (strncmp("Volume", magic, 6) == 0) return EDGE_FILE; + if (strncmp("SSWB", magic, 4) == 0) return DORADE_FILE; if (strncmp("VOLD", magic, 4) == 0) return DORADE_FILE; return UNKNOWN; diff --git a/configure.in b/configure.in index 8697f3e..e0d59af 100644 --- a/configure.in +++ b/configure.in @@ -1,5 +1,5 @@ dnl Process this file with autoconf to produce a configure script. -AC_INIT(rsl, v1.40) +AC_INIT(rsl, v1.40.1) AC_CONFIG_SRCDIR(volume.c) AM_INIT_AUTOMAKE diff --git a/dorade.c b/dorade.c index 6c692bf..bf78a1c 100644 --- a/dorade.c +++ b/dorade.c @@ -38,6 +38,32 @@ void dorade_verbose_off() dorade_verbose = 0; } +static int do_swap = 0; + +/**********************************************************************/ +/* */ +/* read_extra_bytes */ +/* */ +/**********************************************************************/ +int read_extra_bytes(int nbytes, FILE *in) +{ + char *extra; + int nread; + /* + * Read and discard nbytes bytes. The return value is the byte count + * returned by fread, unless there is an error, then it is 0. + */ + extra = (char *) malloc(nbytes); + if (!extra) { + perror("\nError: read_extra_bytes"); + fprintf(stderr,"Tried to allocate %d bytes\n", nbytes); + return 0; + } + nread = fread(extra, sizeof(char), nbytes, in); + free(extra); + return nread; +} + /**********************************************************************/ /* */ /* dorade_read_comment_block */ @@ -52,13 +78,25 @@ Comment_block *dorade_read_comment_block(FILE *in) return NULL; } fread(cb->code, sizeof(cb->code), 1, in); - fread(&cb->len, sizeof(cb->len), 1, in); cb->len = ntohl(cb->len); - cb->comment = (char *) calloc(cb->len, sizeof(char)); + fread(&cb->len, sizeof(cb->len), 1, in); + + /* Check for big endian data on little endian platform. The smallest value + * cb->len could have is 8 (length of cb->code + cb->len), so we put that in + * first byte of test value, and also test for negative, since sign bit could + * be set in a larger value. + */ + if (cb->len > 0x08000000 || cb->len < 0) do_swap = 1; + if (do_swap) cb->len = ntohl(cb->len); + + /* Length of cb->comment is cb->len-8 since cb->code and cb->len have + * already been read. + */ + cb->comment = (char *) calloc(cb->len-8, sizeof(char)); if (cb->comment == NULL) { perror("dorade_read_comment_block: cb->comment"); return cb; } - fread(cb->comment, sizeof(char), cb->len, in); + fread(cb->comment, sizeof(char), cb->len-8, in); return cb; } @@ -79,20 +117,22 @@ Volume_desc *dorade_read_volume_desc (FILE *in) fread(vd, sizeof(Volume_desc), 1, in); /* Now, convert from Big Endian. */ - vd->len = ntohl(vd->len); - vd->version = ntohs(vd->version); - vd->volume_number = ntohs(vd->volume_number); - vd->max_bytes = ntohl(vd->max_bytes); - vd->year = ntohs(vd->year); - vd->month = ntohs(vd->month); - vd->day = ntohs(vd->day); - vd->hour = ntohs(vd->hour); - vd->minute = ntohs(vd->minute); - vd->second = ntohs(vd->second); - vd->gen_year = ntohs(vd->gen_year); - vd->gen_month = ntohs(vd->gen_month); - vd->gen_day = ntohs(vd->gen_day); - vd->nsensors = ntohs(vd->nsensors); + if (do_swap) { + vd->len = ntohl(vd->len); + vd->version = ntohs(vd->version); + vd->volume_number = ntohs(vd->volume_number); + vd->max_bytes = ntohl(vd->max_bytes); + vd->year = ntohs(vd->year); + vd->month = ntohs(vd->month); + vd->day = ntohs(vd->day); + vd->hour = ntohs(vd->hour); + vd->minute = ntohs(vd->minute); + vd->second = ntohs(vd->second); + vd->gen_year = ntohs(vd->gen_year); + vd->gen_month = ntohs(vd->gen_month); + vd->gen_day = ntohs(vd->gen_day); + vd->nsensors = ntohs(vd->nsensors); + } return vd; } @@ -120,7 +160,7 @@ Radar_desc *dorade_read_radar_desc (FILE *in) fread(rd, sizeof(Radar_desc), 1, in); /* Now, convert from Big Endian. */ - if (little_endian()) { + if (do_swap) { swap_4_bytes(&rd->len); swap_4_bytes(&rd->radar_constant); /* Yes, even the ieee floating values. */ swap_4_bytes(&rd->peak_power); @@ -153,6 +193,13 @@ Radar_desc *dorade_read_radar_desc (FILE *in) swap_4_bytes(&rd->period[i]); } } + /* If RADD block is longer than structure, read through extra bytes. + * This sometimes happens. + */ + if (rd->len > sizeof(Radar_desc)) { + if (read_extra_bytes(rd->len - sizeof(Radar_desc), in) <= 0) + fprintf(stderr,"Called from %s, line: %d\n", __FILE__, __LINE__ - 1); + } return rd; } /**********************************************************************/ @@ -172,7 +219,7 @@ Parameter_desc *dorade_read_parameter_desc (FILE *in) fread(pd, sizeof(Parameter_desc), 1, in); /* Now, convert from Big Endian. */ - if (little_endian()) { + if (do_swap) { swap_4_bytes(&pd->len); swap_2_bytes(&pd->ipp); swap_2_bytes(&pd->xmit_freq); @@ -186,6 +233,12 @@ Parameter_desc *dorade_read_parameter_desc (FILE *in) swap_4_bytes(&pd->offset_factor); swap_4_bytes(&pd->missing_data_flag); } + /* If the descriptor block is longer than the structure, read past the extra bytes. + */ + if (pd->len > sizeof(Parameter_desc)) { + if (read_extra_bytes(pd->len - sizeof(Parameter_desc), in) <= 0) + fprintf(stderr,"Called from %s, line: %d\n", __FILE__, __LINE__ - 1); + } return pd; } @@ -209,7 +262,7 @@ Cell_range_vector *dorade_read_cell_range_vector (FILE *in) fread(&cv->code, sizeof(cv->code), 1, in); fread(&cv->len, sizeof(cv->len), 1, in); fread(&cv->ncells, sizeof(cv->ncells), 1, in); - if (little_endian()) { + if (do_swap) { swap_4_bytes(&cv->len); swap_4_bytes(&cv->ncells); } @@ -220,7 +273,7 @@ Cell_range_vector *dorade_read_cell_range_vector (FILE *in) } fread(cv->range_cell, sizeof(float), cv->ncells, in); - if (little_endian()) { + if (do_swap) { for (i=0; incells; i++) swap_4_bytes(&cv->range_cell[i]); } @@ -249,6 +302,8 @@ Cell_range_vector *dorade_read_cell_range_vector (FILE *in) Correction_factor_desc *dorade_read_correction_factor_desc(FILE *in) { Correction_factor_desc *cf; + char *remaining; + int is_cfac = 0; cf = (Correction_factor_desc *) calloc(1, sizeof(Correction_factor_desc)); if(!cf) { @@ -256,9 +311,27 @@ Correction_factor_desc *dorade_read_correction_factor_desc(FILE *in) return NULL; } - fread(cf, sizeof(Correction_factor_desc), 1, in); + /* Make sure we have Correction Factor Descriptor. */ + while (!is_cfac) { + fread(cf->code, sizeof(cf->code), 1, in); + if (strncmp(cf->code, "CFAC", 4) == 0) + is_cfac = 1; + else { + fread(&cf->len, sizeof(cf->len), 1, in); + if (do_swap) swap_4_bytes(&cf->len); + remaining = (char *) malloc(cf->len-8); + if (!remaining) { + perror("\ndorade_read_correction_factor_desc"); + fprintf(stderr,"cf->len = %d\n\n", cf->len); + return NULL; + } + fread(remaining, sizeof(char), cf->len-8, in); + free(remaining); + } + } + fread(&cf->len, sizeof(Correction_factor_desc)-4, 1, in); /* Now, convert from Big Endian. */ - if (little_endian()) { + if (do_swap) { swap_4_bytes(&cf->len); swap_4_bytes(&cf->azimuth); swap_4_bytes(&cf->elevation); @@ -343,7 +416,7 @@ Sweep_info *dorade_read_sweep_info(FILE *in) } /* Now, convert from Big Endian. */ - if (little_endian()) { + if (do_swap) { swap_4_bytes(&si->len); swap_4_bytes(&si->sweep_num); swap_4_bytes(&si->nrays); @@ -375,7 +448,7 @@ Ray_info *dorade_read_ray_info (FILE *in) fread(ri, sizeof(Ray_info), 1, in); /* Now, convert from Big Endian. */ - if (little_endian()) { + if (do_swap) { swap_4_bytes(&ri->len); swap_4_bytes(&ri->sweep_num); swap_4_bytes(&ri->jday); @@ -401,6 +474,7 @@ Ray_info *dorade_read_ray_info (FILE *in) Platform_info *dorade_read_platform_info (FILE *in) { Platform_info *pi; + int len_first_two; pi = (Platform_info *) calloc(1, sizeof(Platform_info)); if(!pi) { @@ -408,35 +482,61 @@ Platform_info *dorade_read_platform_info (FILE *in) return NULL; } - fread(pi, sizeof(Platform_info), 1, in); - /* Now, convert from Big Endian. */ - if (little_endian()) { - swap_4_bytes(&pi->len); - swap_4_bytes(&pi->longitude); - swap_4_bytes(&pi->latitude); - swap_4_bytes(&pi->altitude); - swap_4_bytes(&pi->height); - swap_4_bytes(&pi->ew_speed); - swap_4_bytes(&pi->ns_speed); - swap_4_bytes(&pi->v_speed); - swap_4_bytes(&pi->heading); - swap_4_bytes(&pi->roll); - swap_4_bytes(&pi->pitch); - swap_4_bytes(&pi->drift); - swap_4_bytes(&pi->rotation); - swap_4_bytes(&pi->tilt); - swap_4_bytes(&pi->ew_wind_speed); - swap_4_bytes(&pi->ns_wind_speed); - swap_4_bytes(&pi->v_wind_speed); - swap_4_bytes(&pi->heading_rate); - swap_4_bytes(&pi->pitch_rate); + /* Read the id code to make sure we have "ASIB" for platform info. If + * id is ASIB, then read data into the Platform_info structure. If it is + * XSTF, read and discard remainder of block, which will have a different + * size than Platform_info. XSTF is undocumented, but apparently it takes + * the place of ASIB when radar is grounded. + */ + + fread(pi->code, sizeof(pi->code), 1, in); + fread(&pi->len, sizeof(pi->len), 1, in); + if (do_swap) swap_4_bytes(&pi->len); + len_first_two = sizeof(pi->code) + sizeof(pi->len); + + if (strncmp(pi->code, "ASIB", 4) == 0) { + fread(&pi->longitude, sizeof(Platform_info)-len_first_two, 1, in); + /* Read past any extra bytes. */ + if (pi->len > sizeof(Platform_info)) { + if (read_extra_bytes(pi->len - sizeof(Platform_info), in) <= 0) + fprintf(stderr,"Called from %s, line: %d\n",__FILE__,__LINE__-1); + } + /* Now, convert from Big Endian. */ + if (do_swap) { + swap_4_bytes(&pi->longitude); + swap_4_bytes(&pi->latitude); + swap_4_bytes(&pi->altitude); + swap_4_bytes(&pi->height); + swap_4_bytes(&pi->ew_speed); + swap_4_bytes(&pi->ns_speed); + swap_4_bytes(&pi->v_speed); + swap_4_bytes(&pi->heading); + swap_4_bytes(&pi->roll); + swap_4_bytes(&pi->pitch); + swap_4_bytes(&pi->drift); + swap_4_bytes(&pi->rotation); + swap_4_bytes(&pi->tilt); + swap_4_bytes(&pi->ew_wind_speed); + swap_4_bytes(&pi->ns_wind_speed); + swap_4_bytes(&pi->v_wind_speed); + swap_4_bytes(&pi->heading_rate); + swap_4_bytes(&pi->pitch_rate); + } + } else if (strncmp(pi->code, "XSTF", 4) == 0) { + /* Read to end of XSTF block. */ + if (read_extra_bytes(pi->len - len_first_two, in) <= 0) + fprintf(stderr,"Called from %s, line: %d\n", __FILE__, __LINE__ - 1); + } else { + fprintf(stderr,"Unexpected block id: \"%s\"." + " Expected \"ASIB\" or \"XSTF\"\n", pi->code); } return pi; } + /**********************************************************************/ /* */ -/* dorade_read_parameter_info */ +/* dorade_read_parameter_data */ /* */ /**********************************************************************/ @@ -454,7 +554,7 @@ Parameter_data *dorade_read_parameter_data(FILE *in) fread(&pd->code, sizeof(pd->code), 1, in); fread(&pd->len, sizeof(pd->len), 1, in); fread(&pd->name, sizeof(pd->name), 1, in); - if (little_endian()) swap_4_bytes(&pd->len); + if (do_swap) swap_4_bytes(&pd->len); /* Length is in parameter data block? or calculate if from pd->len. */ len = pd->len /* Use pd->len for now. */ @@ -553,7 +653,7 @@ Sweep_record *dorade_read_sweep(FILE *fp, Sensor_desc **sd) parameter_desc[j]->parameter_type == 4) sr->data_ray[i]->word_size[j] = 4; /* 4 bytes per word */ - if (little_endian()) { /* Numbers were read big-endian. */ + if (do_swap) { /* Numbers were read big-endian. */ if (sr->data_ray[i]->word_size[j] == 2) for (k=0; kdata[k]); @@ -566,13 +666,6 @@ Sweep_record *dorade_read_sweep(FILE *fp, Sensor_desc **sd) return sr; } -/**********************************************************************/ -/* */ -/* dorade_read_ray */ -/* */ -/**********************************************************************/ -Data_ray *dorade_read_ray (FILE *in); - /* MEMORY MANAGEMENT ROUTINES */ diff --git a/dorade.h b/dorade.h index 941a25f..2e87cc3 100644 --- a/dorade.h +++ b/dorade.h @@ -44,7 +44,7 @@ typedef struct { short minute; /* Minute 0-59 (2 byte int) */ short second; /* Second 0-59 (2 byte int) */ char flight_num[8]; /* Flight number (8 ASCII) for airborne raar or IOP number - for ground based radars. */ + for ground based radars. */ char facility_name[8]; /* Generation facility. (8 ASCII) */ short gen_year; /* Generation year YYYY. (2 byte int) */ short gen_month; /* Generation month 1-12. (2 byte int) */ @@ -65,61 +65,61 @@ typedef struct { float horizontal_beam_width; /* [deg] */ float vertical_beam_width; /* [deg] */ short radar_type; /* 0--Ground - 1--Airborne fore - 2--Airborne aft - 3--Airborne tail - 4--Airborne lower fuselage - 5--Shipborne - */ + 1--Airborne fore + 2--Airborne aft + 3--Airborne tail + 4--Airborne lower fuselage + 5--Shipborne + */ short scan_mode; /* 0--Calibration - 1--PPI (Constant elevation) - 2--Coplane - 3--RHI (Constant azimuth) - 4--Vertical pointing - 5--Target (Stationary, not vertical pointing) - 6--Manual - 7--Idle (out of control) - 8--Surveillance - 9--Vertical sweep (rotation axis parallels the fuselage) - */ + 1--PPI (Constant elevation) + 2--Coplane + 3--RHI (Constant azimuth) + 4--Vertical pointing + 5--Target (Stationary, not vertical pointing) + 6--Manual + 7--Idle (out of control) + 8--Surveillance + 9--Vertical sweep (rotation axis parallels the fuselage) + */ float scan_rate; /* Nominal scan rate. [deg/sec] */ float start_angle; /* Nominal start angle. [deg] */ float stop_angle; /* Nominal stop angle. [deg] */ short nparam_desc; /* Total number of parameter descriptors for this radar. (2 byte int) */ short ndesc; /* Total number of descriptors for this radar. (2 byte int) */ short compress_code; /* Data compression format code. (2 byte int): - 0--no compression - 1--data compression (compression algorithm described - in the ASCII file athe the begining of the file. - */ + 0--no compression + 1--data compression (compression algorithm described + in the ASCII file athe the begining of the file. + */ short compress_algo; /* Data reduction algorithm: - 0--No data reduction. - 1--Data recorded between two rotation angles. - 2--Data recorded between two concentic angles. - 3--Data recorded between two altitudes. - 4-N--Other types of data reduction. - */ + 0--No data reduction. + 1--Data recorded between two rotation angles. + 2--Data recorded between two concentic angles. + 3--Data recorded between two altitudes. + 4-N--Other types of data reduction. + */ float data_reduction_param1; /* Data reduction specific parameter #1 - 1--Smallest positive angle [deg]. - 2--Inner circle diameter [km]. - 3--Minimum altitude [km]. - 4-N--Will be defined if other types created. - */ + 1--Smallest positive angle [deg]. + 2--Inner circle diameter [km]. + 3--Minimum altitude [km]. + 4-N--Will be defined if other types created. + */ float data_reduction_param2; /* Data reduction specific parameter #2 - 1--Largest positive angle [deg]. - 2--Outer circle diameter [km]. - 3--Maximum altitude [km]. - 4-N--Will be defined if other types created. - */ + 1--Largest positive angle [deg]. + 2--Outer circle diameter [km]. + 3--Maximum altitude [km]. + 4-N--Will be defined if other types created. + */ float longitude; /* Radar longitude [deg]. If airborne, airport longitude. */ float latitude; /* Radar latitude [deg]. If airborne, airport latitude. */ float altitude; /* Radar altitude of mean sea level (msl) [km]. - If airborne, airport altitude. */ + If airborne, airport altitude. */ float unambiguous_velocity; /* Effective unambiguous velocity [m/s]. */ float unambiguous_range; /* Effective unambiguous range [km]. */ short nfreq; /* Number of freqencies transmitted (2 byte int). */ short npulse_periods; /* Number of different inter-pulse periods (IPP's) transmitted. - (2 byte int). */ + (2 byte int). */ float freq[5]; /* Frequency 1..5 [GHz] (float) */ float period[5]; /* Interpulse Period (IPP) 1..5 [ms] (float) */ } Radar_desc; @@ -131,32 +131,32 @@ typedef struct { char description[40]; /* Description of the parameter. */ char units[8]; /* Units (8 ASCII) */ short ipp; /* Inter-pulse periods. - Bit 0 set to 1 indicates IPP#1 is used in this parameter. - Similiarly for bits 1,2,3 and 4 and IPP's 2,3,4 and 5. - */ + Bit 0 set to 1 indicates IPP#1 is used in this parameter. + Similiarly for bits 1,2,3 and 4 and IPP's 2,3,4 and 5. + */ short xmit_freq; /* Transmittd frequencies. - Bit 0 set to 1 indicates Frequency#1 is used in this parameter. - Similiarly for bits 1,2,3 and 4 and Frequencies 2,3,4 and 5. - */ + Bit 0 set to 1 indicates Frequency#1 is used in this parameter. + Similiarly for bits 1,2,3 and 4 and Frequencies 2,3,4 and 5. + */ float rcvr_bandwidth;/* [MHz] */ short pulse_width; /* [m] */ short polarization; /* 0--Horizontal - 1--Vertical - 2--Circular, Right handed - 3--Elliptical - 4--Circular, Left handed - 5--Dual polarization - */ + 1--Vertical + 2--Circular, Right handed + 3--Elliptical + 4--Circular, Left handed + 5--Dual polarization + */ short nsamp_in_dwell_time; /* Number of samples in dwell time. */ short parameter_type; /* 1--8 bit integer - 2--16 bit integer - 3--32 bit integer - 4--floating point (32 bit IEEE) - */ + 2--16 bit integer + 3--32 bit integer + 4--floating point (32 bit IEEE) + */ char threshold_field[8]; float threshold_value; /* Units depend on the threshold field. */ float scale_factor; /* Scale factor. */ - float offset_factor; /* meteorological val = (recorded val - offset factor / scale factor */ + float offset_factor; /* meteorological val = (recorded val - offset factor) / scale factor */ int missing_data_flag; /* Deleted or missing data flag. 256for bytes, -999 for all others. */ } Parameter_desc; @@ -198,9 +198,9 @@ typedef struct { float stop_angle; /* True stop angle [deg]. */ float fixed_angle; /* Fixed angle [deg]. */ int filter_flag; /* Filter flag: - 0--No filtering in use. - 1--ON (Algorithm described in ASCII file at beginning of the file.) - */ + 0--No filtering in use. + 1--ON (Algorithm described in ASCII file at beginning of the file.) + */ } Sweep_info; typedef struct { @@ -214,14 +214,14 @@ typedef struct { short msec; /* Millisecond. */ float azimuth; /* [deg] */ float elevation; /* [deg] */ - float peak_power;/* [kw] */ - float scan_rate; /* Peak transmitted power [deg/sec]. */ + float peak_power;/* Peak transmitted power [kw] */ + float scan_rate; /* [deg/sec]. */ int status; /* Ray status: - 0--Normal - 1--Transition (antenna repositioning) - 2--Bad - 3--Questionable - */ + 0--Normal + 1--Transition (antenna repositioning) + 2--Bad + 3--Questionable + */ } Ray_info; typedef struct { /* Especially for moving radars. */ @@ -231,8 +231,8 @@ typedef struct { /* Especially for moving radars. */ float latitude; /* Radar latitude [deg]. */ float altitude; /* Radar pressure altitude (msl) [km]. */ float height; /* Radar above ground altitude (agl) [km]. */ - float ew_speed; /* Platform ground speed (E->W) [m/s]. */ - float ns_speed; /* Platform ground speed (N->S) [m/s]. */ + float ew_speed; /* Platform ground speed (E (positive) or W) [m/s]. */ + float ns_speed; /* Platform ground speed (N (positive) or S) [m/s]. */ float v_speed; /* Platform vertical velocity [m/s]. */ float heading; /* Platform heading [deg]. */ float roll; /* Platform roll [deg]. */ @@ -240,9 +240,9 @@ typedef struct { /* Especially for moving radars. */ float drift; /* Platform drift [deg]. */ float rotation; /* Platform rotation angle [deg]. */ float tilt; /* Platform tilt [deg]. */ - float ew_wind_speed; /* Horizontal wind speed at radar (East->West) [m/s]. */ - float ns_wind_speed; /* Horizontal wind speed at radar (North->South) [m/s]. */ - float v_wind_speed; /* Vertical wind speed at radar (East->West) [m/s]. */ + float ew_wind_speed; /* Horizontal wind speed at radar (toward East positive) [m/s]. */ + float ns_wind_speed; /* Horizontal wind speed at radar (toward North positive) [m/s]. */ + float v_wind_speed; /* Vertical wind speed at radar (up is positive) [m/s]. */ float heading_rate; /* Heading change rate [deg/sec]. */ float pitch_rate; /* Pitch change rate [deg/sec]. */ } Platform_info; @@ -268,12 +268,12 @@ typedef struct { Platform_info *platform_info; int nparam; int *data_len; /* 0..nparam-1 - * Length of *parameter_data[i] in words. - * Words can be 2 or 4 bytes. - */ + * Length of *parameter_data[i] in bytes. + * This is length of data portion. + */ int *word_size; /* 0..nparam-1 - * Size of each word in *parameter_data[i]. - */ + * Size of each word in *parameter_data[i]. + */ Parameter_data **parameter_data; /* 0..nparam-1 */ } Data_ray; diff --git a/dorade_print.c b/dorade_print.c index 85e28cf..8b20407 100644 --- a/dorade_print.c +++ b/dorade_print.c @@ -109,8 +109,10 @@ void dorade_print_cell_range_vector(Cell_range_vector *d) printf(">\n"); printf("len = %d\n", d->len); printf("d->ncells = %d\n", d->ncells); + /* for (i=0; incells; i++) printf("d->range_cell[%d] = %f\n", i, d->range_cell[i]); + */ } void dorade_print_parameter_desc(Parameter_desc *d) @@ -240,6 +242,6 @@ void dorade_print_sensor(Sensor_desc *s) dorade_print_parameter_desc(s->p_desc[i]); printf("=====================================================\n"); } - /* dorade_print_cell_range_vector(s->cell_range_vector); */ + dorade_print_cell_range_vector(s->cell_range_vector); dorade_print_correction_factor_desc(s->correction_factor_desc); } diff --git a/nsig.c b/nsig.c index 038b7de..7e5c3c1 100644 --- a/nsig.c +++ b/nsig.c @@ -69,11 +69,11 @@ FILE *nsig_open(char *file_name) /* Open input file */ if (file_name == NULL) { /* Use stdin */ - save_fd = dup(0); - fp = fdopen(save_fd, "r"); + save_fd = dup(0); + fp = fdopen(save_fd, "r"); } else if((fp = fopen(file_name,"r")) == NULL) { - perror(file_name); - return fp; + perror(file_name); + return fp; } fp = uncompress_pipe(fp); /* Transparently gunzip. */ @@ -98,8 +98,8 @@ int nsig_read_record(FILE *fp, char *nsig_rec) if (feof(fp)) return -1; nbytes = 0; - while((n = fread(&buf[nbytes], sizeof(char),NSIG_BLOCK-nbytes, fp)) > 0) { - nbytes += n; + while((n = fread(&buf[nbytes], sizeof(char), NSIG_BLOCK-nbytes, fp)) > 0) { + nbytes += n; } return nbytes; } @@ -110,7 +110,7 @@ int nsig_read_record(FILE *fp, char *nsig_rec) *********************************************************/ void nsig_close(FILE *fp) { - rsl_pclose(fp); + rsl_pclose(fp); } static int do_swap; @@ -125,15 +125,15 @@ int nsig_endianess(NSIG_Record1 *rec1) printf("id = %d %d\n", (int)rec1->struct_head.id[0], (int)rec1->struct_head.id[1]); */ if (rec1->struct_head.id[0] == 0) { /* Possible little-endian */ - if (rec1->struct_head.id[1] >= 20) - /* This is a big-endian file. Typically, version 2. */ - do_swap = little_endian(); - else - do_swap = big_endian(); + if (rec1->struct_head.id[1] >= 20) + /* This is a big-endian file. Typically, version 2. */ + do_swap = little_endian(); + else + do_swap = big_endian(); } else if ((rec1->struct_head.id[1] == 0)) { /* Possible big-endian */ - if (rec1->struct_head.id[0] <= 7) - /* This is a little-endian file. Version 1. */ - do_swap = big_endian(); + if (rec1->struct_head.id[0] <= 7) + /* This is a little-endian file. Version 1. */ + do_swap = big_endian(); } /* printf("DO SWAP = %d\n", do_swap); @@ -170,16 +170,16 @@ void nsig_free_sweep(NSIG_Sweep **s) int i=0,itype; if (s == NULL) return; for (itype=0; itypenparams; itype++) { - if (s[itype] == NULL) continue; - - if (s[itype]->idh.data_type == NSIG_DTB_EXH) - free(s[itype]->ray[i]); - else { - for (i=0; iidh.num_rays_act); i++) - nsig_free_ray(s[itype]->ray[i]); - } - free(s[itype]->ray); - free(s[itype]); + if (s[itype] == NULL) continue; + + if (s[itype]->idh.data_type == NSIG_DTB_EXH) + free(s[itype]->ray[i]); + else { + for (i=0; iidh.num_rays_act); i++) + nsig_free_ray(s[itype]->ray[i]); + } + free(s[itype]->ray); + free(s[itype]); } free(s); } @@ -207,102 +207,102 @@ int nsig_read_chunk(FILE *fp, char *chunk) #define Vprint #undef Vprint while(the_code != 1) { - if (feof(fp)) return -1; - if (ipos == sizeof(data)) { /* the_code is in the next chunk */ + if (feof(fp)) return -1; + if (ipos == sizeof(data)) { /* the_code is in the next chunk */ #ifdef Vprint - printf("Exceeded block size looking for the_code. Get it from next buffer.\n"); + printf("Exceeded block size looking for the_code. Get it from next buffer.\n"); #endif - n = nsig_read_record(fp, (char *)data); - if (n <= 0) return n; /* Problem. */ + n = nsig_read_record(fp, (char *)data); + if (n <= 0) return n; /* Problem. */ #ifdef Vprint - printf("Read %d bytes.\n", n); - printf("Resetting ipos\n"); + printf("Read %d bytes.\n", n); + printf("Resetting ipos\n"); #endif - ipos = sizeof(NSIG_Raw_prod_bhdr); - } - + ipos = sizeof(NSIG_Raw_prod_bhdr); + } + #ifdef Vprint - printf("ipos = %d -- ", ipos); + printf("ipos = %d -- ", ipos); #endif - the_code = NSIG_I2(&data[ipos]); + the_code = NSIG_I2(&data[ipos]); #ifdef Vprint - printf("the_code = %d (%d) -- ", the_code, (unsigned short)the_code); + printf("the_code = %d (%d) -- ", the_code, (unsigned short)the_code); #endif - ipos += sizeof(twob); - - if (the_code < 0) { /* THIS IS DATA */ - nwords = the_code & 0x7fff; + ipos += sizeof(twob); + + if (the_code < 0) { /* THIS IS DATA */ + nwords = the_code & 0x7fff; #ifdef Vprint - printf("#data words (2-bytes) is %d\n", nwords); + printf("#data words (2-bytes) is %d\n", nwords); #endif - if (ipos + sizeof(twob)*nwords > sizeof(data)) { + if (ipos + sizeof(twob)*nwords > sizeof(data)) { #ifdef Vprint - printf("Exceeded block size... transferring and reading new (i=%d).\n", i); + printf("Exceeded block size... transferring and reading new (i=%d).\n", i); #endif - /* Need another phyical block. */ - - /* But, first transfer the remainder to the output chunk. */ - /* And, transfer begining of next block */ - /* Transfer end of current buffer. */ - end_nwords = (NSIG_BLOCK - ipos)/sizeof(twob); - memmove(&chunk[i], &data[ipos], sizeof(twob)*end_nwords); - i += end_nwords * sizeof(twob); - - n = nsig_read_record(fp, (char *)data); - if (n <= 0) return n; /* Problem. */ - /* New ipos */ - nwords -= end_nwords; - ipos = sizeof(NSIG_Raw_prod_bhdr); - - /* Transfer beginning of new buffer */ - if (i+nwords * sizeof(twob) > NSIG_BLOCK) return -1; - memmove(&chunk[i], &data[ipos], sizeof(twob) * nwords); - i += nwords * sizeof(twob); + /* Need another phyical block. */ + + /* But, first transfer the remainder to the output chunk. */ + /* And, transfer begining of next block */ + /* Transfer end of current buffer. */ + end_nwords = (NSIG_BLOCK - ipos)/sizeof(twob); + memmove(&chunk[i], &data[ipos], sizeof(twob)*end_nwords); + i += end_nwords * sizeof(twob); + + n = nsig_read_record(fp, (char *)data); + if (n <= 0) return n; /* Problem. */ + /* New ipos */ + nwords -= end_nwords; + ipos = sizeof(NSIG_Raw_prod_bhdr); + + /* Transfer beginning of new buffer */ + if (i+nwords * sizeof(twob) > NSIG_BLOCK) return -1; + memmove(&chunk[i], &data[ipos], sizeof(twob) * nwords); + i += nwords * sizeof(twob); ipos += nwords * sizeof(twob); - + #ifdef Vprint - printf("Words to transfer (at end of block) is %d\n", end_nwords); - printf("Transfer %d words from beginning of next buffer.\n", nwords); - printf("ipos in new buffer is %d\n", ipos); + printf("Words to transfer (at end of block) is %d\n", end_nwords); + printf("Transfer %d words from beginning of next buffer.\n", nwords); + printf("ipos in new buffer is %d\n", ipos); #endif - } else { /* Normal situation. Position to end of data. - * But, first transfer it to the chunk. - */ - if (i+nwords * sizeof(twob) > NSIG_BLOCK) return -1; - memmove(&chunk[i], &data[ipos], sizeof(twob) * nwords); - i += nwords * sizeof(twob); - ipos += sizeof(twob) * nwords; - } - - } else if (the_code == 1) { /* END OF THE RAY. */ + } else { /* Normal situation. Position to end of data. + * But, first transfer it to the chunk. + */ + if (i+nwords * sizeof(twob) > NSIG_BLOCK) return -1; + memmove(&chunk[i], &data[ipos], sizeof(twob) * nwords); + i += nwords * sizeof(twob); + ipos += sizeof(twob) * nwords; + } + + } else if (the_code == 1) { /* END OF THE RAY. */ #ifdef Vprint - printf("------------------------------> Reached end of ray.\n"); + printf("------------------------------> Reached end of ray.\n"); #endif - break; /* or continue; */ + break; /* or continue; */ - } else if (the_code == 0) { /* UNKNOWN */ - break; - } else { /* NUMBER OF ZERO's */ + } else if (the_code == 0) { /* UNKNOWN */ + break; + } else { /* NUMBER OF ZERO's */ #ifdef Vprint - printf("#000000000000 to skip is %d (i=%d)\n", the_code, i); + printf("#000000000000 to skip is %d (i=%d)\n", the_code, i); #endif - if (i+the_code * sizeof(twob) > NSIG_BLOCK) return -1; - memset(&chunk[i], 0, the_code*sizeof(twob)); - i += the_code * sizeof(twob); - } - - if (ipos >= sizeof(data)) { + if (i+the_code * sizeof(twob) > NSIG_BLOCK) return -1; + memset(&chunk[i], 0, the_code*sizeof(twob)); + i += the_code * sizeof(twob); + } + + if (ipos >= sizeof(data)) { #ifdef Vprint - printf("Exceeded block size ... ipos = %d\n", ipos); - printf("This should be right at the end of the block.\n"); + printf("Exceeded block size ... ipos = %d\n", ipos); + printf("This should be right at the end of the block.\n"); #endif - n = nsig_read_record(fp, (char *)data); - if (n <= 0) return n; /* Problem. */ - ipos = sizeof(NSIG_Raw_prod_bhdr); - } - - } /* End while. */ + n = nsig_read_record(fp, (char *)data); + if (n <= 0) return n; /* Problem. */ + ipos = sizeof(NSIG_Raw_prod_bhdr); + } + + } /* End while. */ return i; } @@ -354,14 +354,14 @@ NSIG_Ray *nsig_read_ray(FILE *fp) if (n == 0) return NULL; /* Silent error. */ if (n < 0) { - fprintf(stderr, "nsig_read_ray: chunk return code = %d.\n", n); - return NULL; + fprintf(stderr, "nsig_read_ray: chunk return code = %d.\n", n); + return NULL; } if (n > NSIG_BLOCK) { /* Whoa! */ - fprintf(stderr, "nsig_read_ray: chunk bigger than buffer. n = %d,\ + fprintf(stderr, "nsig_read_ray: chunk bigger than buffer. n = %d,\ maximum block size allowed is %d\n", n, NSIG_BLOCK); - return NULL; + return NULL; } ray = (NSIG_Ray *) calloc(1, sizeof(NSIG_Ray)); @@ -376,8 +376,13 @@ NSIG_Ray *nsig_read_ray(FILE *fp) #ifdef Vprint printf(" rayh.num_bins = %d (nbins %d, n %d)\n", NSIG_I2(rayh.num_bins), nbins, n); #endif + ray->range = (unsigned char *)calloc(n, sizeof(unsigned char)); + /* Changed calloc nbins to calloc n for 2-byte data. ray->range = (unsigned char *)calloc(nbins, sizeof(unsigned char)); memmove(ray->range, &chunk[sizeof(NSIG_Ray_header)], nbins); + Can remove this commented-out code once we know changes work. + */ + memmove(ray->range, &chunk[sizeof(NSIG_Ray_header)], n); return ray; } @@ -386,11 +391,12 @@ NSIG_Ray *nsig_read_ray(FILE *fp) NSIG_Sweep **nsig_read_sweep(FILE *fp, NSIG_Product_file *prod_file) { NSIG_Sweep **s; - int i, n; + int i, j, n; static NSIG_Ingest_data_header **idh = NULL; static NSIG_Raw_prod_bhdr *bhdr = NULL; NSIG_Ray *nsig_ray; int data_mask, iray, nrays[12], max_rays; + int masks[5]; int nparams; int is_new_ray; int idtype[12]; @@ -441,12 +447,25 @@ NSIG_Sweep **nsig_read_sweep(FILE *fp, NSIG_Product_file *prod_file) (void)nsig_endianess(&prod_file->rec1); /* Setup the array of ingest data headers [0..nparams-1] */ +#ifdef NSIG_VER2 + memmove(&masks[0], prod_file->rec2.task_config.dsp_info.data_mask_cur.mask_word_0, + sizeof(fourb)); + memmove(&masks[1], &prod_file->rec2.task_config.dsp_info.data_mask_cur.mask_word_1, + 4*sizeof(fourb)); + nparams = 0; + for (j=0; j < 5; j++) { + data_mask = masks[j]; + for (i=0; i<32; i++) + nparams += (data_mask >> i) & 0x1; + } +#else memmove(&data_mask, prod_file->rec2.task_config.dsp_info.data_mask, sizeof(fourb)); + for (nparams=i=0; i<32; i++) + nparams += (data_mask >> i) & 0x1; #ifdef Vprint printf("data_mask %x\n", data_mask); #endif - for (nparams=i=0; i<32; i++) - nparams += (data_mask >> i) & 0x1; +#endif /* Number of sweeps */ #ifdef Vprint {int nsweeps; @@ -457,13 +476,13 @@ NSIG_Sweep **nsig_read_sweep(FILE *fp, NSIG_Product_file *prod_file) if (idh == NULL) { - - idh = (NSIG_Ingest_data_header **)calloc(nparams, sizeof(NSIG_Ingest_data_header *)); - ipos = 0; - for (i=0; idata; + + idh = (NSIG_Ingest_data_header **)calloc(nparams, sizeof(NSIG_Ingest_data_header *)); + ipos = 0; + for (i=0; idata; } xh_size = NSIG_I2(prod_file->rec2.ingest_head.size_ext_ray_headers); @@ -503,32 +522,32 @@ NSIG_Sweep **nsig_read_sweep(FILE *fp, NSIG_Product_file *prod_file) /* Now pointers to all possible rays. */ for (i=0; inparams = nparams; - memmove(&s[i]->bhdr, &bhdr, sizeof(NSIG_Raw_prod_bhdr)); - memmove(&s[i]->idh, idh[i], sizeof(NSIG_Ingest_data_header)); - s[i]->ray = (NSIG_Ray **) calloc (max_rays, sizeof(NSIG_Ray *)); - } + s[i] = (NSIG_Sweep *) calloc (nparams, sizeof(NSIG_Sweep)); + s[i]->nparams = nparams; + memmove(&s[i]->bhdr, &bhdr, sizeof(NSIG_Raw_prod_bhdr)); + memmove(&s[i]->idh, idh[i], sizeof(NSIG_Ingest_data_header)); + s[i]->ray = (NSIG_Ray **) calloc (max_rays, sizeof(NSIG_Ray *)); + } /* Process this sweep. Keep track of the end of the ray. */ ipos = sizeof(NSIG_Raw_prod_bhdr); /* Position in the 'data' array */ max_rays = 0; for (i=0; idata_type); - nrays[i] = (int)NSIG_I2(idh[i]->num_rays_act); - if (nrays[i] > max_rays) max_rays = nrays[i]; + idtype[i] = NSIG_I2(idh[i]->data_type); + nrays[i] = (int)NSIG_I2(idh[i]->num_rays_act); + if (nrays[i] > max_rays) max_rays = nrays[i]; #ifdef Vprint - printf("New ray: parameter %d has idtype=%d\n", i, idtype[i]); - printf("Number of expected rays in sweep %d is %d\n", isweep, (int)NSIG_I2(idh[i]->num_rays_exp)); - printf("Number of actual rays in sweep %d is %d\n", isweep, (int)NSIG_I2(idh[i]->num_rays_act)); + printf("New ray: parameter %d has idtype=%d\n", i, idtype[i]); + printf("Number of expected rays in sweep %d is %d\n", isweep, (int)NSIG_I2(idh[i]->num_rays_exp)); + printf("Number of actual rays in sweep %d is %d\n", isweep, (int)NSIG_I2(idh[i]->num_rays_act)); #endif } if (is_new_sweep) - ipos += nparams * sizeof(NSIG_Ingest_data_header); + ipos += nparams * sizeof(NSIG_Ingest_data_header); - /* ipos = sizeof(NSIG_Raw_prod_bhdr) + nparams*sizeof(NSIG_Ingest_data_header); */ + /* ipos = sizeof(NSIG_Raw_prod_bhdr) + nparams*sizeof(NSIG_Ingest_data_header); */ /* 'iray' is the true ray index into 's', whereas, 'nsig_iray' is what * the NSIG file says it is. I'll trust 'iray' * @@ -564,44 +583,44 @@ NSIG_Sweep **nsig_read_sweep(FILE *fp, NSIG_Product_file *prod_file) do { #ifdef Vprint - printf("---------------------- New Ray <%d> --------------------\n", iray); + printf("---------------------- New Ray <%d> --------------------\n", iray); #endif - if (feof(fp)) { /* Premature eof */ - return NULL; /* This will have to do. */ - } - /* For all parameters present. */ - is_new_ray = 0; - for (i=0; irange = (unsigned char *)exh0; - } - } else { - exh1 = nsig_read_ext_header_ver1(fp); - if (exh1) { - nsig_ray = (NSIG_Ray *)calloc(1, sizeof(NSIG_Ray)); - nsig_ray->range = (unsigned char *)exh1; - } - } - } - if (nsig_ray) is_new_ray = 1; - if (iray > nrays[i]) break; - s[i]->ray[iray] = nsig_ray; - - } /* End for */ - if (is_new_ray) iray++; - + */ + nsig_ray = NULL; + if (idtype[i] != 0) { /* Not an extended header. */ + nsig_ray = nsig_read_ray(fp); + + } else { /* Check extended header version. */ + if (xh_size <= 20) { + exh0 = nsig_read_ext_header_ver0(fp); + if (exh0) { + nsig_ray = (NSIG_Ray *)calloc(1, sizeof(NSIG_Ray)); + nsig_ray->range = (unsigned char *)exh0; + } + } else { + exh1 = nsig_read_ext_header_ver1(fp); + if (exh1) { + nsig_ray = (NSIG_Ray *)calloc(1, sizeof(NSIG_Ray)); + nsig_ray->range = (unsigned char *)exh1; + } + } + } + if (nsig_ray) is_new_ray = 1; + if (iray > nrays[i]) break; + s[i]->ray[iray] = nsig_ray; + + } /* End for */ + if (is_new_ray) iray++; + } while (iray < max_rays); #ifdef Vprint printf("iray = %d\n", iray); diff --git a/nsig.h b/nsig.h index 63e44a5..5844ddf 100644 --- a/nsig.h +++ b/nsig.h @@ -22,10 +22,23 @@ #define NSIG_DTB_VEL 3 #define NSIG_DTB_WID 4 #define NSIG_DTB_ZDR 5 +#define NSIG_DTB_UCR2 8 +#define NSIG_DTB_CR2 9 +#define NSIG_DTB_VEL2 10 +#define NSIG_DTB_WID2 11 +#define NSIG_DTB_ZDR2 12 #define NSIG_DTB_KDP 14 +#define NSIG_DTB_KDP2 15 #define NSIG_DTB_PHIDP 16 +#define NSIG_DTB_VELC 17 #define NSIG_DTB_SQI 18 #define NSIG_DTB_RHOHV 19 +#define NSIG_DTB_RHOHV2 20 +#define NSIG_DTB_VELC2 22 +#define NSIG_DTB_SQI2 23 +#define NSIG_DTB_PHIDP2 24 +#define NSIG_DTB_HCLASS 55 +#define NSIG_DTB_HCLASS2 56 /* Product type code ,value for byte 12 in product configuration * struct, III-35 @@ -376,7 +389,7 @@ typedef struct { * Used for warning, shearline, and track only */ twob iresults_count; -#define PROD_END_PAD (20) +#define PROD_END_PAD 20 char ipad_end[PROD_END_PAD]; #else char color_unscan; /* Color used for unscanned area. */ @@ -597,18 +610,50 @@ typedef struct { /*============================================================*/ /*============================================================*/ +#ifdef NSIG_VER2 +typedef struct { + twob low_prf; /* Hertz */ + twob low_prf_frac; /* Fraction part, scaled by 2**-16 */ + twob low_prf_sample_size; + twob low_prf_range_averaging; /* In bins */ + twob thresh_refl_unfolding; /* Threshold for reflectivity unfolding in 1/100 dB */ + twob thresh_vel_unfolding; /* Threshold for velocity unfolding in 1/100 dB */ + twob thresh_sw_unfolding; /* Threshold for width unfolding in 1/100 dB */ + char spare[18]; +} NSIG_Task_dsp_mode_batch; +#endif +/*============================================================*/ +/*============================================================*/ + +#ifdef NSIG_VER2 +typedef struct { + fourb mask_word_0; + fourb ext_hdr_type; + fourb mask_word_1; + fourb mask_word_2; + fourb mask_word_3; + fourb mask_word_4; +} NSIG_Dsp_data_mask; +#endif +/*============================================================*/ +/*============================================================*/ + /* Task dsp info 3.3.2.2, page III-22 */ typedef struct { #ifdef NSIG_VER2 twob dsp_num; fourb dsp_type; + NSIG_Dsp_data_mask data_mask_cur; + NSIG_Dsp_data_mask data_mask_orig; + NSIG_Task_dsp_mode_batch task_dsp_mode; + char spare[52]; #else twob dsp_num; twob dsp_type; -#endif fourb data_mask; fourb aux_data_def[32]; +#endif fourb prf; fourb pwid; twob prf_mode; @@ -625,8 +670,14 @@ typedef struct { #endif twob atten_gain; #ifdef NSIG_VER2 - twob igas_atten ; /* 100000 * db/km */ -#define TASK_DSP_INFO_PAD (148) + twob igas_atten; /* 100000 * db/km */ + twob clutter_map_flag; + twob xmt_phase_seq; + fourb ray_hdr_mask; + twob time_series_playback; + twob spare2; + char custom_ray_hdr_name[16]; +#define TASK_DSP_INFO_PAD 120 char ipad_end[TASK_DSP_INFO_PAD]; #else char spare[150]; diff --git a/nsig_to_radar.c b/nsig_to_radar.c index 1e8b390..8454ab1 100644 --- a/nsig_to_radar.c +++ b/nsig_to_radar.c @@ -66,6 +66,7 @@ extern int rsl_qfield[]; /* See RSL_select_fields */ #define MAX_NSIG_SWEEPS 30 #define MAX_NSIG_RAYS 400 #define NSIG_NO_ECHO -32.0 +#define NSIG_NO_ECHO2 -999.0 static float (*f)(Range x); static Range (*invf)(float x); @@ -73,13 +74,14 @@ static Range (*invf)(float x); FILE *file; void get_extended_header_info(NSIG_Sweep **nsig_sweep, int xh_size, int iray, - int nparams, - int *msec, float *azm, float *elev, - float *pitch, float *roll, float *heading, - float *azm_rate, float *elev_rate, - float *pitch_rate, float *roll_rate, float *heading_rate, - float *lat, float *lon, int *alt, float *rvc, - float *vel_east, float *vel_north, float *vel_up) + int nparams, + int *msec, float *azm, float *elev, + float *pitch, float *roll, float *heading, + float *azm_rate, float *elev_rate, + float *pitch_rate, float *roll_rate, + float *heading_rate, + float *lat, float *lon, int *alt, float *rvc, + float *vel_east, float *vel_north, float *vel_up) { static NSIG_Ext_header_ver1 xh; int data_type, itype; @@ -91,7 +93,7 @@ void get_extended_header_info(NSIG_Sweep **nsig_sweep, int xh_size, int iray, /* Determine where 'itype' for extended header is. */ for (itype = 0; itypeidh.data_type); - if (data_type == NSIG_DTB_EXH) break; + if (data_type == NSIG_DTB_EXH) break; } /* printf("...extended header itype=%d, nparams=%d\n", itype, nparams); */ if (itype == nparams) return; /* No extended header. */ @@ -103,7 +105,7 @@ void get_extended_header_info(NSIG_Sweep **nsig_sweep, int xh_size, int iray, *msec = NSIG_I4(xh.msec); /* printf("...extended header msec= %d\n", *msec); */ if (xh_size <= 20) /* Stop, only have version 0. */ - return; + return; /* Version 1 processing. */ *azm = nsig_from_bang(xh.azm); @@ -141,62 +143,65 @@ RSL_nsig_to_radar #endif (char *filename) { - FILE *fp; - /* RSL structures */ - Radar *radar; - Ray *ray; - - int i, j, k, n; - int year, month, day; - int hour, minute, sec; - int numbins, numsweep; - int num_rays, sea_lvl_hgt; - int radar_number, num_samples; - int latd, latm, lats, lond, lonm, lons; - int data_type; - int bin_num; - int sweep_year, sweep_day, sweep_month; - int sweep_hour, sweep_minute, sweep_second; - int sweep_sec; - int z_flag_unc, z_flag_cor, v_flag, w_flag, speckle; - int ant_scan_mode; - float second; - float pw; - float bin_space; - float prf, wave, beam_width; - float vert_half_bw, horz_half_bw; - float rng_last_bin; - float rng_first_bin, freq; - float max_vel, sweep_rate, azim_rate; - float ray_data; - float az1, az2; - double tmp; - float sqi, log, csr, sig, cal_dbz; - char radar_type[50], state[2], city[15]; - char site_name[16]; + FILE *fp; + /* RSL structures */ + Radar *radar; + Ray *ray; + + int i, j, k, n; + int year, month, day; + int hour, minute, sec; + int numbins, numsweep; + int num_rays, sea_lvl_hgt; + int radar_number, num_samples; + int latd, latm, lats, lond, lonm, lons; + int data_type; + int bin_num; + int sweep_year, sweep_day, sweep_month; + int sweep_hour, sweep_minute, sweep_second; + int sweep_sec; + int z_flag_unc, z_flag_cor, v_flag, w_flag, speckle; + int ant_scan_mode; + float second; + float pw; + float bin_space; + float prf, wave, beam_width; + float vert_half_bw, horz_half_bw; + float rng_last_bin; + float rng_first_bin, freq; + float max_vel, sweep_rate, azim_rate; + float ray_data; + float az1, az2; + double tmp; + float sqi, log, csr, sig, cal_dbz; + char radar_type[50], state[2], city[15]; + char site_name[16]; NSIG_Product_file *prod_file; short id; int data_mask, nrays; + int masks[5]; int nparams, nsweeps; NSIG_Sweep **nsig_sweep; NSIG_Ray *ray_p; int itype, ifield; + unsigned short nsig_u2byte; /* New for 2-byte data types, Aug 2009 */ Sweep *sweep; int msec; float azm, elev, pitch, roll, heading, azm_rate, elev_rate, - pitch_rate, roll_rate, heading_rate, - lat, lon; + pitch_rate, roll_rate, heading_rate, + lat, lon; int alt; /* Altitude */ float rvc; /* Radial correction velocity m/s */ float vel_east, vel_north, vel_up; /* Platform velocity vectors m/sec */ int xh_size; + float incr; extern int *rsl_qsweep; /* See RSL_read_these_sweeps in volume.c */ extern int rsl_qsweep_max; extern float rsl_kdp_wavelen; radar = NULL; if (radar_verbose_flag) - fprintf(stderr, "open file: %s\n", filename); + fprintf(stderr, "open file: %s\n", filename); /** Opening nsig file **/ if((fp = nsig_open(filename)) == NULL) return NULL; @@ -218,30 +223,30 @@ RSL_nsig_to_radar n = nsig_read_record(fp, (char *)&prod_file->rec1); nsig_endianess(&prod_file->rec1); if (radar_verbose_flag) - fprintf(stderr, "Read %d bytes for rec1.\n", n); + fprintf(stderr, "Read %d bytes for rec1.\n", n); id = NSIG_I2(prod_file->rec1.struct_head.id); if (radar_verbose_flag) - fprintf(stderr, "ID = %d\n", (int)id); + fprintf(stderr, "ID = %d\n", (int)id); if (id != 7 && id != 27) { /* testing: Use 27 for Version 2 data */ - fprintf(stderr, "File is not a SIGMET version 1 nor version 2 raw product file.\n"); - return NULL; + fprintf(stderr, "File is not a SIGMET version 1 nor version 2 raw product file.\n"); + return NULL; } n = nsig_read_record(fp, (char *)&prod_file->rec2); if (radar_verbose_flag) - fprintf(stderr, "Read %d bytes for rec2.\n", n); + fprintf(stderr, "Read %d bytes for rec2.\n", n); /** Test for scan mode -- If scan is a RHI will return NULL **/ /** because RSL can't handle RHI's. In the future, replace **/ /** NULL will a routine to convert RHI's to RSL Format **/ ant_scan_mode =NSIG_I2(prod_file->rec2.task_config.scan_info.ant_scan_mode); if(ant_scan_mode == 2) - { - if (radar_verbose_flag) - fprintf(stderr, "RHI scan detected. Unable to process, returning NULL.\n"); - /* return NULL; */ - } + { + if (radar_verbose_flag) + fprintf(stderr, "RHI scan detected. Unable to process, returning NULL.\n"); + /* return NULL; */ + } /* Count the bits set in 'data_mask' to determine the number * of parameters present. @@ -249,10 +254,23 @@ RSL_nsig_to_radar xh_size = NSIG_I2(prod_file->rec2.ingest_head.size_ext_ray_headers); nrays = NSIG_I2(prod_file->rec2.ingest_head.num_rays); if (radar_verbose_flag) - fprintf(stderr, "Expecting %d rays in each sweep.\n", nrays); + fprintf(stderr, "Expecting %d rays in each sweep.\n", nrays); +#ifdef NSIG_VER2 + memmove(&masks[0], prod_file->rec2.task_config.dsp_info.data_mask_cur.mask_word_0, + sizeof(fourb)); + memmove(&masks[1], &prod_file->rec2.task_config.dsp_info.data_mask_cur.mask_word_1, + 4*sizeof(fourb)); + nparams = 0; + for (j=0; j < 5; j++) { + data_mask = masks[j]; + for (i=0; i<32; i++) + nparams += (data_mask >> i) & 0x1; + } +#else memmove(&data_mask, prod_file->rec2.task_config.dsp_info.data_mask, sizeof(fourb)); for (nparams=i=0; i<32; i++) - nparams += (data_mask >> i) & 0x1; + nparams += (data_mask >> i) & 0x1; +#endif /* Number of sweeps */ nsweeps = NSIG_I2(prod_file->rec2.task_config.scan_info.num_swp); @@ -262,8 +280,8 @@ RSL_nsig_to_radar memmove(site_name, prod_file->rec1.prod_end.site_name, sizeof(prod_file->rec1.prod_end.site_name)); site_name[sizeof(site_name)-1] = '\0'; if (radar_verbose_flag) { - fprintf(stderr, "nparams = %d, nsweeps = %d\n", nparams, nsweeps); - fprintf(stderr, "Site name = <%s>\n", site_name); + fprintf(stderr, "nparams = %d, nsweeps = %d\n", nparams, nsweeps); + fprintf(stderr, "Site name = <%s>\n", site_name); } /* nsig_sweep = nsig_read_sweep(fp, prod_file) @@ -319,22 +337,22 @@ RSL_nsig_to_radar sea_lvl_hgt = NSIG_I2(prod_file->rec1.prod_end.grnd_sea_ht); if (radar_verbose_flag) - fprintf(stderr, "sea: %d\n", sea_lvl_hgt); + fprintf(stderr, "sea: %d\n", sea_lvl_hgt); if (radar_verbose_flag) - fprintf(stderr, "site_name: %s", site_name); + fprintf(stderr, "site_name: %s", site_name); /** Determine beamwidth from input variables (not saved in nsig file) **/ if(strncmp(site_name,"mit",3) == 0 || strncmp(site_name,"MIT",3) == 0) - beam_width = MIT_BEAMWIDTH; + beam_width = MIT_BEAMWIDTH; else if(strncmp(site_name,"tog",3) == 0 || strncmp(site_name,"TOG",3) == 0) - beam_width = TOG_BEAMWIDTH; + beam_width = TOG_BEAMWIDTH; else if(strncmp(site_name,"kwa",3) == 0 || strncmp(site_name,"KWA",3) == 0) - beam_width = KWA_BEAMWIDTH; + beam_width = KWA_BEAMWIDTH; else - beam_width = DEFAULT_BEAMWIDTH; + beam_width = DEFAULT_BEAMWIDTH; if (radar_verbose_flag) - fprintf(stderr, "beamwidth: %f\n", beam_width); + fprintf(stderr, "beamwidth: %f\n", beam_width); vert_half_bw = beam_width/2.0; horz_half_bw = beam_width/2.0; @@ -358,8 +376,8 @@ RSL_nsig_to_radar prf = NSIG_I4(prod_file->rec1.prod_end.prf); /* pulse repetition frequency */ wave = (NSIG_I4(prod_file->rec1.prod_end.wavelen))/100.0; /* wavelength (cm) */ rsl_kdp_wavelen = wave; /* EXTERNAL (volume.c) This sets KD_F and KD_INVF - * to operate with the proper wavelength. - */ + * to operate with the proper wavelength. + */ numbins = NSIG_I4(prod_file->rec1.prod_end.num_bin); /* # bins in ray */ rng_first_bin = (float)NSIG_I4(prod_file->rec1.prod_end.rng_f_bin)/100.0; rng_last_bin = (float)NSIG_I4(prod_file->rec1.prod_end.rng_l_bin)/100.0; @@ -385,25 +403,25 @@ RSL_nsig_to_radar /** Verbose calibration information **/ if (radar_verbose_flag) - { - fprintf(stderr, "LOG = %5.2f\n", log); - fprintf(stderr, "SQI = %5.2f\n", sqi); - fprintf(stderr, "CSR = %5.2f\n", csr); - fprintf(stderr, "SIG = %5.2f\n", sig); - fprintf(stderr, "Calibration reflectivity: %5.2f dBZ\n", cal_dbz); - fprintf(stderr, "ZT flags: %d\n", z_flag_unc); /** can find these **/ - fprintf(stderr, "DZ flags: %d\n", z_flag_cor); /** defn in the **/ - fprintf(stderr, "VR flags: %d\n", v_flag); /** SIGMET Doc **/ - fprintf(stderr, "SW flags: %d\n", w_flag); - fprintf(stderr, "Flags: -3856 = SQI thresholding\n"); - fprintf(stderr, " -21846 = LOG thresholding\n"); - fprintf(stderr, " -24416 = LOG & SQI thresholding\n"); - fprintf(stderr, " -24516 = LOG & SQI & SIG thresholding\n"); - fprintf(stderr, "speckle remover: %d\n", speckle); - } + { + fprintf(stderr, "LOG = %5.2f\n", log); + fprintf(stderr, "SQI = %5.2f\n", sqi); + fprintf(stderr, "CSR = %5.2f\n", csr); + fprintf(stderr, "SIG = %5.2f\n", sig); + fprintf(stderr, "Calibration reflectivity: %5.2f dBZ\n", cal_dbz); + fprintf(stderr, "ZT flags: %d\n", z_flag_unc); /** can find these **/ + fprintf(stderr, "DZ flags: %d\n", z_flag_cor); /** defn in the **/ + fprintf(stderr, "VR flags: %d\n", v_flag); /** SIGMET Doc **/ + fprintf(stderr, "SW flags: %d\n", w_flag); + fprintf(stderr, "Flags: -3856 = SQI thresholding\n"); + fprintf(stderr, " -21846 = LOG thresholding\n"); + fprintf(stderr, " -24416 = LOG & SQI thresholding\n"); + fprintf(stderr, " -24516 = LOG & SQI & SIG thresholding\n"); + fprintf(stderr, "speckle remover: %d\n", speckle); + } if (radar_verbose_flag) - fprintf(stderr, "vel: %f prf: %f\n", max_vel, prf); + fprintf(stderr, "vel: %f prf: %f\n", max_vel, prf); /** Extracting Latitude and Longitude from nsig file **/ lat = nsig_from_fourb_ang(prod_file->rec2.ingest_head.lat_rad); @@ -411,7 +429,7 @@ RSL_nsig_to_radar if(lat > 180.0) lat -= 360.0; if(lon > 180.0) lon -= 360.0; if (radar_verbose_flag) - fprintf(stderr, "nsig_to_radar: lat %f, lon %f\n", lat, lon); + fprintf(stderr, "nsig_to_radar: lat %f, lon %f\n", lat, lon); /** Latitude deg, min, sec **/ latd = (int)lat; tmp = (lat - latd) * 60.0; @@ -426,10 +444,10 @@ RSL_nsig_to_radar /** Allocating memory for radar structure **/ radar = RSL_new_radar(MAX_RADAR_VOLUMES); if (radar == NULL) - { - fprintf(stderr, "nsig_to_radar: radar is NULL\n"); - return NULL; - } + { + fprintf(stderr, "nsig_to_radar: radar is NULL\n"); + return NULL; + } /** Filling Radar Header **/ radar->h.month = month; @@ -456,327 +474,436 @@ RSL_nsig_to_radar if (radar_verbose_flag) { #ifdef NSIG_VER2 - fprintf(stderr, "\nSIGMET version 2 raw product file.\n"); + fprintf(stderr, "\nSIGMET version 2 raw product file.\n"); #else - fprintf(stderr, "\nSIGMET version 1 raw product file.\n"); + fprintf(stderr, "\nSIGMET version 1 raw product file.\n"); #endif - fprintf(stderr, "Date: %2.2d/%2.2d/%4.4d %2.2d:%2.2d:%f\n", - radar->h.month, radar->h.day, radar->h.year, - radar->h.hour, radar->h.minute, radar->h.sec); - fprintf(stderr, "Name: "); - for (i=0; ih.name); i++) - fprintf(stderr, "%c", radar->h.name[i]); - fprintf(stderr, "\n"); - fprintf(stderr, "Lat/lon (%d %d' %d'', %d %d' %d'')\n", - radar->h.latd, radar->h.latm, radar->h.lats, - radar->h.lond, radar->h.lonm, radar->h.lons); + fprintf(stderr, "Date: %2.2d/%2.2d/%4.4d %2.2d:%2.2d:%f\n", + radar->h.month, radar->h.day, radar->h.year, + radar->h.hour, radar->h.minute, radar->h.sec); + fprintf(stderr, "Name: "); + for (i=0; ih.name); i++) + fprintf(stderr, "%c", radar->h.name[i]); + fprintf(stderr, "\n"); + fprintf(stderr, "Lat/lon (%d %d' %d'', %d %d' %d'')\n", + radar->h.latd, radar->h.latm, radar->h.lats, + radar->h.lond, radar->h.lonm, radar->h.lons); } /** Converting data **/ if (radar_verbose_flag) fprintf(stderr, "Expecting %d sweeps.\n", numsweep); for(i = 0; i < numsweep; i++) - { - nsig_sweep = nsig_read_sweep(fp, prod_file); - if (nsig_sweep == NULL) { /* EOF possibility */ - if (feof(fp)) break; - else continue; - } - if (rsl_qsweep != NULL) { - if (i > rsl_qsweep_max) break; - if (rsl_qsweep[i] == 0) continue; - } - if (radar_verbose_flag) - fprintf(stderr, "Read sweep # %d\n", i); - /* The whole sweep is 'nsig_sweep' ... pretty slick. + { + nsig_sweep = nsig_read_sweep(fp, prod_file); + if (nsig_sweep == NULL) { /* EOF possibility */ + if (feof(fp)) break; + else continue; + } + if (rsl_qsweep != NULL) { + if (i > rsl_qsweep_max) break; + if (rsl_qsweep[i] == 0) continue; + } + if (radar_verbose_flag) + fprintf(stderr, "Read sweep # %d\n", i); + /* The whole sweep is 'nsig_sweep' ... pretty slick. * * nsig_sweep[itype] -- [0..nparams], if non-null. */ - for (itype=0; itypeidh.time.month); - sweep_year = NSIG_I2(nsig_sweep[itype]->idh.time.year); - sweep_day = NSIG_I2(nsig_sweep[itype]->idh.time.day); - sweep_sec = NSIG_I4(nsig_sweep[itype]->idh.time.sec); + for (itype=0; itypeidh.time.month); + sweep_year = NSIG_I2(nsig_sweep[itype]->idh.time.year); + sweep_day = NSIG_I2(nsig_sweep[itype]->idh.time.day); + sweep_sec = NSIG_I4(nsig_sweep[itype]->idh.time.sec); #ifdef NSIG_VER2 - msec = NSIG_I2(nsig_sweep[itype]->idh.time.msec); - /* printf("....... msec == %d\n", msec); */ + msec = NSIG_I2(nsig_sweep[itype]->idh.time.msec); + /* printf("....... msec == %d\n", msec); */ #endif - /* converting seconds since mid to time of day */ - tmp = sweep_sec/3600.0; - sweep_hour = (int)tmp; - tmp = (tmp - sweep_hour) * 60.0; - sweep_minute = (int)tmp; - sweep_second = sweep_sec - (sweep_hour*3600 + sweep_minute*60); - - num_rays = NSIG_I2(nsig_sweep[itype]->idh.num_rays_exp); - - data_type = NSIG_I2(nsig_sweep[itype]->idh.data_type); - - ifield = 0; - switch (data_type) { - case NSIG_DTB_EXH: - ifield = -1; + /* converting seconds since mid to time of day */ + tmp = sweep_sec/3600.0; + sweep_hour = (int)tmp; + tmp = (tmp - sweep_hour) * 60.0; + sweep_minute = (int)tmp; + sweep_second = sweep_sec - (sweep_hour*3600 + sweep_minute*60); + + num_rays = NSIG_I2(nsig_sweep[itype]->idh.num_rays_exp); + + data_type = NSIG_I2(nsig_sweep[itype]->idh.data_type); + + ifield = 0; + switch (data_type) { + case NSIG_DTB_EXH: + ifield = -1; + break; + case NSIG_DTB_UCR: + case NSIG_DTB_UCR2: + ifield = ZT_INDEX; + f = ZT_F; + invf = ZT_INVF; + break; + case NSIG_DTB_CR: + case NSIG_DTB_CR2: + ifield = DZ_INDEX; + f = DZ_F; + invf = DZ_INVF; + break; + case NSIG_DTB_VEL: + case NSIG_DTB_VEL2: + ifield = VR_INDEX; + f = VR_F; + invf = VR_INVF; + break; + case NSIG_DTB_WID: + case NSIG_DTB_WID2: + ifield = SW_INDEX; + f = SW_F; + invf = SW_INVF; + break; + case NSIG_DTB_ZDR: + case NSIG_DTB_ZDR2: + ifield = DR_INDEX; + f = DR_F; + invf = DR_INVF; + break; + case NSIG_DTB_KDP: + ifield = KD_INDEX; + f = KD_F; + invf = KD_INVF; + break; + case NSIG_DTB_PHIDP: /* SRB 990127 */ + ifield = PH_INDEX; + f = PH_F; + invf = PH_INVF; + break; + case NSIG_DTB_RHOHV: /* SRB 000414 */ + ifield = RH_INDEX; + f = RH_F; + invf = RH_INVF; + break; + case NSIG_DTB_VELC: + case NSIG_DTB_VELC2: + ifield = VC_INDEX; + f = VC_F; + invf = VC_INVF; + break; + case NSIG_DTB_KDP2: + ifield = KD_INDEX; + f = KD_F; + invf = KD_INVF; + break; + case NSIG_DTB_PHIDP2: + ifield = PH_INDEX; + f = PH_F; + invf = PH_INVF; + break; + case NSIG_DTB_RHOHV2: + ifield = RH_INDEX; + f = RH_F; + invf = RH_INVF; + break; + case NSIG_DTB_SQI: + case NSIG_DTB_SQI2: + ifield = SQ_INDEX; + f = SQ_F; + invf = SQ_INVF; + break; + case NSIG_DTB_HCLASS: + case NSIG_DTB_HCLASS2: + ifield = HC_INDEX; + f = HC_F; + invf = HC_INVF; + break; + default: + fprintf(stderr,"Unknown field type: %d Skipping it.\n", data_type); + continue; + } + + if (radar_verbose_flag) + fprintf(stderr, " nsig_sweep[%d], data_type = %d, rays(expected) = %d, nrays(actual) = %d\n", itype, data_type, num_rays, NSIG_I2(nsig_sweep[itype]->idh.num_rays_act)); + + if (data_type != NSIG_DTB_EXH) { + if ((radar->v[ifield] == NULL)) { + if (rsl_qfield[ifield]) { + radar->v[ifield] = RSL_new_volume(numsweep); + radar->v[ifield]->h.f = f; + radar->v[ifield]->h.invf = invf; + } else { + /* Skip this field, because, the user does not want it. */ + continue; + } + } + if (radar->v[ifield]->sweep[i] == NULL) + radar->v[ifield]->sweep[i] = RSL_new_sweep(num_rays); + } + else + continue; /* Skip the actual extended header processing. + * This is different than getting it, so that + * the information is available for the other + * fields when filling the RSL ray headers. + */ + + /** DATA conversion time **/ + sweep = radar->v[ifield]->sweep[i]; + sweep->h.f = f; + sweep->h.invf = invf; + sweep->h.sweep_num = i; + sweep->h.beam_width = beam_width; + sweep->h.vert_half_bw = vert_half_bw; + sweep->h.horz_half_bw = horz_half_bw; + elev = nsig_from_bang(nsig_sweep[itype]->idh.fix_ang); + sweep->h.elev = elev; + + for(j = 0; j < num_rays; j++) + { + ray_p = nsig_sweep[itype]->ray[j]; + if (ray_p == NULL) continue; + bin_num = NSIG_I2(ray_p->h.num_bins); + + /* Load extended header information, if available. + * We need to pass the entire nsig_sweep and search for + * the extended header field (it may not be data_type==0). + */ + get_extended_header_info(nsig_sweep, xh_size, j, nparams, + &msec, &azm, &elev, + &pitch, &roll, &heading, + &azm_rate, &elev_rate, + &pitch_rate, &roll_rate, &heading_rate, + &lat, &lon, &alt, &rvc, + &vel_east, &vel_north, &vel_up); + + + if (radar->v[ifield]->sweep[i]->ray[j] == NULL) + radar->v[ifield]->sweep[i]->ray[j] = RSL_new_ray(bin_num); + ray = radar->v[ifield]->sweep[i]->ray[j]; + ray->h.f = f; + ray->h.invf = invf; + /** Ray is at nsig_sweep[itype].ray->... **/ + /** Loading nsig data into data structure **/ + + ray->h.month = sweep_month; + ray->h.day = sweep_day; + ray->h.year = sweep_year; /* Year 2000 compliant. */ + ray->h.hour = sweep_hour; + ray->h.minute = sweep_minute; + if (msec == 0) { /* No extended header */ + ray->h.sec = NSIG_I2(ray_p->h.sec) + sweep_second; + elev = sweep->h.elev; + } else + ray->h.sec = sweep_second + msec/1000.0; + + /* add time ... handles end of min,hour,month,year and century. */ + if (ray->h.sec >= 60) /* Should I fix the time no matter what? */ + RSL_fix_time(ray); /* Repair second overflow. */ + + ray->h.ray_num = j; + ray->h.elev_num = i; + ray->h.range_bin1 = (int)rng_first_bin; + ray->h.gate_size = (int)(bin_space+.5); /* Nearest int */ + ray->h.vel_res = bin_space; + ray->h.sweep_rate = sweep_rate; + ray->h.prf = (int)prf; + if (prf != 0) + ray->h.unam_rng = 299793000.0 / (2.0 * prf * 1000.0); /* km */ + else + ray->h.unam_rng = 0.0; + ray->h.fix_angle = (float)sweep->h.elev; + ray->h.azim_rate = azim_rate; + ray->h.pulse_count = (float)num_samples; + ray->h.pulse_width = pw; + ray->h.beam_width = beam_width; + ray->h.frequency = freq / 1000.0; /* GHz */ + ray->h.wavelength = wave/100.0; /* meters */ + ray->h.nyq_vel = max_vel; /* m/s */ + if (elev == 0.) elev = sweep->h.elev; + ray->h.elev = elev; + /* Compute mean azimuth angle for ray. */ + az1 = nsig_from_bang(ray_p->h.beg_azm); + az2 = nsig_from_bang(ray_p->h.end_azm); + /* printf("az1, %f, az2 %f\n", az1, az2); */ + if(az1 > az2) + if((az1 - az2) > 180.0) az2 += 360.0; + else + ; + else + if((az2 - az1) > 180.0) az1 += 360.0; + + az1 = (az1 + az2) / 2.0; + if (az1 > 360) az1 -= 360; + ray->h.azimuth = az1; + + /* From the extended header information, we learn the following. */ + ray->h.pitch = pitch; + ray->h.roll = roll; + ray->h.heading = heading; + ray->h.pitch_rate = pitch_rate; + ray->h.roll_rate = roll_rate; + ray->h.heading_rate = heading_rate; + ray->h.lat = lat; + ray->h.lon = lon; + ray->h.alt = alt; + ray->h.rvc = rvc; + ray->h.vel_east = vel_east; + ray->h.vel_north = vel_north; + ray->h.vel_up = vel_up; + + /* printf("Processing sweep[%d]->ray[%d]: %d %f %f %f %f %f %f %f %f %d nbins=%d, bin1=%d gate=%d\n", + i, j, msec, ray->h.sec, ray->h.azimuth, ray->h.elev, ray->h.pitch, ray->h.roll, ray->h.heading, ray->h.lat, ray->h.lon, ray->h.alt, ray->h.nbins, ray->h.range_bin1, ray->h.gate_size); + */ + /* TODO: ingest data header contains a value for bits-per-bin. + * This might be of use to allocate an array for ray->range with + * either 1-byte or 2-byte elements. Then there's no need for + * memmove() whenever we need 2 bytes. + */ + + if (data_type == NSIG_DTB_EXH) continue; + ray_data = 0; + for(k = 0; k < bin_num; k++) { + switch(data_type) { + case NSIG_DTB_UCR: + case NSIG_DTB_CR: + if (ray_p->range[k] == 0) ray_data = NSIG_NO_ECHO; + else ray_data = (float)((ray_p->range[k]-64.0)/2.0); + break; + /* Simplified the velocity conversion for NSIG_DTB_VEL, using + * formula from IRIS Programmer's Manual. BLK, Oct 9 2009. + */ + case NSIG_DTB_VEL: + if (ray_p->range[k] == 0) ray_data = NSIG_NO_ECHO; + else ray_data = (float)((ray_p->range[k]-128.0)/127.0)*max_vel; + break; + + case NSIG_DTB_WID: + if (ray_p->range[k] == 0) ray_data = NSIG_NO_ECHO; + else ray_data =(float)((ray_p->range[k])/256.0)*max_vel; + break; + + case NSIG_DTB_ZDR: + if (ray_p->range[k] == 0) ray_data = NSIG_NO_ECHO; + else ray_data = (float)((ray_p->range[k]-128.0)/16.0); + break; + + case NSIG_DTB_KDP: + if (ray_p->range[k] == 0 || ray_p->range[k] == 255 || + rsl_kdp_wavelen == 0.0) { + ray_data = NSIG_NO_ECHO; break; - case NSIG_DTB_UCR: - ifield = ZT_INDEX; - f = ZT_F; - invf = ZT_INVF; - break; - case NSIG_DTB_CR: - ifield = DZ_INDEX; - f = DZ_F; - invf = DZ_INVF; - break; - case NSIG_DTB_VEL: - ifield = VR_INDEX; - f = VR_F; - invf = VR_INVF; - break; - case NSIG_DTB_WID: - ifield = SW_INDEX; - f = SW_F; - invf = SW_INVF; - break; - case NSIG_DTB_ZDR: - ifield = DR_INDEX; - f = DR_F; - invf = DR_INVF; - break; - case NSIG_DTB_KDP: - ifield = KD_INDEX; - f = KD_F; - invf = KD_INVF; - break; - case NSIG_DTB_PHIDP: /* SRB 990127 */ - ifield = PH_INDEX; - f = PH_F; - invf = PH_INVF; - break; - case NSIG_DTB_RHOHV: /* SRB 000414 */ - ifield = RH_INDEX; - f = RH_F; - invf = RH_INVF; - break; - case NSIG_DTB_SQI: - ifield = SQ_INDEX; - f = SQ_F; - invf = SQ_INVF; - break; - default: - fprintf(stderr,"Unknown field type: %d Skipping it.\n", data_type); - continue; - } - - if (radar_verbose_flag) - fprintf(stderr, " nsig_sweep[%d], data_type = %d, rays(expected) = %d, nrays(actual) = %d\n", itype, data_type, num_rays, NSIG_I2(nsig_sweep[itype]->idh.num_rays_act)); - - if (data_type != NSIG_DTB_EXH) { - if ((radar->v[ifield] == NULL)) { - if (rsl_qfield[ifield]) { - radar->v[ifield] = RSL_new_volume(numsweep); - radar->v[ifield]->h.f = f; - radar->v[ifield]->h.invf = invf; - } else { - /* Skip this field, because, the user does not want it. */ - continue; - } - } - if (radar->v[ifield]->sweep[i] == NULL) - radar->v[ifield]->sweep[i] = RSL_new_sweep(num_rays); - } - else - continue; /* Skip the actual extended header processing. - * This is different than getting it, so that - * the information is available for the other - * fields when filling the RSL ray headers. - */ - - /** DATA conversion time **/ - sweep = radar->v[ifield]->sweep[i]; - sweep->h.f = f; - sweep->h.invf = invf; - sweep->h.sweep_num = i; - sweep->h.beam_width = beam_width; - sweep->h.vert_half_bw = vert_half_bw; - sweep->h.horz_half_bw = horz_half_bw; - elev = nsig_from_bang(nsig_sweep[itype]->idh.fix_ang); - sweep->h.elev = elev; - - for(j = 0; j < num_rays; j++) - { - ray_p = nsig_sweep[itype]->ray[j]; - if (ray_p == NULL) continue; - bin_num = NSIG_I2(ray_p->h.num_bins); - - /* Load extended header information, if available. - * We need to pass the entire nsig_sweep and search for - * the extended header field (it may not be data_type==0). - */ - get_extended_header_info(nsig_sweep, xh_size, j, nparams, - &msec, &azm, &elev, - &pitch, &roll, &heading, - &azm_rate, &elev_rate, - &pitch_rate, &roll_rate, &heading_rate, - &lat, &lon, &alt, &rvc, - &vel_east, &vel_north, &vel_up); - - - if (radar->v[ifield]->sweep[i]->ray[j] == NULL) - radar->v[ifield]->sweep[i]->ray[j] = RSL_new_ray(bin_num); - ray = radar->v[ifield]->sweep[i]->ray[j]; - ray->h.f = f; - ray->h.invf = invf; - /** Ray is at nsig_sweep[itype].ray->... **/ - /** Loading nsig data into data structure **/ - - ray->h.month = sweep_month; - ray->h.day = sweep_day; - ray->h.year = sweep_year; /* Year 2000 compliant. */ - ray->h.hour = sweep_hour; - ray->h.minute = sweep_minute; - if (msec == 0) { /* No extended header */ - ray->h.sec = NSIG_I2(ray_p->h.sec) + sweep_second; - elev = sweep->h.elev; - } else - ray->h.sec = sweep_second + msec/1000.0; - - /* add time ... handles end of min,hour,month,year and century. */ - if (ray->h.sec >= 60) /* Should I fix the time no matter what? */ - RSL_fix_time(ray); /* Repair second overflow. */ - - ray->h.ray_num = j; - ray->h.elev_num = i; - ray->h.range_bin1 = (int)rng_first_bin; - ray->h.gate_size = (int)(bin_space+.5); /* Nearest int */ - ray->h.vel_res = bin_space; - ray->h.sweep_rate = sweep_rate; - ray->h.prf = (int)prf; - if (prf != 0) - ray->h.unam_rng = 299793000.0 / (2.0 * prf * 1000.0); /* km */ - else - ray->h.unam_rng = 0.0; - ray->h.fix_angle = (float)sweep->h.elev; - ray->h.azim_rate = azim_rate; - ray->h.pulse_count = (float)num_samples; - ray->h.pulse_width = pw; - ray->h.beam_width = beam_width; - ray->h.frequency = freq / 1000.0; /* GHz */ - ray->h.wavelength = wave/100.0; /* meters */ - ray->h.nyq_vel = max_vel; /* m/s */ - if (elev == 0.) elev = sweep->h.elev; - ray->h.elev = elev; - /* Compute mean azimuth angle for ray. */ - az1 = nsig_from_bang(ray_p->h.beg_azm); - az2 = nsig_from_bang(ray_p->h.end_azm); - /* printf("az1, %f, az2 %f\n", az1, az2); */ - if(az1 > az2) - if((az1 - az2) > 180.0) az2 += 360.0; - else - ; - else - if((az2 - az1) > 180.0) az1 += 360.0; - - az1 = (az1 + az2) / 2.0; - if (az1 > 360) az1 -= 360; - ray->h.azimuth = az1; - - /* From the extended header information, we learn the following. */ - ray->h.pitch = pitch; - ray->h.roll = roll; - ray->h.heading = heading; - ray->h.pitch_rate = pitch_rate; - ray->h.roll_rate = roll_rate; - ray->h.heading_rate = heading_rate; - ray->h.lat = lat; - ray->h.lon = lon; - ray->h.alt = alt; - ray->h.rvc = rvc; - ray->h.vel_east = vel_east; - ray->h.vel_north = vel_north; - ray->h.vel_up = vel_up; - - /* printf("Processing sweep[%d]->ray[%d]: %d %f %f %f %f %f %f %f %f %d nbins=%d, bin1=%d gate=%d\n", - i, j, msec, ray->h.sec, ray->h.azimuth, ray->h.elev, ray->h.pitch, ray->h.roll, ray->h.heading, ray->h.lat, ray->h.lon, ray->h.alt, ray->h.nbins, ray->h.range_bin1, ray->h.gate_size); - */ - if (data_type == NSIG_DTB_EXH) continue; - ray_data = 0; - for(k = 0; k < bin_num; k++) { - switch(data_type) { - case NSIG_DTB_UCR: - case NSIG_DTB_CR: - if (ray_p->range[k] == 0) ray_data = NSIG_NO_ECHO; - else ray_data = (float)((ray_p->range[k]-64.0)/2.0); - break; - - case NSIG_DTB_VEL: - if (ray_p->range[k] == 0) ray_data = NSIG_NO_ECHO; - else ray_data = (float)((ray_p->range[k]*max_vel/127.0)+ - max_vel*(1.0-255.0/127.0)); - break; - - case NSIG_DTB_WID: - if (ray_p->range[k] == 0) ray_data = NSIG_NO_ECHO; - else ray_data =(float)((ray_p->range[k])/256.0)*max_vel; - break; - - case NSIG_DTB_ZDR: - if (ray_p->range[k] == 0) ray_data = NSIG_NO_ECHO; - else ray_data = (float)((ray_p->range[k]-128.0)/16.0); - break; - - /* - * Special optimization note: - * For KDP, PHIDP, RHOHV we skip the float conversion, - * and carry the native sigmet data values into RSL storage. - */ - case NSIG_DTB_KDP: - ray_data = ray_p->range[k]; - /* F_OFFSET *must* match the definition in volume.c */ -#define F_OFFSET 4 - if (ray_data == 0 || ray_data == 255) ray_data = NSIG_NO_ECHO; - else ray_data += F_OFFSET; - break; - - case NSIG_DTB_RHOHV: - case NSIG_DTB_PHIDP: - if (ray_p->range[k] <= 0 || ray_p->range[k] >= 255) - ray_data = NSIG_NO_ECHO; - else - ray_data = ray_p->range[k]; - break; - case NSIG_DTB_SQI: - if (ray_p->range[k] == 0) ray_data = NSIG_NO_ECHO; - else ray_data = - (float)sqrt((ray_p->range[k]-1.0)/253.0); - } - - if (ray_data == NSIG_NO_ECHO) - ray->range[k] = ray->h.invf(BADVAL); - else - if ( (data_type == NSIG_DTB_KDP) || - (data_type == NSIG_DTB_RHOHV) || - (data_type == NSIG_DTB_PHIDP) ) - ray->range[k] = (Range)ray_data; - else - ray->range[k] = ray->h.invf(ray_data); - - /* - if (data_type == NSIG_DTB_KDP) - printf("v[%d]->sweep[%d]->ray[%d]->range[%d] = %f, %d, %f\n", - ifield, i, j, k, ray->h.f(ray->range[k]), - (int)ray_p->range[k], ray_data); - */ - } } - } + if (ray_p->range[k] < 128) + ray_data = (-0.25 * + pow((double)600.0,(double)((127-ray_p->range[k])/126.0))) / + rsl_kdp_wavelen; + else if (ray_p->range[k] > 128) + ray_data = (0.25 * + pow((double)600.0,(double)((ray_p->range[k]-129)/126.0))) / + rsl_kdp_wavelen; + else + ray_data = 0.0; + break; + + case NSIG_DTB_PHIDP: + if (ray_p->range[k] == 0 || ray_p->range[k] == 255) + ray_data = NSIG_NO_ECHO; + else + ray_data = 180.0*((ray_p->range[k]-1.0)/254.0); + break; + + case NSIG_DTB_RHOHV: + if (ray_p->range[k] == 0 || ray_p->range[k] == 255) + ray_data = NSIG_NO_ECHO; + else + ray_data = sqrt((double)((ray_p->range[k]-1.0)/253.0)); + break; + + case NSIG_DTB_HCLASS: + if (ray_p->range[k] == 0 || ray_p->range[k] == 255) + ray_data = NSIG_NO_ECHO; + else + ray_data = ray_p->range[k]; + break; + + case NSIG_DTB_SQI: + if (ray_p->range[k] == 0) ray_data = NSIG_NO_ECHO; + else ray_data = (float)sqrt((ray_p->range[k]-1.0)/253.0); + break; + + case NSIG_DTB_VELC: + if (ray_p->range[k] == 0) ray_data = NSIG_NO_ECHO; + else { + incr=75./127.; /* (+|- 75m/s) / 254 values */ + ray_data = (float)(ray_p->range[k]-128)*incr; + } + break; + + case NSIG_DTB_UCR2: + case NSIG_DTB_CR2: + case NSIG_DTB_VEL2: + case NSIG_DTB_VELC2: + case NSIG_DTB_ZDR2: + case NSIG_DTB_KDP2: + memmove(&nsig_u2byte, &ray_p->range[2*k], 2); + nsig_u2byte = NSIG_I2(&nsig_u2byte); + if (nsig_u2byte == 0 || nsig_u2byte == 65535) + ray_data = NSIG_NO_ECHO2; + else ray_data = (float)(nsig_u2byte-32768)/100.; + break; + + case NSIG_DTB_WID2: + memmove(&nsig_u2byte, &ray_p->range[2*k], 2); + nsig_u2byte = NSIG_I2(&nsig_u2byte); + if (nsig_u2byte == 0 || nsig_u2byte == 65535) + ray_data = NSIG_NO_ECHO2; + else ray_data = (float)nsig_u2byte/100.; + break; + + case NSIG_DTB_PHIDP2: + memmove(&nsig_u2byte, &ray_p->range[2*k], 2); + nsig_u2byte = NSIG_I2(&nsig_u2byte); + if (nsig_u2byte == 0 || nsig_u2byte == 65535) + ray_data = NSIG_NO_ECHO; + else + ray_data = 360.*(nsig_u2byte-1)/65534.; + break; + + case NSIG_DTB_SQI2: + case NSIG_DTB_RHOHV2: + memmove(&nsig_u2byte, &ray_p->range[2*k], 2); + nsig_u2byte = NSIG_I2(&nsig_u2byte); + if (nsig_u2byte == 0 || nsig_u2byte == 65535) + ray_data = NSIG_NO_ECHO2; + else ray_data = (float)(nsig_u2byte-1)/65533.; + break; + + case NSIG_DTB_HCLASS2: + memmove(&nsig_u2byte, &ray_p->range[2*k], 2); + nsig_u2byte = NSIG_I2(&nsig_u2byte); + if (nsig_u2byte == 0 || nsig_u2byte == 65535) + ray_data = NSIG_NO_ECHO2; + else + ray_data = nsig_u2byte; + } + + if (ray_data == NSIG_NO_ECHO || ray_data == NSIG_NO_ECHO2) + ray->range[k] = ray->h.invf(BADVAL); + else + ray->range[k] = ray->h.invf(ray_data); + + /* + if (data_type == NSIG_DTB_KDP) + printf("v[%d]->sweep[%d]->ray[%d]->range[%d] = %f, %d, %f\n", + ifield, i, j, k, ray->h.f(ray->range[k]), + (int)ray_p->range[k], ray_data); + */ + } + } + } nsig_free_sweep(nsig_sweep); - } + } /* Do not reset radar->h.nvolumes. It is already set properly. */ if (radar_verbose_flag) - fprintf(stderr, "Max index of radar->v[0..%d]\n", radar->h.nvolumes); + fprintf(stderr, "Max index of radar->v[0..%d]\n", radar->h.nvolumes); /** close nsig file **/ diff --git a/prune.c b/prune.c index 6ee0459..1a23913 100644 --- a/prune.c +++ b/prune.c @@ -36,6 +36,21 @@ #include "rsl.h" extern int radar_verbose_flag; +/* Define global variable for pruning and the functions to set or unset it. + * Added by Bart Kelley, SSAI, August 26, 2009 + */ +int prune_radar = 1; + +void RSL_prune_radar_on() +{ + prune_radar = 1; +} + +void RSL_prune_radar_off() +{ + prune_radar = 0; +} + Ray *RSL_prune_ray(Ray *ray) { if (ray == NULL) return NULL; @@ -50,20 +65,20 @@ Sweep *RSL_prune_sweep(Sweep *s) if (s == NULL) return NULL; if (s->h.nrays == 0) { - RSL_free_sweep(s); - return NULL; + RSL_free_sweep(s); + return NULL; } /* * Squash out all dataless rays. 'j' is the index for the squashed (pruned) * rays. */ for (i=0,j=0; ih.nrays; i++) - if ((s->ray[i] = RSL_prune_ray(s->ray[i]))) - s->ray[j++] = s->ray[i]; /* Keep this ray. */ + if ((s->ray[i] = RSL_prune_ray(s->ray[i]))) + s->ray[j++] = s->ray[i]; /* Keep this ray. */ if (j==0) { - RSL_free_sweep(s); - return NULL; /* All rays were pruned. */ + RSL_free_sweep(s); + return NULL; /* All rays were pruned. */ } for (i=j; ih.nrays; i++) s->ray[i] = NULL; s->h.nrays = j; @@ -76,19 +91,19 @@ Volume *RSL_prune_volume(Volume *v) if (v == NULL) return NULL; if (v->h.nsweeps == 0) { - RSL_free_volume(v); - return NULL; + RSL_free_volume(v); + return NULL; } /* * Squash out all dataless sweeps. 'j' is the index for sweep containing data. */ for (i=0,j=0; ih.nsweeps; i++) - if ((v->sweep[i] = RSL_prune_sweep(v->sweep[i]))) - v->sweep[j++] = v->sweep[i]; /* Keep this sweep. */ + if ((v->sweep[i] = RSL_prune_sweep(v->sweep[i]))) + v->sweep[j++] = v->sweep[i]; /* Keep this sweep. */ if (j==0) { - RSL_free_volume(v); - return NULL; /* All sweeps were pruned. */ + RSL_free_volume(v); + return NULL; /* All sweeps were pruned. */ } for (i=j; ih.nsweeps; i++) v->sweep[i] = NULL; v->h.nsweeps = j; @@ -100,8 +115,11 @@ Radar *RSL_prune_radar(Radar *radar) int i; /* Volume indexes are fixed so we just prune the substructures. */ if (radar == NULL) return NULL; - for (i=0; ih.nvolumes; i++) - radar->v[i] = RSL_prune_volume(radar->v[i]); + if (prune_radar) + for (i=0; ih.nvolumes; i++) + radar->v[i] = RSL_prune_volume(radar->v[i]); + else if (radar_verbose_flag) fprintf(stderr, + "RSL_prune_radar: No pruning done. prune_radar = %d\n", prune_radar); return radar; } diff --git a/volume.c b/volume.c index 4c17f95..2ff6aa9 100644 --- a/volume.c +++ b/volume.c @@ -99,6 +99,10 @@ extern int radar_verbose_flag; #define F_DZ_RANGE_OFFSET 32 #endif +/* #define F_VR_OFFSET 63.5 */ +#define F_VR_OFFSET 127.0 +#define F_DR_OFFSET 12.0 + /* IMPORTANT: This is the offset from reserved values. This * number must be exactly (or >=) the number of * reserved values in XX_F and XX_INVF. @@ -121,7 +125,7 @@ float DZ_F(Range x) { float VR_F(Range x) { float val; if (x >= F_OFFSET) { /* This test works when Range is unsigned. */ - val = (((float)x-F_OFFSET)/F_FACTOR - 63.5); /* Default wsr88d coding. */ + val = (((float)x-F_OFFSET)/F_FACTOR - F_VR_OFFSET); /* Default wsr88d coding. */ /* fprintf(stderr, "x=%d, val=%f\n", x, val); */ return val; } @@ -135,7 +139,7 @@ float VR_F(Range x) { float DR_F(Range x) { /* Differential reflectivity */ float val; if (x >= F_OFFSET) { /* This test works when Range is unsigned. */ - val = (((float)x-F_OFFSET)/F_DR_FACTOR - 12.0); + val = (((float)x-F_OFFSET)/F_DR_FACTOR - F_DR_OFFSET); return val; } if (x == 0) return BADVAL; @@ -155,6 +159,11 @@ float LR_F(Range x) {/* From MCTEX */ return BADVAL; } +float HC_F(Range x) { /* HydroClass (Sigmet) */ + if (x == 0) return BADVAL; + return (float)x; +} + /**************************** Sigmet RhoHV : one_byte values > RohHV = sqrt((N-1)/253) @@ -167,8 +176,9 @@ float LR_F(Range x) {/* From MCTEX */ *******************************/ float RH_F(Range x) { if (x == 0) return BADVAL; - return (float)(sqrt((double)((x-1.0)/253.0))); - } + /* return (float)(sqrt((double)((x-1.0)/253.0))); */ + return (float)(x-1) / 65533.; +} /***************************** Sigmet PhiDP : one_byte values @@ -182,42 +192,51 @@ float RH_F(Range x) { ******************************/ float PH_F(Range x) { if (x == 0) return BADVAL; - return (float)(180.0*((x-1.0)/254.0)); + /*return (float)(180.0*((x-1.0)/254.0));*/ + return (360.*(x-1.))/65534.; } +/* TODO: Should this be 5. cm. instead of 0.5? Or maybe 10. cm.? */ float rsl_kdp_wavelen = 0.5; /* Default radar wavelen = .5 cm. See * nsig_to_radar.c for the code that sets this. - */ + */ +/* KD_F for 1 or 2 byte. */ float KD_F(Range x) { +/****** Commented-out code for 1-byte Sigmet native data format. if (x >= F_OFFSET) { - x -= F_OFFSET; - if (rsl_kdp_wavelen == 0.0) return BADVAL; - if (x < 128) - return (float)( - -0.25 * pow((double)600.0,(double)((127-x)/126.0)) - )/rsl_kdp_wavelen; - else if (x > 128) - return (float)( - 0.25 * pow((double)600.0,(double)((x-129)/126.0)) - )/rsl_kdp_wavelen; - else - return 0.0; + x -= F_OFFSET; + if (rsl_kdp_wavelen == 0.0) return BADVAL; + if (x < 128) + return (float)( + -0.25 * pow((double)600.0,(double)((127-x)/126.0)) + )/rsl_kdp_wavelen; + else if (x > 128) + return (float)( + 0.25 * pow((double)600.0,(double)((x-129)/126.0)) + )/rsl_kdp_wavelen; + else + return 0.0; } - if (x == 0) return BADVAL; if (x == 1) return RFVAL; if (x == 2) return APFLAG; if (x == 3) return NOECHO; return BADVAL; +******/ + + if (x == 0) return BADVAL; + return (x-32768.)/100.; +} + +float NP_F(Range x) { /* Normalized Coherent Power (DORADE) */ + if (x == 0) return BADVAL; + return (float)x / 100.; } +/* Signal Quality Index */ float SQ_F(Range x) { - if (x >= F_OFFSET) return (float)((x - F_OFFSET)/10000.); if (x == 0) return BADVAL; - if (x == 1) return RFVAL; - if (x == 2) return APFLAG; - if (x == 3) return NOECHO; - return BADVAL; /* Can't get here, but quiets the compiler. */ + return (float)(x-1) / 65533.; } float TI_F(Range x) @@ -246,6 +265,11 @@ float CH_F(Range x) { return DZ_F(x); } float AH_F(Range x) { return DZ_F(x); } float CV_F(Range x) { return DZ_F(x); } float AV_F(Range x) { return DZ_F(x); } +float VS_F(Range x) { return VR_F(x); } +float VL_F(Range x) { return VR_F(x); } +float VG_F(Range x) { return VR_F(x); } +float VT_F(Range x) { return VR_F(x); } +float VC_F(Range x) { return VR_F(x); } @@ -260,6 +284,7 @@ Range DZ_INVF(float x) if (x == RFVAL) return (Range)1; if (x == APFLAG) return (Range)2; if (x == NOECHO) return (Range)3; + if (x < -F_DZ_RANGE_OFFSET) return (Range)0; return (Range)(F_FACTOR*(x+F_DZ_RANGE_OFFSET)+.5 + F_OFFSET); /* Default wsr88d. */ } @@ -269,7 +294,8 @@ Range VR_INVF(float x) if (x == RFVAL) return (Range)1; if (x == APFLAG) return (Range)2; if (x == NOECHO) return (Range)3; - return (Range)(F_FACTOR*(x+63.5)+.5 + F_OFFSET); /* Default wsr88d coding. */ + if (x < -F_VR_OFFSET) return (Range)0; + return (Range)(F_FACTOR*(x+F_VR_OFFSET)+.5 + F_OFFSET); /* Default wsr88d coding. */ } Range DR_INVF(float x) /* Differential reflectivity */ @@ -278,7 +304,14 @@ Range DR_INVF(float x) /* Differential reflectivity */ if (x == RFVAL) return (Range)1; if (x == APFLAG) return (Range)2; if (x == NOECHO) return (Range)3; - return (Range)(F_DR_FACTOR*(x + 12.0) + F_OFFSET + 0.5); + if (x < -F_DR_OFFSET) return (Range)0; + return (Range)(F_DR_FACTOR*(x + F_DR_OFFSET) + F_OFFSET + 0.5); +} + +Range HC_INVF(float x) /* HydroClass (Sigmet) */ +{ + if (x == BADVAL) return (Range)0; + return (Range)(x + 0.5); /* Round */ } Range LR_INVF(float x) /* MCTEX */ @@ -300,9 +333,11 @@ Range LR_INVF(float x) /* MCTEX */ > 253 0.9980 > 254 1.0000 ****************************/ +/* RH_INVF for 1 or 2 byte data. */ Range RH_INVF(float x) { if (x == BADVAL) return (Range)0; - return (Range)(x * x * 253.0 + 1.0 + 0.5); + /* return (Range)(x * x * 253.0 + 1.0 + 0.5); */ + return (Range)(x * 65533. + 1. +.5); } /****************************** @@ -317,11 +352,15 @@ Range RH_INVF(float x) { *******************************/ Range PH_INVF(float x) { if (x == BADVAL) return (Range)0; - return (Range)((x / 180.0) * 254.0 + 1.0 + 0.5); + /* return (Range)((x / 180.0) * 254.0 + 1.0 + 0.5); */ + return (Range)(x*65534./360. + 1.0 + 0.5); } + +/* KD_INVF for 1 or 2 byte data. */ Range KD_INVF(float x) { if (x == BADVAL) return (Range)0; +/****** Commented-out code for 1-byte Sigmet native data format. if (x == RFVAL) return (Range)1; if (x == APFLAG) return (Range)2; if (x == NOECHO) return (Range)3; @@ -340,17 +379,23 @@ Range KD_INVF(float x) { x = 128; } x += F_OFFSET; - return (Range) x; +******/ + return (Range)(x * 100. + 32768. + 0.5); } -Range SQ_INVF(float x) /* Signal Quality Index */ +/* Signal Quality Index */ +Range SQ_INVF(float x) { if (x == BADVAL) return (Range)0; - if (x == RFVAL) return (Range)1; - if (x == APFLAG) return (Range)2; - if (x == NOECHO) return (Range)3; - return (Range)(x * 10000. + F_OFFSET); + return (Range)(x * 65533. + 1. +.5); +} + + +Range NP_INVF(float x) /* Normalized Coherent Power (DORADE) */ +{ + if (x == BADVAL) return (0); + return (Range)(x * 100.); } Range TI_INVF(float x) /* MCTEX */ @@ -379,6 +424,11 @@ Range CH_INVF(float x) { return DZ_INVF(x); } Range AH_INVF(float x) { return DZ_INVF(x); } Range CV_INVF(float x) { return DZ_INVF(x); } Range AV_INVF(float x) { return DZ_INVF(x); } +Range VS_INVF(float x) { return VR_INVF(x); } +Range VL_INVF(float x) { return VR_INVF(x); } +Range VG_INVF(float x) { return VR_INVF(x); } +Range VT_INVF(float x) { return VR_INVF(x); } +Range VC_INVF(float x) { return VR_INVF(x); } @@ -1332,7 +1382,9 @@ int rsl_qfield[MAX_RADAR_VOLUMES] = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, - 1, 1, 1, 1, 1 + 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, + 1, 1 };