#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);
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;
/* Determine where 'itype' for extended header is. */
for (itype = 0; itype<nparams; itype++) {
data_type = NSIG_I2(nsig_sweep[itype]->idh.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. */
*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);
#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;
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.
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);
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)
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;
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;
/** 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);
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;
/** 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;
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; i<sizeof(radar->h.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; i<sizeof(radar->h.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; itype<nparams; itype++) {
- if (nsig_sweep[itype] == NULL) continue;
-
- /** Reading date and time **/
- sweep_month = NSIG_I2(nsig_sweep[itype]->idh.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; itype<nparams; itype++) {
+ if (nsig_sweep[itype] == NULL) continue;
+
+ /** Reading date and time **/
+ sweep_month = NSIG_I2(nsig_sweep[itype]->idh.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 **/
projects / ~andy / rsl / blobdiff