X-Git-Url: http://pileus.org/git/?p=~andy%2Frsl;a=blobdiff_plain;f=nsig_to_radar.c;h=8454ab1771c5f91734ca10354934484c5b978a2d;hp=1e8b39095afa2adadb9633d1829151d6338d079c;hb=e953e14aeb8c5b83435ee0103f3b069cc1875c96;hpb=012916676d26251849e408aaf574458e196df2c4 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 **/