From d08a7f8a699a044bc4ac5f93917aa7f6c463923b Mon Sep 17 00:00:00 2001 From: Andy Spencer Date: Wed, 2 Feb 2011 03:57:21 +0000 Subject: [PATCH] Changes from Bart (2011-02-01) --- CHANGES | 17 +- Makefile.am | 2 +- configure.in | 2 +- doc/RSL_clear.html | 3 +- doc/RSL_copy.html | 3 +- doc/RSL_get_win.html | 9 +- doc/RSL_new.html | 3 +- doc/RSL_read.html | 3 +- dorade_to_radar.c | 414 +++++++++++--- examples/Makefile.am | 4 +- examples/qlook.c | 1079 ++++++++++++++++++----------------- examples/qlook_usage.c | 139 +++-- examples/wsr_hist_uf_test.c | 20 +- mkinstalldirs | 40 ++ nsig_to_radar.c | 71 ++- prune.c | 46 +- radar.c | 1 + radar_to_uf.c | 647 +++++++++++---------- rsl.h | 137 ++++- uf_to_radar.c | 670 +++++++++++----------- volume.c | 380 ++++++------ wsr88d.c | 3 + wsr88d_locations.dat | 1 + wsr88d_m31.c | 592 +++++++------------ wsr88d_to_radar.c | 22 +- 25 files changed, 2319 insertions(+), 1989 deletions(-) create mode 100644 mkinstalldirs diff --git a/CHANGES b/CHANGES index 9266df3..fb2b03f 100644 --- a/CHANGES +++ b/CHANGES @@ -4,16 +4,15 @@ * v1.41 In progress * * 1. wsr88d_m31.c: Simplified the WSR-88D ray structure and supporting code. - * Removed function read_data_moment. + * Removed function read_data_moment. Added support for dual-polarization + * data fields. * 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() + * 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. + * Modified nsig_to_radar.c and rsl.h to handle Sigmet dual PRF. Much + * thanks to Fabio Sato and Cesar Beneti for the dual PRF code. + * 3. Completed DORADE ingest. Files involved: dorade_to_radar.c, dorade.c, + * volume.c. *--------------------------------------------------------------------- * v1.40 Released 10/10/2008 * diff --git a/Makefile.am b/Makefile.am index 361f169..d859d71 100644 --- a/Makefile.am +++ b/Makefile.am @@ -8,7 +8,7 @@ colordir = $(libdir)/colors lib_LTLIBRARIES = librsl.la -librsl_la_LDFLAGS = -version-info 1:40 +librsl_la_LDFLAGS = -version-info 1:41 librsl_la_SOURCES = \ $(rapic_c) $(radtec_c)\ dorade.c dorade_print.c dorade_to_radar.c\ diff --git a/configure.in b/configure.in index e0d59af..62a6ae5 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.1) +AC_INIT(rsl, v1.40.3) AC_CONFIG_SRCDIR(volume.c) AM_INIT_AUTOMAKE diff --git a/doc/RSL_clear.html b/doc/RSL_clear.html index ab8e55d..4503182 100644 --- a/doc/RSL_clear.html +++ b/doc/RSL_clear.html @@ -26,9 +26,8 @@ Upon successful completion, a pointer to the appropriate structure is returned.

See also

-Example
RSL_free_volume, RSL_free_sweep, RSL_free_ray
-

Author: John H. Merritt +

Author: John H. Merritt diff --git a/doc/RSL_copy.html b/doc/RSL_copy.html index 7a182a8..6f1ef1c 100644 --- a/doc/RSL_copy.html +++ b/doc/RSL_copy.html @@ -26,9 +26,8 @@ Upon successful completion, a pointer to the appropriate structure is returned.

See also

-Example
RSL_free_volume, RSL_free_sweep, RSL_free_ray
-

Author: John H. Merritt +

Author: John H. Merritt diff --git a/doc/RSL_get_win.html b/doc/RSL_get_win.html index e5019df..a43b5fc 100644 --- a/doc/RSL_get_win.html +++ b/doc/RSL_get_win.html @@ -19,7 +19,10 @@

Description

-RSL_get_window_from_radar calls RSL_get_window_from_volume for the number of volumes. RSL_get_window_from_volume calls RSL_get_window_from_sweep for the number of sweeps. RSL_get_window_from_sweep calls RSL_get_window_from_ray for the number of rays. And RSL_get_window_from_ray gets window, between low_azim and hi_azim, from min_range to max_range. +RSL_get_window_from_radar calls RSL_get_window_from_volume for the number of volumes.
+RSL_get_window_from_volume calls RSL_get_window_from_sweep for the number of sweeps.
+RSL_get_window_from_sweep calls RSL_get_window_from_ray for the number of rays.
+RSL_get_window_from_ray gets window, between low_azim and hi_azim (degrees), from min_range to max_range (km).

Return value

@@ -27,8 +30,8 @@ Upon successful completion, a pointer to the appropriate structure is returned.

See also

-Example +rsl/examples/test_get_win.c
-

Author: John H. Merritt +

Author: John H. Merritt diff --git a/doc/RSL_new.html b/doc/RSL_new.html index cb15e9b..5b53589 100644 --- a/doc/RSL_new.html +++ b/doc/RSL_new.html @@ -39,10 +39,9 @@ Upon successful completion, a pointer to the appropriate structure is returned.

See also

-Example
RSL_free_volume, RSL_free_sweep, RSL_free_ray
RSL_copy_volume, RSL_copy_sweep, RSL_copy_ray
-

Author: John H. Merritt +

Author: John H. Merritt diff --git a/doc/RSL_read.html b/doc/RSL_read.html index e016545..6f3ac82 100644 --- a/doc/RSL_read.html +++ b/doc/RSL_read.html @@ -30,12 +30,11 @@ Upon successful completion, a pointer to the appropriate structure is returned.

See also

-Example
RSL_new_volume, RSL_new_sweep, RSL_new_ray,
RSL_write_volume, RSL_write_sweep, RSL_write_ray,
RSL_read_radar, RSL_write_radar,
File format.
-

Author: John H. Merritt +

Author: John H. Merritt diff --git a/dorade_to_radar.c b/dorade_to_radar.c index ff963a8..b66b9a8 100644 --- a/dorade_to_radar.c +++ b/dorade_to_radar.c @@ -25,6 +25,7 @@ #include #include #include +#include #define USE_RSL_VARS #include "rsl.h" #include "dorade.h" @@ -39,30 +40,112 @@ extern int radar_verbose_flag; int find_rsl_field_index(char *dorade_field_name) { /* - * Dorade: VE, DM, SW, DZ, ZDR, PHI, RHO, LDR, DX, CH, AH, CV, AV - * RSL: VR, DM, SW, DZ, ZD, PH, RH, LR, *DX, *CH, *AH, *CV, *AV. + * Dorade: VE, DM, SW, DBZ, ZDR, PHI, RHO, LDR, DX, CH, AH, CV, AV + * RSL: VR, DM, SW, DZ, ZD, PH, RH, LR, *DX, *CH, *AH, *CV, *AV. */ - if (strncasecmp(dorade_field_name, "ve", 2) == 0) return VR_INDEX; - if (strncasecmp(dorade_field_name, "dm", 2) == 0) return DM_INDEX; - if (strncasecmp(dorade_field_name, "sw", 2) == 0) return SW_INDEX; - if (strncasecmp(dorade_field_name, "dz", 2) == 0) return DZ_INDEX; + if (strncasecmp(dorade_field_name, "ve", 2) == 0) return VR_INDEX; + if (strncasecmp(dorade_field_name, "vr", 2) == 0) return VR_INDEX; + if (strncasecmp(dorade_field_name, "dm", 2) == 0) return DM_INDEX; + if (strncasecmp(dorade_field_name, "sw", 2) == 0) return SW_INDEX; + if (strncasecmp(dorade_field_name, "dbz", 3) == 0) return DZ_INDEX; if (strncasecmp(dorade_field_name, "zdr", 3) == 0) return ZD_INDEX; if (strncasecmp(dorade_field_name, "phi", 3) == 0) return PH_INDEX; if (strncasecmp(dorade_field_name, "rho", 3) == 0) return RH_INDEX; if (strncasecmp(dorade_field_name, "ldr", 3) == 0) return LR_INDEX; - if (strncasecmp(dorade_field_name, "dx", 2) == 0) return DX_INDEX; - if (strncasecmp(dorade_field_name, "ch", 2) == 0) return CH_INDEX; - if (strncasecmp(dorade_field_name, "ah", 2) == 0) return AH_INDEX; - if (strncasecmp(dorade_field_name, "cv", 2) == 0) return CV_INDEX; - if (strncasecmp(dorade_field_name, "av", 2) == 0) return AV_INDEX; + if (strncasecmp(dorade_field_name, "dx", 2) == 0) return DX_INDEX; + if (strncasecmp(dorade_field_name, "ch", 2) == 0) return CH_INDEX; + if (strncasecmp(dorade_field_name, "ah", 2) == 0) return AH_INDEX; + if (strncasecmp(dorade_field_name, "cv", 2) == 0) return CV_INDEX; + if (strncasecmp(dorade_field_name, "av", 2) == 0) return AV_INDEX; + if (strncasecmp(dorade_field_name, "vs", 2) == 0) return VS_INDEX; + if (strncasecmp(dorade_field_name, "vl", 2) == 0) return VL_INDEX; + if (strncasecmp(dorade_field_name, "vg", 2) == 0) return VG_INDEX; + if (strncasecmp(dorade_field_name, "vt", 2) == 0) return VT_INDEX; + if (strncasecmp(dorade_field_name, "ncp", 2) == 0) return NP_INDEX; - fprintf(stderr, "Unknown DORADE type <%s>\n", dorade_field_name); return -1; } +void prt_skipped_field_msg(char *dorade_field_name) +{ + char prtname[9]; + int i, already_printed; +#define MAXFIELDS 20 + static int nskipped = 0; + static char skipped_list[MAXFIELDS][9]; + + /* Make sure name is a properly formed string. */ + strncpy(prtname, dorade_field_name, 8); + prtname[8] = '\0'; + + /* We don't want to repeat messages for the same field, so keep a list of + * fields already printed. + */ + already_printed = 0; + i = 0; + while (!already_printed && i < nskipped) { + if (strncmp(prtname, skipped_list[i], 2) == 0) already_printed = 1; + i++; + } + if (!already_printed) { + fprintf(stderr, "Unknown DORADE parameter type <%s> -- skipping.\n", prtname); + strcpy(skipped_list[nskipped], prtname); + nskipped++; + } +} + +/* For date conversion routines. */ +static int daytab[2][13] = { + {0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334, 365}, + {0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335, 366} +}; + +/*************************************************************/ +/* */ +/* julian */ +/* */ +/*************************************************************/ +static int julian(int year, int mo, int day) +{ +/* Converts a calendar date (month, day, year) to a Julian date. + Returns: + Julian day. +*/ + int leap; + + leap = (year%4 == 0 && year%100 != 0) || year%400 == 0; + return(day + daytab[leap][mo-1]); +} + +/*************************************************************/ +/* */ +/* ymd */ +/* */ +/*************************************************************/ +static void ymd(int jday, int year, int *mm, int *dd) +{ + /* Input: jday, yyyy */ + /* Output: mm, dd */ + /* Copied from hdf_to_radar.c, written by Mike Kolander. */ + + int leap; + int i; + + leap = (year%4 == 0 && year%100 != 0) || year%400 == 0; + for (i=0; daytab[leap][i]nsensors, sizeof(Sensor_desc *)); for (i=0; insensors; i++) { - sd[i] = dorade_read_sensor(fp); + sd[i] = dorade_read_sensor(fp); } /* P R I N T */ if (radar_verbose_flag) { - for (i=0; insensors; i++) { - fprintf(stderr, "============ S E N S O R # %d =====================\n", i); - dorade_print_sensor(sd[i]); - } + for (i=0; insensors; i++) { + fprintf(stderr, "============ S E N S O R # %d =====================\n", i); + dorade_print_sensor(sd[i]); + } } /* R E A D sweeps. */ if (vd->nsensors > 1) { - fprintf(stderr, "RSL_dorade_to_radar: Unable to process for more than 1 sensor.\n"); - fprintf(stderr, "RSL_dorade_to_radar: Number of sensors is %d\n", vd->nsensors); - return NULL; + fprintf(stderr, "RSL_dorade_to_radar: Unable to process for more than 1 sensor.\n"); + fprintf(stderr, "RSL_dorade_to_radar: Number of sensors is %d\n", vd->nsensors); + return NULL; } /* Use sensor 0 for vitals. */ rd = sd[0]->radar_desc; + range_bin1 = sd[0]->cell_range_vector->range_cell[0]; + gate_size = sd[0]->cell_range_vector->range_cell[1] - range_bin1; radar = RSL_new_radar(MAX_RADAR_VOLUMES); radar->h.month = vd->month; @@ -159,7 +262,32 @@ Radar *RSL_dorade_to_radar(char *infile) radar->h.height = rd->altitude * 1000.0; radar->h.spulse = 0; /* FIXME */ radar->h.lpulse = 0; /* FIXME */ - + + year = vd->year; + month = vd->month; + day = vd->day; + jday_vol = julian(year, month, day); + + /* TODO: + Get any threshold values present in parameter descriptors. + + Note: The parameter descriptor contains an 8-character string which may + contain the, "Name of parameter upon which this parameter is thresholded". + According to documentation, if there is no thresholding parameter, the + string will be "NONE". In practice, this string has occasionally been + seen to contain only spaces for this condition. The corresponding + missing-data-flag for threshold may also vary, the nominal value being + -999, but has also been -32768. + + Pseudo code: + nparam = rd->nparam_desc; [<--useable code, i.e., not pseudo] + create string array for threshold parameter names, size nparams. + for each param: + strcpy thresh param name from param desc to threshold param array. + if string is all blanks, replace with 'NONE'. + endfor + */ + /* Begin volume code. */ /* We don't know how many sweeps per volume exist, until we read * the file. So allocate a large number of pointers, hope we don't @@ -169,21 +297,21 @@ Radar *RSL_dorade_to_radar(char *infile) */ if (radar_verbose_flag) - fprintf(stderr, "Number of parameters: %d\n", rd->nparam_desc); + fprintf(stderr, "Number of parameters: %d\n", rd->nparam_desc); /* All the parameters are together, however, their order within * the ray is not guarenteed. For instance, VE could appear after * DM. For this we'll keep a list of parameter names and perform * a (linear) search. The result will be an index into the RSL - * volume array (radar->v[i]). It is likely that the order will + * volume array (radar->v[i]). It is likely that the order will be * consistant within a file, therefore, we'll keep track of the index of * our previous parameter type and begin the search from there; the next * index should be a match. * * The dorade parameter names and the rsl mapping is: * - * Dorade: VE, DM, SW, DZ, ZDR, PHI, RHO, LDR, DX, CH, AH, CV, AV - * RSL: VR, DM, SW, DZ, ZD, PH, RH, LR, *DX, *CH, *AH, *CV, *AV. + * Dorade: VE, DM, SW, DBZ, ZDR, PHI, RHO, LDR, DX, CH, AH, CV, AV + * RSL: VR, DM, SW, DZ, ZD, PH, RH, LR, *DX, *CH, *AH, *CV, *AV. * * * means this is a new RSL name. */ @@ -191,56 +319,189 @@ Radar *RSL_dorade_to_radar(char *infile) #define DORADE_MAX_SWEEP 20 nsweep = 0; while((sr = dorade_read_sweep(fp, sd))) { - for(iray = 0; iray < sr->nrays; iray++) { - dray = sr->data_ray[iray]; - - /* Now, loop through the parameters and fill the rsl structures. */ - for (iparam = 0; iparam < dray->nparam; iparam++) { - pd = dray->parameter_data[iparam]; - iv = find_rsl_field_index(pd->name); - if (radar->v[iv] == NULL) { - radar->v[iv] = RSL_new_volume(DORADE_MAX_SWEEP); /* Expandable */ - } else if (nsweep >= radar->v[iv]->h.nsweeps) { - /* Must expand the number of sweeps. */ - /* Expand by another DORADE_MAX_SWEEP. */ - if (radar_verbose_flag) { - fprintf(stderr, "nsweeps (%d) exceeds radar->v[%d]->h.nsweeps (%d).\n", nsweep, iv, radar->v[iv]->h.nsweeps); - fprintf(stderr, "Increading it to %d sweeps\n", radar->v[iv]->h.nsweeps+DORADE_MAX_SWEEP); - } - new_volume = RSL_new_volume(radar->v[iv]->h.nsweeps+DORADE_MAX_SWEEP); - /* Look in uf_to_radar.c for 'copy_sweeps_into_volume' */ - new_volume = copy_sweeps_into_volume(new_volume, radar->v[iv]); - radar->v[iv] = new_volume; - } - - /* Allocate the ray and load the parameter data. */ - if ((sweep = radar->v[iv]->sweep[nsweep]) == NULL) { - sweep = radar->v[iv]->sweep[nsweep] = RSL_new_sweep(sr->s_info->nrays); - sweep->h.sweep_num = sr->s_info->sweep_num; - sweep->h.elev = sr->s_info->fixed_angle; + for(iray = 0; iray < sr->nrays; iray++) { + dray = sr->data_ray[iray]; + + /* Now, loop through the parameters and fill the rsl structures. */ + for (iparam = 0; iparam < dray->nparam; iparam++) { + pd = dray->parameter_data[iparam]; + iv = find_rsl_field_index(pd->name); + if (iv < 0) { + prt_skipped_field_msg(pd->name); + continue; + } + if (radar->v[iv] == NULL) { + radar->v[iv] = RSL_new_volume(DORADE_MAX_SWEEP); /* Expandable */ + } else if (nsweep >= radar->v[iv]->h.nsweeps) { + /* Must expand the number of sweeps. */ + /* Expand by another DORADE_MAX_SWEEP. */ + if (radar_verbose_flag) { + fprintf(stderr, "nsweeps (%d) exceeds radar->v[%d]->h.nsweeps (%d)." + "\n", nsweep, iv, radar->v[iv]->h.nsweeps); + fprintf(stderr, "Increasing it to %d sweeps\n", + radar->v[iv]->h.nsweeps+DORADE_MAX_SWEEP); + } + new_volume = RSL_new_volume(radar->v[iv]->h.nsweeps+DORADE_MAX_SWEEP); + /* Look in uf_to_radar.c for 'copy_sweeps_into_volume' */ + new_volume = copy_sweeps_into_volume(new_volume, radar->v[iv]); + radar->v[iv] = new_volume; + } + /* Assign f and invf + switch (iv) . . . + * Or just use RSL_f_list[iv] and RSL_invf_list[iv] as in sweep.h below. + */ + + /* Allocate the ray and load the parameter data. */ + if ((sweep = radar->v[iv]->sweep[nsweep]) == NULL) { + sweep = radar->v[iv]->sweep[nsweep] = RSL_new_sweep(sr->s_info->nrays); + sweep->h.sweep_num = sr->s_info->sweep_num; + sweep->h.elev = sr->s_info->fixed_angle; sweep->h.beam_width = rd->horizontal_beam_width; - sweep->h.vert_half_bw = radar->v[iv]->sweep[nsweep]->h.beam_width / 2.0; - sweep->h.horz_half_bw = rd->horizontal_beam_width / 2.0; - sweep->h.f = RSL_f_list[iv]; - sweep->h.invf = RSL_invf_list[iv]; - } - - - if ((ray = sweep->ray[iray]) == NULL) { - if (radar_verbose_flag) - fprintf(stderr, "Allocating %d bins for ray %d\n", dray->data_len[iparam], iray); - ray = sweep->ray[iray] = RSL_new_ray(dray->data_len[iparam]); - } - - /* Copy the ray data into the RSL ray. */ + sweep->h.vert_half_bw = rd->vertical_beam_width / 2.0; + /* + sweep->h.vert_half_bw = radar->v[iv]->sweep[nsweep]->h.beam_width / 2.0; + */ + sweep->h.horz_half_bw = rd->horizontal_beam_width / 2.0; + sweep->h.f = RSL_f_list[iv]; + sweep->h.invf = RSL_invf_list[iv]; + } + + + data_len = dray->data_len[iparam]; + word_size = dray->word_size[iparam]; + if (word_size != 2 && word_size != 4) { + fprintf(stderr,"RSL_dorade_to_radar: dray->word_size[%d] = %d\n", + iparam, dray->word_size[iparam]); + fprintf(stderr,"Expected value is 2 or 4.\n"); + return NULL; + } + nbins = data_len / word_size; + + if ((ray = sweep->ray[iray]) == NULL) { + if (radar_verbose_flag) + fprintf(stderr, "Allocating %d bins for ray %d\n", + dray->data_len[iparam], iray); + ray = sweep->ray[iray] = RSL_new_ray(nbins); + } + + /* Load ray header. */ + + ray_info = dray->ray_info; + jday = ray_info->jday; + if (jday != jday_vol) { + if (jday > 0 && jday < 367) { + /* Recompute year month day */ + if (jday < jday_vol) year++; + ymd(jday, year, &month, &day); + jday_vol = jday; + } + else { /* Invalid jday */ + } + } + ray->h.year = year; + ray->h.month = month; + ray->h.day = day; + ray->h.hour = ray_info->hour; + ray->h.minute = ray_info->minute; + ray->h.sec = ray_info->second + ray_info->msec / 1000.; + ray->h.azimuth = ray_info->azimuth; + ray->h.ray_num = iray + 1; + ray->h.elev = ray_info->elevation; + ray->h.elev_num = ray_info->sweep_num; + ray->h.unam_rng = rd->unambiguous_range; + ray->h.nyq_vel = rd->unambiguous_velocity; + ray->h.azim_rate = ray_info->scan_rate; + /* TODO + ray->h.fix_angle = ; + */ + ray->h.range_bin1 = range_bin1; + ray->h.gate_size = gate_size; + platform_info = dray->platform_info; + ray->h.pitch = platform_info->pitch; + ray->h.roll = platform_info->roll; + ray->h.heading = platform_info->heading; + ray->h.pitch_rate = platform_info->pitch_rate; + ray->h.heading_rate = platform_info->heading_rate; + ray->h.lat = platform_info->latitude; + ray->h.lon = platform_info->longitude; + ray->h.alt = platform_info->altitude; + ray->h.vel_east = platform_info->ew_speed; + ray->h.vel_north = platform_info->ns_speed; + ray->h.vel_up = platform_info->v_speed; + /* TODO: Does DORADE have Doppler velocity res? + A: No. + ray->h.vel_res = N/A + */ + ray->h.beam_width = rd->horizontal_beam_width; + /* TODO + ray->h.frequency = @Radar Descriptor + Get parm_desc.xmit_freq + Find first set-bit, use frequency from corresponding index in rad_desc frequenies. + + ray->h.pulse_count = Is it number samples used? @Param desc. Probably not. + + * The following are DONE * + ray->h.pulse_width = @Parameter Descriptor--not same--it's in meters-- + we use micro-sec. They're using pulse length. + pulse width = pulse length/c + ray->h.wavelength = Can compute using Nyquist velocity and PRF or PRT: + wl = 4*vmax/PRF or wl = 4*vmax*PRT + ray->h.prf = Can compute if nothing else: prf = c/(2*Rmax) + */ + /* DORADE is actually using pulse length. Convert to pulse width. */ + ray->h.pulse_width = (sd[0]->p_desc[iparam]->pulse_width/ + RSL_SPEED_OF_LIGHT) * 1e6; + prf = RSL_SPEED_OF_LIGHT/(2.*rd->unambiguous_range*1000.); + ray->h.prf = prf; + ray->h.wavelength = (4.*rd->unambiguous_velocity)/prf; + ray->h.nbins = nbins; + ray->h.f = RSL_f_list[iv]; + ray->h.invf = RSL_invf_list[iv]; + + /* Copy the ray data into the RSL ray. */ /* .... fill here .... */ + + /* Assign pointer to data. + * Get datum using word size and proper cast. + * Convert and store in rsl ray->range. + * Increment pointer to data based on word size. + */ + + ptr2bytes = (short *) pd->data; + ptr4bytes = (int *) pd->data; + scale = sd[0]->p_desc[iparam]->scale_factor; + offset = sd[0]->p_desc[iparam]->offset_factor; + missing_data_flag = sd[0]->p_desc[iparam]->missing_data_flag; + + for (i=0; irange[i] = ray->h.invf(value); } - nsweep++; - if (radar_verbose_flag) fprintf(stderr, "______NEW SWEEP__<%d>____\n", nsweep); - /* Save for loading into volume structure. */ - dorade_free_sweep(sr); + + if (iray == 0) { + radar->v[iv]->h.nsweeps = nsweep + 1; + radar->v[iv]->h.f = RSL_f_list[iv]; + radar->v[iv]->h.invf = RSL_invf_list[iv]; + } + } + } + nsweep++; + if (radar_verbose_flag) fprintf(stderr, "______NEW SWEEP__<%d>____\n", nsweep); + /* Save for loading into volume structure. */ + dorade_free_sweep(sr); } /* The following avoids a broken pipe message, since a VOLD at the end @@ -252,4 +513,3 @@ Radar *RSL_dorade_to_radar(char *infile) return radar; } - diff --git a/examples/Makefile.am b/examples/Makefile.am index 3fdf629..7bcf745 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -5,14 +5,14 @@ AUTOMAKE_OPTIONS = foreign INCLUDES = -I$(prefix)/include LOCAL_LIB = ../.libs/librsl.a LDADD = @LIBS@ $(LOCAL_LIB) -bin_PROGRAMS = any_to_gif any_to_uf +bin_PROGRAMS = any_to_gif any_to_uf qlook any_to_gif_LDFLAGS = -static any_to_uf_LDFLAGS = -static # Additional program to build but not install noinst_PROGRAMS = any_to_ppm any_to_ufgz bscan \ cappi_image dorade_main killer_sweep \ kwaj_subtract_one_day lassen_to_gif print_hash_table \ - print_header_info qlook sector test_get_win \ + print_header_info sector test_get_win \ wsr88d_to_gif wsr_hist_uf_test diff --git a/examples/qlook.c b/examples/qlook.c index dc767f6..5fa6c17 100644 --- a/examples/qlook.c +++ b/examples/qlook.c @@ -5,25 +5,42 @@ #include #include #include +#include +#define USE_RSL_VARS #include "rsl.h" #define ZDR_WIDTH 10 int verbose; -void rebin_ray(Ray *r, int width); -void rebin_sweep(Sweep *s, int width); -void rebin_volume(Volume *v, int width); +void RSL_rebin_ray(Ray *r, int width); +void RSL_rebin_sweep(Sweep *s, int width); +void RSL_rebin_volume(Volume *v, int width); void adjust_gate_size(Radar *radar, float gate_size_adjustment); void process_args(int argc, char **argv, char *in_file, int *verbose, - char *site_id, char *tape_id, - int *qc_reflectivity, int *total_reflectivity, - int *differential_reflectivity, - int *velocity, int *spectral_width, - int *make_gif, int *make_pgm, int *make_bscan, int *make_uf, - int *num_sweeps, float *dbz_offset, - int *xdim, int *ydim, float *range, - float *gate_size_adjustment, int *print_azim); + char *site_id, char *tape_id, + int *qc_reflectivity, int *total_reflectivity, + int *differential_reflectivity, + int *velocity, int *spectral_width, + int *make_gif, int *make_pgm, int *make_bscan, int *make_uf, + int *num_sweeps, float *dbz_offset, + int *xdim, int *ydim, float *range, + float *gate_size_adjustment, int *print_azim, + char *gifdir, char *pgmdir, char *ufdir); + +void make_pathname(char *filename, char *dir, char *pathname) +{ + /* Make pathname by combining directory name, if given, with filename. */ + + if (*dir == NULL) { + strcpy(pathname, filename); + } + else { + strcpy(pathname, dir); + strcat(pathname, "/"); + strcat(pathname, filename); + } +} /***********************************************************************/ /* This program uses the NASA TRMM Office Radar Software Library (RSL) */ @@ -32,77 +49,90 @@ void process_args(int argc, char **argv, char *in_file, int *verbose, /***********************************************************************/ main(int argc, char **argv) { - Radar *radar; - Volume *dz_volume, *vr_volume; - Volume *dr_volume, *zt_volume; - Volume *sw_volume, *qc_volume, *volume; - - Sweep *sweep; - Sweep *dz_sweep, *qc_sweep; - Sweep *dr_sweep, *zt_sweep; - Sweep *vr_sweep, *sw_sweep; - - Ray *ray; - int reflectivity, qc_reflectivity, nvolumes; - int differential_reflectivity, total_reflectivity; - int velocity, spectral_width; - int make_catalog, make_gif, make_pgm; - int make_uf, make_bscan; - int xdim, ydim, num_sweeps; - int i,j,k,l,n, verbose, kk; - int month, day, year, hour, min; - int ncbins, width; - int print_azim; - float scale, dbz_offset; - float latitude, longitude; - - float sec; - float maxr; - float nyquist, max_range, gate_size_adjustment; - - char tape_id[100]; - char in_file[100], site_id[100]; - char filename[100], outfile[100], nexfile[100]; - char command[100], file_prefix[100], file_suffix[3]; - char dir_string[100], red[120], grn[120], blu[120]; - char time_string[14], site_string[10],img_base[20]; - FILE *fp; + Radar *radar; + Volume *dz_volume, *vr_volume; + Volume *dr_volume, *zt_volume; + Volume *sw_volume, *qc_volume, *volume; + + Sweep *sweep; + Sweep *dz_sweep, *qc_sweep; + Sweep *dr_sweep, *zt_sweep; + Sweep *vr_sweep, *sw_sweep; + + Ray *ray; + int reflectivity, qc_reflectivity, nvolumes; + int differential_reflectivity, total_reflectivity; + int velocity, spectral_width; + int make_catalog, make_gif, make_pgm; + int make_uf, make_bscan; + int xdim, ydim, num_sweeps; + int i,j,k,l,n, verbose, kk; + int month, day, year, hour, min; + int ncbins, width; + int print_azim; + float scale, dbz_offset; + float latitude, longitude; + + float sec; + float maxr; + float nyquist, max_range, gate_size_adjustment; + + char tape_id[100]; + char in_file[100], site_id[100]; + char filename[100], outfile[100], nexfile[100]; + char command[100], file_prefix[100], file_suffix[3]; + char dir_string[100], red[120], grn[120], blu[120]; + char time_string[14], site_string[10],img_base[20]; + char pathname[256], gifdir[100], pgmdir[100], ufdir[100]; + char *inpath; + FILE *fp; /* Set default values */ - strcpy(site_id, "KMLB"); - strcpy(tape_id, "WSR88D"); - - verbose = FALSE; - reflectivity = TRUE; - qc_reflectivity = FALSE; - total_reflectivity = FALSE; - differential_reflectivity = FALSE; - velocity = FALSE; - spectral_width = FALSE; - make_gif = FALSE; - make_pgm = FALSE; - make_catalog = FALSE; - make_uf = FALSE; - print_azim = FALSE; - num_sweeps = 1; - xdim = ydim = 400; - maxr = 200.; - dbz_offset = 0.0; - gate_size_adjustment = 1.0; - + /*strcpy(site_id, "KMLB");*/ + strcpy(tape_id, "WSR88D"); + + verbose = FALSE; + reflectivity = TRUE; + qc_reflectivity = FALSE; + total_reflectivity = FALSE; + differential_reflectivity = FALSE; + velocity = FALSE; + spectral_width = FALSE; + make_gif = FALSE; + make_pgm = FALSE; + make_catalog = FALSE; + make_uf = FALSE; + print_azim = FALSE; + num_sweeps = 1; + xdim = ydim = 400; + maxr = 200.; + dbz_offset = 0.0; + gate_size_adjustment = 1.0; + *gifdir = *pgmdir = *ufdir = '\0'; + *site_id = '\0'; + /* Process command_line arguments */ - process_args(argc, argv, in_file, &verbose, - site_id, tape_id, - &qc_reflectivity, &total_reflectivity, - &differential_reflectivity, - &velocity, &spectral_width, - &make_gif, &make_pgm, &make_bscan, &make_uf, - &num_sweeps, &dbz_offset, - &xdim, &ydim, &maxr, &gate_size_adjustment, - &print_azim); - + process_args(argc, argv, in_file, &verbose, + site_id, tape_id, + &qc_reflectivity, &total_reflectivity, + &differential_reflectivity, + &velocity, &spectral_width, + &make_gif, &make_pgm, &make_bscan, &make_uf, + &num_sweeps, &dbz_offset, + &xdim, &ydim, &maxr, &gate_size_adjustment, + &print_azim, gifdir, pgmdir, ufdir); + +/* If site_id empty, use first 4 characters of file name. */ + + if (*site_id == 0) { + inpath = strdup(in_file); + strncpy(site_id, basename(inpath), 4); + site_id[4] = '\0'; + if (strcmp(site_id, "KWA0") == 0) strcpy(site_id, "KWAJ"); + } + /* Be a chatty Kathy? */ if(verbose) @@ -112,128 +142,131 @@ main(int argc, char **argv) { /* Read data into radar */ - if(verbose) printf("Calling any_format_to_radar\n"); - radar = RSL_anyformat_to_radar(in_file, site_id); - if(verbose) printf("Called any_format_to_radar\n"); - if(radar==NULL) { - if (verbose) - printf("No radar loaded, bye\n"); - exit(-1); - } + if(verbose) printf("Calling any_format_to_radar\n"); + radar = RSL_anyformat_to_radar(in_file, site_id); + if(verbose) printf("Called any_format_to_radar\n"); + if(radar==NULL) { + if (verbose) + printf("No radar loaded, bye\n"); + exit(-1); + } /* Print command line parameters */ - if(verbose) { - printf("Site ID = %s\n",site_id); - printf("Tape ID = %s\n",tape_id); - printf("Do reflectivity = %d\n",reflectivity); - printf("Do qc_reflectivity = %d\n",qc_reflectivity); - printf("Do differential_reflectivity = %d\n", - differential_reflectivity); - printf("Do total_reflectivity = %d\n", - total_reflectivity); - printf("Do qc_reflectivity = %d\n",qc_reflectivity); - printf("Do velocity = %d\n",velocity); - printf("Do spectral_width = %d\n",spectral_width); - printf("Make gif = %d\n",make_gif); - printf("Make pgm = %d\n",make_pgm); - printf("Make UF file = %d\n",make_uf); - printf("dBZ Offset = %.2f\n",dbz_offset); - printf("Gate Size Adjust = %.2f\n",gate_size_adjustment); - printf("Print Azimuths = %d\n",print_azim); - } + if(verbose) { + printf("Site ID = %s\n",site_id); + printf("Tape ID = %s\n",tape_id); + printf("Do reflectivity = %d\n",reflectivity); + printf("Do qc_reflectivity = %d\n",qc_reflectivity); + printf("Do differential_reflectivity = %d\n", + differential_reflectivity); + printf("Do total_reflectivity = %d\n", + total_reflectivity); + printf("Do qc_reflectivity = %d\n",qc_reflectivity); + printf("Do velocity = %d\n",velocity); + printf("Do spectral_width = %d\n",spectral_width); + printf("Make gif = %d\n",make_gif); + printf("Make pgm = %d\n",make_pgm); + printf("Make UF file = %d\n",make_uf); + if (ufdir != NULL) printf("UF output dir = %s\n",ufdir); + if (gifdir != NULL) printf("GIF output dir = %s\n",gifdir); + if (pgmdir != NULL) printf("PGM output dir = %s\n",pgmdir); + printf("dBZ Offset = %.2f\n",dbz_offset); + printf("Gate Size Adjust = %.2f\n",gate_size_adjustment); + printf("Print Azimuths = %d\n",print_azim); + } /* If Gate Size Adjustment is not unity, then we must change the following: old_gs = radar->v[i]->sweep[sweepIndex]->ray[rayIndex=]->h.gate_size - radar->v[i]->sweep[sweepIndex]->ray[rayIndex=]->h.gate_size = - old_gs*gate_size_adjustment + radar->v[i]->sweep[sweepIndex]->ray[rayIndex=]->h.gate_size = + old_gs*gate_size_adjustment Here are some comments from Sandra Yuter on the necessity of this fix. > I dug into the revelant code and it looks like we can do a relatively - > simple workaround for the SIGMET raw product file range bin size - > errors for the RHB data pulses widths of 0.5 usec and 2.0 usec as follows. - > - > Since we are all converting from sigmet to UF I suggest we resize - > the range bin size values in the ray headers in the qlook step - > where the sigmet to UF conversion occurs. - > - > The resize requires only 1 additional line of code (I have included - > a few others for context) in qlook.c - > - > rangeToFirstGate = 0.001 * - > radar->v[i]->sweep[sweepIndex]->ray[rayIndex]->h.range_bin1; - > gateSize = 0.001 * - > radar->v[i]->sweep[sweepIndex]->ray[rayIndex]->h.gate_size; - > radar->v[i]->sweep[sweepIndex]->ray[rayIndex]->h.gate_size= - > gateSize*0.6*1000; - > - > I have used 0.6 adjustment factor since that is 75/125 which corresponds - > to the error in my 0.5 usec data, for the SUR scans, this adjustment - > factor is 133.33/125 or 1.067. - - The following is from Joe Holmes from SIGMET - - > - > I think you have experienced a problem with the RVP7 range resolution - > configuration. Both in IRIS and in the RVP7 you manually type in - > the range resolution. The RVP7 allows a separate resolution for - > each pulsewidth, while IRIS only allows 1 value. There is no feedback - > if these values are not typed in the same. Based on setup information - > we have here from the RH Brown from Oct 23, 1998, you had the following - > configuration: - > - > RVP7: - > 0 0.50 usec 75.0m - > 1 0.80 usec 125.0m - > 2 1.39 usec 125.0m - > 3 2.00 usec 133.3m - > - > IRIS: 125.0 meters - > - > I think the error at PW#0 was corrected a while back, but - > the error in PW#3 was never corrected. Next time someone is - > at the ship, they should check this, fix the long pulse, and remake - > the bandpass filter for the long pulse. - > - > In the short term, you can correct the error by taking all your - > long pulse data and changing the header to correctly document the - > range resolution. We have a program to do this, it is called "change_raw". - > The source is on any IRIS system, which was installed with the - > source, headers, and objects turned on. It is in the - > ${IRIS_ROOT}utils/examples directory. We can supply you with - > a compiled version of this program, if you want. Available platforms - > are Linux, HP-UX, and IRIX. + > simple workaround for the SIGMET raw product file range bin size + > errors for the RHB data pulses widths of 0.5 usec and 2.0 usec as follows. + > + > Since we are all converting from sigmet to UF I suggest we resize + > the range bin size values in the ray headers in the qlook step + > where the sigmet to UF conversion occurs. + > + > The resize requires only 1 additional line of code (I have included + > a few others for context) in qlook.c + > + > rangeToFirstGate = 0.001 * + > radar->v[i]->sweep[sweepIndex]->ray[rayIndex]->h.range_bin1; + > gateSize = 0.001 * + > radar->v[i]->sweep[sweepIndex]->ray[rayIndex]->h.gate_size; + > radar->v[i]->sweep[sweepIndex]->ray[rayIndex]->h.gate_size= + > gateSize*0.6*1000; + > + > I have used 0.6 adjustment factor since that is 75/125 which corresponds + > to the error in my 0.5 usec data, for the SUR scans, this adjustment + > factor is 133.33/125 or 1.067. + + The following is from Joe Holmes from SIGMET + + > + > I think you have experienced a problem with the RVP7 range resolution + > configuration. Both in IRIS and in the RVP7 you manually type in + > the range resolution. The RVP7 allows a separate resolution for + > each pulsewidth, while IRIS only allows 1 value. There is no feedback + > if these values are not typed in the same. Based on setup information + > we have here from the RH Brown from Oct 23, 1998, you had the following + > configuration: + > + > RVP7: + > 0 0.50 usec 75.0m + > 1 0.80 usec 125.0m + > 2 1.39 usec 125.0m + > 3 2.00 usec 133.3m + > + > IRIS: 125.0 meters + > + > I think the error at PW#0 was corrected a while back, but + > the error in PW#3 was never corrected. Next time someone is + > at the ship, they should check this, fix the long pulse, and remake + > the bandpass filter for the long pulse. + > + > In the short term, you can correct the error by taking all your + > long pulse data and changing the header to correctly document the + > range resolution. We have a program to do this, it is called "change_raw". + > The source is on any IRIS system, which was installed with the + > source, headers, and objects turned on. It is in the + > ${IRIS_ROOT}utils/examples directory. We can supply you with + > a compiled version of this program, if you want. Available platforms + > are Linux, HP-UX, and IRIX. */ - - if(gate_size_adjustment != 1.0) { - printf("Adjusting Gate Size by Factor: %.3f\n",gate_size_adjustment); - adjust_gate_size(radar, gate_size_adjustment); - } - + + if(gate_size_adjustment != 1.0) { + printf("Adjusting Gate Size by Factor: %.3f\n",gate_size_adjustment); + adjust_gate_size(radar, gate_size_adjustment); + } + /* Create the filename prefix. Consists of the Site ID, and time string (YYMMDD_hhmm). The file suffix is like MIME type (e.g, .gif, .pgm, .uf, etc.) */ - sprintf(time_string,"%4d/%2.2d%2.2d %2.2d:%2.2d UTC", - radar->h.year, radar->h.month, radar->h.day, - radar->h.hour, radar->h.minute); + sprintf(time_string,"%4d/%2.2d%2.2d %2.2d:%2.2d UTC", + radar->h.year, radar->h.month, radar->h.day, + radar->h.hour, radar->h.minute); /* Determine the location (lat/lon) of the radar. */ - latitude = radar->h.latd + radar->h.latm/60. + radar->h.lats/3600.; - longitude = radar->h.lond + radar->h.lonm/60. + radar->h.lons/3600.; + latitude = radar->h.latd + radar->h.latm/60. + radar->h.lats/3600.; + longitude = radar->h.lond + radar->h.lonm/60. + radar->h.lons/3600.; - printf("%s %s %s %.6f %.6f \n", - in_file, radar->h.radar_name, time_string, longitude, latitude); + printf("%s %s %s %.6f %.6f \n", + in_file, radar->h.radar_name, time_string, longitude, latitude); - sprintf(time_string,"%4d_%2.2d%2.2d_%2.2d%2.2d", - radar->h.year, radar->h.month, radar->h.day, - radar->h.hour, radar->h.minute); + sprintf(time_string,"%4d_%2.2d%2.2d_%2.2d%2.2d", + radar->h.year, radar->h.month, radar->h.day, + radar->h.hour, radar->h.minute); /* Print the radar/volume info. @@ -317,153 +350,153 @@ main(int argc, char **argv) { */ - if(verbose) { - for(i=0; i< radar->h.nvolumes; i++) { - if(radar->v[i] != NULL) { - printf("Vol[%2.2d] has %d sweeps\n",i,radar->v[i]->h.nsweeps); - } else { - printf("Vol[%2.2d] == NULL\n",i); - } - } - printf("--------------------------------------------\n"); - printf("Number of volumes in radar: %d\n",radar->h.nvolumes); - } - - /* DZ_INDEX */ - if(radar->v[DZ_INDEX] == NULL) { - printf("DZ_INDEX == NULL\n"); - reflectivity = FALSE; - } else { - dz_volume = radar->v[DZ_INDEX]; - if(verbose) printf("Number of sweeps in dz_volume = %d\n", - dz_volume->h.nsweeps); - } - - /* CZ_INDEX */ - if(radar->v[CZ_INDEX] == NULL) { - if(verbose) printf("CZ_INDEX == NULL\n"); - qc_reflectivity = FALSE; - } else { - qc_volume = radar->v[CZ_INDEX]; - if(verbose) printf("Number of sweeps in qc_volume = %d\n", - qc_volume->h.nsweeps); - } - - /* ZT_INDEX */ - if(radar->v[ZT_INDEX] == NULL) { - if(verbose) printf("ZT_INDEX == NULL\n"); - total_reflectivity = FALSE; - } else { - zt_volume = radar->v[ZT_INDEX]; - if(verbose) printf("Number of sweeps in zt_volume = %d\n", - zt_volume->h.nsweeps); - } - /* ZD_INDEX */ - if(radar->v[ZD_INDEX] == NULL) { - if(verbose) printf("ZD_INDEX == NULL\n"); - differential_reflectivity = FALSE; - } else { - dr_volume = radar->v[ZD_INDEX]; - if(verbose) printf("Number of sweeps in dr_volume = %d\n", - dr_volume->h.nsweeps); - } - /* VR_INDEX */ - if(radar->v[VR_INDEX] == NULL) { - if(verbose) printf("VR_INDEX == NULL\n"); - velocity = FALSE; - } else { - vr_volume = radar->v[VR_INDEX]; - if(verbose) printf("Number of sweeps in vr_volume = %d\n", - vr_volume->h.nsweeps); - } - - /* SW_INDEX */ - if(radar->v[SW_INDEX] == NULL) { - if(verbose) printf("SW_INDEX == NULL\n"); - spectral_width = FALSE; - } else { - sw_volume = radar->v[SW_INDEX]; - if(verbose) printf("Number of sweeps in sw_volume = %d\n", - sw_volume->h.nsweeps); - } - if(verbose) printf("--------------------------------------------\n"); + if(verbose) { + for(i=0; i< radar->h.nvolumes; i++) { + if(radar->v[i] != NULL) { + printf("Vol[%2.2d] has %d sweeps\n",i,radar->v[i]->h.nsweeps); + } else { + printf("Vol[%2.2d] == NULL\n",i); + } + } + printf("--------------------------------------------\n"); + printf("Number of volumes in radar: %d\n",radar->h.nvolumes); + } + + /* DZ_INDEX */ + if(radar->v[DZ_INDEX] == NULL) { + printf("DZ_INDEX == NULL\n"); + reflectivity = FALSE; + } else { + dz_volume = radar->v[DZ_INDEX]; + if(verbose) printf("Number of sweeps in dz_volume = %d\n", + dz_volume->h.nsweeps); + } + + /* CZ_INDEX */ + if(radar->v[CZ_INDEX] == NULL) { + if(verbose) printf("CZ_INDEX == NULL\n"); + qc_reflectivity = FALSE; + } else { + qc_volume = radar->v[CZ_INDEX]; + if(verbose) printf("Number of sweeps in qc_volume = %d\n", + qc_volume->h.nsweeps); + } + + /* ZT_INDEX */ + if(radar->v[ZT_INDEX] == NULL) { + if(verbose) printf("ZT_INDEX == NULL\n"); + total_reflectivity = FALSE; + } else { + zt_volume = radar->v[ZT_INDEX]; + if(verbose) printf("Number of sweeps in zt_volume = %d\n", + zt_volume->h.nsweeps); + } + /* DR_INDEX */ + if(radar->v[DR_INDEX] == NULL) { + if(verbose) printf("DR_INDEX == NULL\n"); + differential_reflectivity = FALSE; + } else { + dr_volume = radar->v[DR_INDEX]; + if(verbose) printf("Number of sweeps in dr_volume = %d\n", + dr_volume->h.nsweeps); + } + /* VR_INDEX */ + if(radar->v[VR_INDEX] == NULL) { + if(verbose) printf("VR_INDEX == NULL\n"); + velocity = FALSE; + } else { + vr_volume = radar->v[VR_INDEX]; + if(verbose) printf("Number of sweeps in vr_volume = %d\n", + vr_volume->h.nsweeps); + } + + /* SW_INDEX */ + if(radar->v[SW_INDEX] == NULL) { + if(verbose) printf("SW_INDEX == NULL\n"); + spectral_width = FALSE; + } else { + sw_volume = radar->v[SW_INDEX]; + if(verbose) printf("Number of sweeps in sw_volume = %d\n", + sw_volume->h.nsweeps); + } + if(verbose) printf("--------------------------------------------\n"); /* Print out the elevation angles */ - if(verbose) { - if(reflectivity) { - printf("Reflectivity Tilts\n"); - printf("----------------------\n"); - if(dz_volume != NULL) { - for(i=0; ih.nsweeps; i++) { - sweep = dz_volume->sweep[i]; - if(sweep == NULL) { - printf("sweep[%d]==NULL\n",i); - continue; - } - printf("Tilt %2d, Elev=%4.1f\n",i,sweep->h.elev); - } - printf("----------------------\n"); - } - } - } - /* - Print out the values of the rays in each sweep requsted - */ - - if(print_azim) { - printf("Ray angles\n"); - if(reflectivity) { - for(i=0; ih.nsweeps; i++) { - if(dz_volume->sweep[i] != NULL) sweep = dz_volume->sweep[i]; - printf("Elevation angle: %f\n",sweep->h.elev); - printf("Number of rays in sweep[%d] = %d\n",i,sweep->h.nrays); - - for(j=0; jh.nrays-1; j++) { - if(sweep->ray[j] != NULL) { - ray = sweep->ray[j]; - printf("%d: %7.2f\n ",j,sweep->ray[j]->h.azimuth); - } - } - } - } - } + if(verbose) { + if(reflectivity) { + printf("Reflectivity Tilts\n"); + printf("----------------------\n"); + if(dz_volume != NULL) { + for(i=0; ih.nsweeps; i++) { + sweep = dz_volume->sweep[i]; + if(sweep == NULL) { + printf("sweep[%d]==NULL\n",i); + continue; + } + printf("Tilt %2d, Elev=%4.1f\n",i,sweep->h.elev); + } + printf("----------------------\n"); + } + } + } + /* + Print out the values of the rays in each sweep requsted + */ + + if(print_azim) { + printf("Ray angles\n"); + if(reflectivity) { + for(i=0; ih.nsweeps; i++) { + if(dz_volume->sweep[i] != NULL) sweep = dz_volume->sweep[i]; + printf("Elevation angle: %f\n",sweep->h.elev); + printf("Number of rays in sweep[%d] = %d\n",i,sweep->h.nrays); + + for(j=0; jh.nrays-1; j++) { + if(sweep->ray[j] != NULL) { + ray = sweep->ray[j]; + printf("%d: %7.2f\n ",j,sweep->ray[j]->h.azimuth); + } + } + } + } + } /* Write out some volume statistics */ - if(verbose) { - printf("******************* Volume Statistics *******************\n"); - if(reflectivity) { - sweep = RSL_get_first_sweep_of_volume(dz_volume); - if(sweep != NULL) { - printf("Number rays = %d\n", sweep->h.nrays); - printf("Beam width = %.2f deg\n", sweep->h.beam_width); - ray = RSL_get_first_ray_of_sweep(sweep); - if(ray!= NULL) { - max_range = ray->h.range_bin1/1000.0 + - ray->h.nbins*ray->h.gate_size/1000.0; - printf("Number of bins = %d\n",ray->h.nbins); - printf("Max range = %.1f km\n", max_range); - printf("Range to first bin = %d m\n",ray->h.range_bin1); - printf("Gate size = %d m\n",ray->h.gate_size); - printf("Pulse Rep. Freq. = %d Hz\n",ray->h.prf); - printf("Pulse width = %.2f us\n",ray->h.pulse_width); - printf("Wavelength = %.2f m\n",ray->h.wavelength); - printf("Frequency = %.2f \n",ray->h.frequency); - } - } - } - } + if(verbose) { + printf("******************* Volume Statistics *******************\n"); + if(reflectivity) { + sweep = RSL_get_first_sweep_of_volume(dz_volume); + if(sweep != NULL) { + printf("Number rays = %d\n", sweep->h.nrays); + printf("Beam width = %.2f deg\n", sweep->h.beam_width); + ray = RSL_get_first_ray_of_sweep(sweep); + if(ray!= NULL) { + max_range = ray->h.range_bin1/1000.0 + + ray->h.nbins*ray->h.gate_size/1000.0; + printf("Number of bins = %d\n",ray->h.nbins); + printf("Max range = %.1f km\n", max_range); + printf("Range to first bin = %d m\n",ray->h.range_bin1); + printf("Gate size = %d m\n",ray->h.gate_size); + printf("Pulse Rep. Freq. = %d Hz\n",ray->h.prf); + printf("Pulse width = %.2f us\n",ray->h.pulse_width); + printf("Wavelength = %.2f m\n",ray->h.wavelength); + printf("Frequency = %.2f \n",ray->h.frequency); + } + } + } + } /* - Add a dBZ offset if requested. The offset can be positive or negative. - if(dbz_offset != 0.0) { - printf("Adding dbz_offset to dz_volume: %.2f\n", dbz_offset); - RSL_add_dbzoffset_to_volume(dz_volume, dbz_offset); - printf("Added dbz_offset to dz_volume: %.2f\n", dbz_offset); - exit(0); - } + Add a dBZ offset if requested. The offset can be positive or negative. + if(dbz_offset != 0.0) { + printf("Adding dbz_offset to dz_volume: %.2f\n", dbz_offset); + RSL_add_dbzoffset_to_volume(dz_volume, dbz_offset); + printf("Added dbz_offset to dz_volume: %.2f\n", dbz_offset); + exit(0); + } */ /* @@ -473,215 +506,235 @@ main(int argc, char **argv) { */ - /* CZ_INDEX */ - if(qc_reflectivity) { - if(verbose) printf("Loading refl colortable...\n"); - RSL_load_refl_color_table(); - for(i=0; isweep[i]; - if(sweep == NULL) { - printf("sweep[%d]==NULL\n",i); - continue; - } - if(make_pgm) { - sprintf(file_suffix,"pgm"); - sprintf(filename,"qc_%s_%2.2d.%s", time_string,i,file_suffix); - printf("Creating: %s\n", filename); - RSL_sweep_to_pgm(sweep, filename, xdim, ydim, maxr); - } - if(make_gif) { - sprintf(file_suffix,"gif"); - sprintf(filename,"qc_%s_%2.2d.%s", time_string,i,file_suffix); - printf("Creating: %s\n", filename); - RSL_sweep_to_gif(sweep,filename,xdim, ydim, maxr); - } - } - } - - /* DZ_INDEX */ - if(reflectivity) { - if(verbose) printf("Loading refl colortable...\n"); - RSL_load_refl_color_table(); - for(i=0; isweep[i]; - if(sweep == NULL) { - printf("sweep[%d]==NULL\n",i); - continue; - } - if(make_pgm) { - sprintf(file_suffix,"pgm"); - sprintf(filename,"dz_%s_%2.2d.%s", - time_string,i,file_suffix); - printf("Creating: %s\n", filename); - RSL_sweep_to_pgm(sweep, filename, xdim, ydim, maxr); - } - if(make_gif) { - sprintf(file_suffix,"gif"); - sprintf(filename,"dz_%s_%2.2d.%s", time_string,i,file_suffix); -/* - sprintf(filename,"dz_%s.%s.%s", time_string,in_file,file_suffix); + /* CZ_INDEX */ + if(qc_reflectivity) { + if(verbose) printf("Loading refl colortable...\n"); + RSL_load_refl_color_table(); + for(i=0; isweep[i]; + if(sweep == NULL) { + printf("sweep[%d]==NULL\n",i); + continue; + } + if(make_pgm) { + sprintf(file_suffix,"pgm"); + sprintf(filename,"qc_%s_%2.2d.%s", time_string,i,file_suffix); + printf("Creating: %s\n", filename); + make_pathname(filename, pgmdir, pathname); + RSL_sweep_to_pgm(sweep, pathname, xdim, ydim, maxr); + } + if(make_gif) { + sprintf(file_suffix,"gif"); + sprintf(filename,"qc_%s_%2.2d.%s", time_string,i,file_suffix); + printf("Creating: %s\n", filename); + make_pathname(filename, gifdir, pathname); + RSL_sweep_to_gif(sweep,pathname,xdim, ydim, maxr); + } + } + } + + /* DZ_INDEX */ + if(reflectivity) { + if(verbose) printf("Loading refl colortable...\n"); + RSL_load_refl_color_table(); + for(i=0; isweep[i]; + if(sweep == NULL) { + printf("sweep[%d]==NULL\n",i); + continue; + } + if(make_pgm) { + sprintf(file_suffix,"pgm"); + sprintf(filename,"dz_%s_%2.2d.%s", + time_string,i,file_suffix); + printf("Creating: %s\n", filename); + make_pathname(filename, pgmdir, pathname); + RSL_sweep_to_pgm(sweep, pathname, xdim, ydim, maxr); + } + if(make_gif) { + sprintf(file_suffix,"gif"); + sprintf(filename,"dz_%s_%2.2d.%s", time_string,i,file_suffix); +/* + sprintf(filename,"dz_%s.%s.%s", time_string,in_file,file_suffix); */ - printf("Creating: %s\n", filename); - RSL_sweep_to_gif(sweep,filename,xdim, ydim, maxr); - } - } - } - - /* ZT_INDEX */ - if(total_reflectivity) { - if(verbose) printf("Loading refl colortable...\n"); - RSL_load_refl_color_table(); - for(i=0; isweep[i]; - if(sweep == NULL) { - printf("sweep[%d]==NULL\n",i); - continue; - } - if(make_pgm) { - sprintf(file_suffix,"pgm"); - sprintf(filename,"zt_%s_%2.2d.%s", - time_string,i,file_suffix); - printf("Creating: %s\n", filename); - RSL_sweep_to_pgm(sweep, filename, xdim, ydim, maxr); - } - if(make_gif) { - sprintf(file_suffix,"gif"); - sprintf(filename,"zt_%s_%2.2d.%s", - time_string,i,file_suffix); - printf("Creating: %s\n", filename); - RSL_sweep_to_gif(sweep,filename,xdim, ydim, maxr); - } - } - } - - /* DR_INDEX */ - if(differential_reflectivity) { - scale = 0.5; - ncbins = 21; - width = 10; - printf("Calling RSL_rebin, %d %d %.2f\n", width); - RSL_rebin_volume(dr_volume, width); - if(verbose) printf("Loading zdr colortable...\n"); - RSL_load_zdr_color_table(); - for(i=0; isweep[i]; - if(sweep == NULL) { - printf("sweep[%d]==NULL\n",i); - continue; - } - if(make_pgm) { - sprintf(file_suffix,"pgm"); - sprintf(filename,"dr_%s_%2.2d.%s", - time_string,i,file_suffix); - printf("Creating: %s\n", filename); - RSL_sweep_to_pgm(sweep, filename, xdim, ydim, maxr); - } - if(make_gif) { - sprintf(file_suffix,"gif"); - sprintf(filename,"dr_%s_%2.2d.%s", - time_string,i,file_suffix); - printf("Creating: %s\n", filename); - RSL_sweep_to_gif(sweep,filename,xdim, ydim, maxr); - } - } - } - - - /* VR_INDEX */ - if(velocity) { - if(verbose) printf("Loading vel colortable...\n"); - RSL_load_vel_color_table(); - for(i=0; isweep[i]; - if(sweep == NULL) { - printf("sweep[%d]==NULL\n",i); - continue; - } - if(make_pgm) { - sprintf(file_suffix,"pgm"); - sprintf(filename,"vr_%s_%2.2d.%s", time_string,i,file_suffix); - printf("Creating: %s\n", filename); - RSL_sweep_to_pgm(sweep, filename, xdim, ydim, maxr); - } - if(make_gif) { - sprintf(file_suffix,"gif"); - sprintf(filename,"vr_%s_%2.2d.%s", time_string,i,file_suffix); - printf("Creating: %s\n", filename); - RSL_sweep_to_gif(sweep,filename,xdim, ydim, maxr); - } - } - } - - /* SW_INDEX */ - if(spectral_width) { - if(verbose) printf("Loading sw colortable...\n"); - RSL_load_sw_color_table(); - for(i=0; isweep[i]; - if(sweep == NULL) { - printf("sweep[%d]==NULL\n",i); - continue; - } - if(make_pgm) { - sprintf(file_suffix,"pgm"); - sprintf(filename,"sw_%s_%2.2d.%s", - time_string,i,file_suffix); - printf("Creating: %s\n", filename); - RSL_sweep_to_pgm(sweep, filename, xdim, ydim, maxr); - } - if(make_gif) { - sprintf(file_suffix,"gif"); - sprintf(filename,"sw_%s_%2.2d.%s", - time_string,i,file_suffix); - printf("Creating: %s\n", filename); - RSL_sweep_to_gif(sweep,filename,xdim, ydim, maxr); - } - } - } - + printf("Creating: %s\n", filename); + make_pathname(filename, gifdir, pathname); + RSL_sweep_to_gif(sweep,pathname,xdim, ydim, maxr); + } + } + } + + /* ZT_INDEX */ + if(total_reflectivity) { + if(verbose) printf("Loading refl colortable...\n"); + RSL_load_refl_color_table(); + for(i=0; isweep[i]; + if(sweep == NULL) { + printf("sweep[%d]==NULL\n",i); + continue; + } + if(make_pgm) { + sprintf(file_suffix,"pgm"); + sprintf(filename,"zt_%s_%2.2d.%s", + time_string,i,file_suffix); + printf("Creating: %s\n", filename); + make_pathname(filename, pgmdir, pathname); + RSL_sweep_to_pgm(sweep, pathname, xdim, ydim, maxr); + } + if(make_gif) { + sprintf(file_suffix,"gif"); + sprintf(filename,"zt_%s_%2.2d.%s", + time_string,i,file_suffix); + printf("Creating: %s\n", filename); + make_pathname(filename, gifdir, pathname); + RSL_sweep_to_gif(sweep,pathname,xdim, ydim, maxr); + } + } + } + + /* DR_INDEX */ + if(differential_reflectivity) { + scale = 0.5; + ncbins = 21; + width = 10; + printf("Calling RSL_rebin, %d %d %.2f\n", width); + RSL_rebin_volume(dr_volume, width); + if(verbose) printf("Loading zdr colortable...\n"); + RSL_load_zdr_color_table(); + for(i=0; isweep[i]; + if(sweep == NULL) { + printf("sweep[%d]==NULL\n",i); + continue; + } + if(make_pgm) { + sprintf(file_suffix,"pgm"); + sprintf(filename,"dr_%s_%2.2d.%s", + time_string,i,file_suffix); + printf("Creating: %s\n", filename); + make_pathname(filename, pgmdir, pathname); + RSL_sweep_to_pgm(sweep, pathname, xdim, ydim, maxr); + } + if(make_gif) { + sprintf(file_suffix,"gif"); + sprintf(filename,"dr_%s_%2.2d.%s", + time_string,i,file_suffix); + printf("Creating: %s\n", filename); + make_pathname(filename, gifdir, pathname); + RSL_sweep_to_gif(sweep,pathname,xdim, ydim, maxr); + } + } + } + + + /* VR_INDEX */ + if(velocity) { + if(verbose) printf("Loading vel colortable...\n"); + RSL_load_vel_color_table(); + for(i=0; isweep[i]; + if(sweep == NULL) { + printf("sweep[%d]==NULL\n",i); + continue; + } + if(make_pgm) { + sprintf(file_suffix,"pgm"); + sprintf(filename,"vr_%s_%2.2d.%s", time_string,i,file_suffix); + printf("Creating: %s\n", filename); + make_pathname(filename, pgmdir, pathname); + RSL_sweep_to_pgm(sweep, pathname, xdim, ydim, maxr); + } + if(make_gif) { + sprintf(file_suffix,"gif"); + sprintf(filename,"vr_%s_%2.2d.%s", time_string,i,file_suffix); + printf("Creating: %s\n", filename); + make_pathname(filename, gifdir, pathname); + RSL_sweep_to_gif(sweep,pathname,xdim, ydim, maxr); + } + } + } + + /* SW_INDEX */ + if(spectral_width) { + if(verbose) printf("Loading sw colortable...\n"); + RSL_load_sw_color_table(); + for(i=0; isweep[i]; + if(sweep == NULL) { + printf("sweep[%d]==NULL\n",i); + continue; + } + if(make_pgm) { + sprintf(file_suffix,"pgm"); + sprintf(filename,"sw_%s_%2.2d.%s", + time_string,i,file_suffix); + printf("Creating: %s\n", filename); + make_pathname(filename, pgmdir, pathname); + RSL_sweep_to_pgm(sweep, pathname, xdim, ydim, maxr); + } + if(make_gif) { + sprintf(file_suffix,"gif"); + sprintf(filename,"sw_%s_%2.2d.%s", + time_string,i,file_suffix); + printf("Creating: %s\n", filename); + make_pathname(filename, gifdir, pathname); + RSL_sweep_to_gif(sweep,pathname,xdim, ydim, maxr); + } + } + } + /* Write uf file if requested */ - if(make_uf) { - sprintf(file_suffix,"uf.gz"); - sprintf(filename,"%s_%s.%s",site_id, time_string,file_suffix); - printf("Creating UF file: %s\n", filename); - RSL_radar_to_uf_gzip(radar, filename); - } - - printf("-->> FIELDS: [ "); - if(radar->v[0] != NULL) printf("DZ "); - if(radar->v[1] != NULL) printf("VR "); - if(radar->v[2] != NULL) printf("SW "); - if(radar->v[3] != NULL) printf("CZ "); - if(radar->v[4] != NULL) printf("ZT "); - if(radar->v[5] != NULL) printf("DR "); - if(radar->v[6] != NULL) printf("LR "); - if(radar->v[7] != NULL) printf("ZD "); - if(radar->v[8] != NULL) printf("DM "); - if(radar->v[9] != NULL) printf("RH "); - if(radar->v[10] != NULL) printf("PH "); - if(radar->v[11] != NULL) printf("XZ "); - if(radar->v[12] != NULL) printf("CR "); - if(radar->v[13] != NULL) printf("MZ "); - if(radar->v[14] != NULL) printf("MR "); - if(radar->v[15] != NULL) printf("ZE "); - if(radar->v[16] != NULL) printf("VE "); - if(radar->v[17] != NULL) printf("KD "); - if(radar->v[18] != NULL) printf("TI "); - if(radar->v[19] != NULL) printf("DX "); - if(radar->v[20] != NULL) printf("CH "); - if(radar->v[21] != NULL) printf("AH "); - if(radar->v[22] != NULL) printf("CV "); - if(radar->v[23] != NULL) printf("AV "); - if(radar->v[24] != NULL) printf("SQ "); - printf("] \n\n"); + if(make_uf) { + sprintf(file_suffix,"uf.gz"); + sprintf(filename,"%s_%s.%s",site_id, time_string,file_suffix); + printf("Creating UF file: %s\n", filename); + make_pathname(filename, ufdir, pathname); + RSL_radar_to_uf_gzip(radar, pathname); + } + + printf("-->> FIELDS: [ "); + /* Modified to use RSL_ftype from rsl.h (#define USE_RSL_VARS) and to + * loop through volumes. 10/16/2009, BLK + */ + for (i=0; i < MAX_RADAR_VOLUMES; i++) + if (radar->v[i] != NULL) printf("%s ", RSL_ftype[i]); + /* + if(radar->v[0] != NULL) printf("DZ "); + if(radar->v[1] != NULL) printf("VR "); + if(radar->v[2] != NULL) printf("SW "); + if(radar->v[3] != NULL) printf("CZ "); + if(radar->v[4] != NULL) printf("ZT "); + if(radar->v[5] != NULL) printf("DR "); + if(radar->v[6] != NULL) printf("LR "); + if(radar->v[7] != NULL) printf("ZD "); + if(radar->v[8] != NULL) printf("DM "); + if(radar->v[9] != NULL) printf("RH "); + if(radar->v[10] != NULL) printf("PH "); + if(radar->v[11] != NULL) printf("XZ "); + if(radar->v[12] != NULL) printf("CR "); + if(radar->v[13] != NULL) printf("MZ "); + if(radar->v[14] != NULL) printf("MR "); + if(radar->v[15] != NULL) printf("ZE "); + if(radar->v[16] != NULL) printf("VE "); + if(radar->v[17] != NULL) printf("KD "); + if(radar->v[18] != NULL) printf("TI "); + if(radar->v[19] != NULL) printf("DX "); + if(radar->v[20] != NULL) printf("CH "); + if(radar->v[21] != NULL) printf("AH "); + if(radar->v[22] != NULL) printf("CV "); + if(radar->v[23] != NULL) printf("AV "); + if(radar->v[24] != NULL) printf("SQ "); + */ + printf("] \n\n"); /* Wrap it up! */ if(verbose) - printf("Finished!\n"); + printf("Finished!\n"); exit (0); } /* End of main */ diff --git a/examples/qlook_usage.c b/examples/qlook_usage.c index 80ca2e5..4eef24b 100644 --- a/examples/qlook_usage.c +++ b/examples/qlook_usage.c @@ -7,32 +7,33 @@ #include void process_args(int argc, char **argv, char *in_file, int *verbose, - char *site_id, char *tape_id, - int *qc_reflectivity, int *total_reflectivity, - int *differential_reflectivity, - int *velocity, int *spectral_width, - int *make_gif, int *make_pgm, int *make_bscan, int *make_uf, - int *num_sweeps, float *dbz_offset, - int *xdim, int *ydim, float *range, - float *gate_size_adjustment, int *print_azim) + char *site_id, char *tape_id, + int *qc_reflectivity, int *total_reflectivity, + int *differential_reflectivity, + int *velocity, int *spectral_width, + int *make_gif, int *make_pgm, int *make_bscan, int *make_uf, + int *num_sweeps, float *dbz_offset, + int *xdim, int *ydim, float *range, + float *gate_size_adjustment, int *print_azim, + char *gifdir, char *pgmdir, char *ufdir) { - extern char *optarg; + extern char *optarg; extern int optind, optopt; - char c; + char c; - while ((c = getopt(argc, argv, "vgpus:t:n:x:y:r:o:a:ADCQTWV")) != -1) { + while ((c = getopt(argc, argv, "vgpus:t:n:x:y:r:o:a:ADCQTWVG:P:U:")) != -1) { - switch(c) { + switch(c) { /* RSL Verbose flag */ - case 'v': *verbose = TRUE; break; + case 'v': *verbose = TRUE; break; /* s: First file or call sign */ - case 's': strcpy(site_id, optarg); break; - case 't': strcpy(tape_id, optarg); break; + case 's': strcpy(site_id, optarg); break; + case 't': strcpy(tape_id, optarg); break; /* x: x dimension @@ -40,83 +41,95 @@ void process_args(int argc, char **argv, char *in_file, int *verbose, r: max range z: zoom factor (km/pixel) */ - case 'x': *xdim = atoi(optarg); break; - case 'y': *ydim = atoi(optarg); break; - case 'r': *range = atof(optarg); break; - case 'a': *gate_size_adjustment = atof(optarg); break; - + case 'x': *xdim = atoi(optarg); break; + case 'y': *ydim = atoi(optarg); break; + case 'r': *range = atof(optarg); break; + case 'a': *gate_size_adjustment = atof(optarg); break; + /* dBZ Offset */ - case 'o': *dbz_offset = atof(optarg); break; + case 'o': *dbz_offset = atof(optarg); break; /* T: Total reflectivity Q: Do qc'd reflectivity V: Do radial velocity W: Do spectral width */ - case 'Q': *qc_reflectivity = TRUE; break; - case 'T': *total_reflectivity = TRUE; break; - case 'V': *velocity = TRUE; break; - case 'W': *spectral_width = TRUE; break; - case 'A': *print_azim = TRUE; break; - case 'D': *differential_reflectivity = TRUE; break; + case 'Q': *qc_reflectivity = TRUE; break; + case 'T': *total_reflectivity = TRUE; break; + case 'V': *velocity = TRUE; break; + case 'W': *spectral_width = TRUE; break; + case 'A': *print_azim = TRUE; break; + case 'D': *differential_reflectivity = TRUE; break; /* g: Make gif images p: Make pgm images u: Make uf files */ - case 'g': *make_gif = TRUE; break; - case 'p': *make_pgm = TRUE; break; - case 'u': *make_uf = TRUE; break; + case 'g': *make_gif = TRUE; break; + case 'p': *make_pgm = TRUE; break; + case 'u': *make_uf = TRUE; break; + +/* + G: gif directory + P: pgm directory + U: uf directory +*/ + case 'G': strcpy(gifdir, optarg); break; + case 'P': strcpy(pgmdir, optarg); break; + case 'U': strcpy(ufdir, optarg); break; /* num_sweeps: Number of sweeps to make images of */ - case 'n': *num_sweeps = atoi(optarg); break; + case 'n': *num_sweeps = atoi(optarg); break; /* Deal with bad input */ - case '?': fprintf(stderr, "ERROR: option -%c is undefined\n", optopt); - goto Usage; - case ':': fprintf(stderr, "ERROR: option -%c requires an argument\n",optopt); - goto Usage; - default: break; - } - } + case '?': fprintf(stderr, "ERROR: option -%c is undefined\n", optopt); + goto Usage; + case ':': fprintf(stderr, "ERROR: option -%c requires an argument\n",optopt); + goto Usage; + default: break; + } + } /* Must have at the least a file listed on the command lines, everything can be defaulted. */ - if (argc - optind != 1) { + if (argc - optind != 1) { Usage: - fprintf(stderr,"ERROR:::\n"); - fprintf(stderr,"%s [options] input_file:",argv[0]); - fprintf(stderr,"\n[options]: "); - fprintf(stderr,"\n\t[-v verbose_flag?] "); - fprintf(stderr,"\n\t[-s First file or call sign?] "); - fprintf(stderr,"\n\t[-t Tape ID] "); - fprintf(stderr,"\n\t[-u Make UF file]"); - fprintf(stderr,"\n\t[-g Make GIF images?]"); - fprintf(stderr,"\n\t[-p Make PGM images?]"); - fprintf(stderr,"\n\t[-x X dimension]"); - fprintf(stderr,"\n\t[-y Y dimension]"); - fprintf(stderr,"\n\t[-r max range]"); - fprintf(stderr,"\n\t[-n Number of sweeps to make images]"); - fprintf(stderr,"\n\t[-Q Do qc reflectivity]"); - fprintf(stderr,"\n\t[-T Do total reflectivity]"); - fprintf(stderr,"\n\t[-V Do velocity]"); - fprintf(stderr,"\n\t[-W Do spectral_width]"); - fprintf(stderr,"\n\t[-D Do differential reflectivity"); - fprintf(stderr,"\n\t[-o Apply dBZ offset"); - fprintf(stderr,":::\n"); - exit(-1); - } - - strcpy(in_file, argv[optind]); + fprintf(stderr,"ERROR:::\n"); + fprintf(stderr,"%s [options] input_file:",argv[0]); + fprintf(stderr,"\n[options]: "); + fprintf(stderr,"\n\t[-v verbose_flag?] "); + fprintf(stderr,"\n\t[-s First file or call sign?] "); + fprintf(stderr,"\n\t[-t Tape ID] "); + fprintf(stderr,"\n\t[-u Make UF file]"); + fprintf(stderr,"\n\t[-g Make GIF images?]"); + fprintf(stderr,"\n\t[-p Make PGM images?]"); + fprintf(stderr,"\n\t[-U Directory for UF output files]"); + fprintf(stderr,"\n\t[-G Directory for GIF output files]"); + fprintf(stderr,"\n\t[-P Directory for PGM output files]"); + fprintf(stderr,"\n\t[-x X dimension]"); + fprintf(stderr,"\n\t[-y Y dimension]"); + fprintf(stderr,"\n\t[-r max range]"); + fprintf(stderr,"\n\t[-n Number of sweeps to make images]"); + fprintf(stderr,"\n\t[-Q Do qc reflectivity]"); + fprintf(stderr,"\n\t[-T Do total reflectivity]"); + fprintf(stderr,"\n\t[-V Do velocity]"); + fprintf(stderr,"\n\t[-W Do spectral_width]"); + fprintf(stderr,"\n\t[-D Do differential reflectivity"); + fprintf(stderr,"\n\t[-o Apply dBZ offset"); + fprintf(stderr,":::\n"); + exit(-1); + } + + strcpy(in_file, argv[optind]); } diff --git a/examples/wsr_hist_uf_test.c b/examples/wsr_hist_uf_test.c index 615b572..c7b129a 100644 --- a/examples/wsr_hist_uf_test.c +++ b/examples/wsr_hist_uf_test.c @@ -20,19 +20,28 @@ #endif #include #include +#include #include "rsl.h" usage() { - fprintf(stderr,"Usage: wsr_hist_uf_test infile\n"); + fprintf(stderr,"Usage: wsr_hist_uf_test infile [-s site_id]\n"); exit(-1); } -process_args(int argc, char **argv, char **in_file) +process_args(int argc, char **argv, char **in_file, char **site) { - if (argc == 2) *in_file = strdup(argv[1]); + int c; + + while ((c = getopt(argc, argv, "s:")) != -1) + switch (c) { + case 's': *site = strdup(optarg); break; + case '?': usage(argv); break; + default: break; + } + if (argc - optind == 1) *in_file = strdup(argv[optind]); else usage(); } @@ -40,14 +49,15 @@ process_args(int argc, char **argv, char **in_file) main(int argc, char **argv) { char *infile; + char *site = NULL; Radar *radar; Histogram *histogram = NULL; - process_args(argc, argv, &infile); + process_args(argc, argv, &infile, &site); RSL_radar_verbose_on(); - if ((radar = RSL_anyformat_to_radar(infile, "KMLB")) == NULL) { + if ((radar = RSL_anyformat_to_radar(infile, site)) == NULL) { /* RSL_wsr88d_to_radar writes an error message to stdout. */ exit(-1); } diff --git a/mkinstalldirs b/mkinstalldirs new file mode 100644 index 0000000..5b16fbe --- /dev/null +++ b/mkinstalldirs @@ -0,0 +1,40 @@ +#! /bin/sh +# mkinstalldirs --- make directory hierarchy +# Author: Noah Friedman +# Created: 1993-05-16 +# Public domain + +# $Id: mkinstalldirs,v 1.1 1999/12/06 05:21:06 merritt Exp $ + +errstatus=0 + +for file +do + set fnord `echo ":$file" | sed -ne 's/^:\//#/;s/^://;s/\// /g;s/^#/\//;p'` + shift + + pathcomp= + for d + do + pathcomp="$pathcomp$d" + case "$pathcomp" in + -* ) pathcomp=./$pathcomp ;; + esac + + if test ! -d "$pathcomp"; then + echo "mkdir $pathcomp" + + mkdir "$pathcomp" || lasterr=$? + + if test ! -d "$pathcomp"; then + errstatus=$lasterr + fi + fi + + pathcomp="$pathcomp/" + done +done + +exit $errstatus + +# mkinstalldirs ends here diff --git a/nsig_to_radar.c b/nsig_to_radar.c index 8454ab1..88788e4 100644 --- a/nsig_to_radar.c +++ b/nsig_to_radar.c @@ -71,7 +71,7 @@ extern int rsl_qfield[]; /* See RSL_select_fields */ static float (*f)(Range x); static Range (*invf)(float x); -FILE *file; +extern FILE *file; void get_extended_header_info(NSIG_Sweep **nsig_sweep, int xh_size, int iray, int nparams, @@ -165,7 +165,9 @@ RSL_nsig_to_radar float second; float pw; float bin_space; - float prf, wave, beam_width; + float prf, prf2, wave, beam_width; + int prf_mode; + float prf_modes[] = {1.0, 2.0/3.0, 3.0/4.0, 4.0/5.0}; float vert_half_bw, horz_half_bw; float rng_last_bin; float rng_first_bin, freq; @@ -184,7 +186,8 @@ RSL_nsig_to_radar NSIG_Sweep **nsig_sweep; NSIG_Ray *ray_p; int itype, ifield; - unsigned short nsig_u2byte; /* New for 2-byte data types, Aug 2009 */ + unsigned short nsig_2byte; /* New for 2-byte data types, Aug 2009 */ + twob nsig_twob; Sweep *sweep; int msec; float azm, elev, pitch, roll, heading, azm_rate, elev_rate, @@ -374,6 +377,8 @@ RSL_nsig_to_radar num_rays = 0; pw = (NSIG_I4(prod_file->rec1.prod_end.pulse_wd))/100.0; /* pulse width */ prf = NSIG_I4(prod_file->rec1.prod_end.prf); /* pulse repetition frequency */ + prf_mode = NSIG_I2(prod_file->rec2.task_config.dsp_info.prf_mode); + prf2 = prf * prf_modes[prf_mode]; 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. @@ -387,7 +392,18 @@ RSL_nsig_to_radar num_samples = NSIG_I2(prod_file->rec1.prod_end.num_samp); sweep_rate = 3.0; /** Approximate value -- info not stored **/ azim_rate = sweep_rate*360.0/60.0; - max_vel = wave*prf/(100.0*4.0); + if (prf_mode != 0) + { + float max_vel1 = wave*prf/(100.0*4.0); + float max_vel2 = wave*prf2/(100.0*4.0); + + max_vel = (max_vel1 * max_vel2)/(max_vel1-max_vel2); + } + else + { + max_vel = wave*prf/(100.0*4.0); + } + freq = (299793000.0/wave)*1.0e-4; /** freq in MHZ **/ sqi = NSIG_I2(prod_file->rec2.task_config.calib_info.sqi)/256.0; @@ -421,7 +437,7 @@ RSL_nsig_to_radar } if (radar_verbose_flag) - fprintf(stderr, "vel: %f prf: %f\n", max_vel, prf); + fprintf(stderr, "vel: %f prf: %f prf2: %f\n", max_vel, prf, prf2); /** Extracting Latitude and Longitude from nsig file **/ lat = nsig_from_fourb_ang(prod_file->rec2.ingest_head.lat_rad); @@ -708,9 +724,10 @@ RSL_nsig_to_radar 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.prf2 = (int) prf2; + ray->h.fix_angle = (float)sweep->h.elev; ray->h.azim_rate = azim_rate; - ray->h.pulse_count = (float)num_samples; + ray->h.pulse_count = num_samples; ray->h.pulse_width = pw; ray->h.beam_width = beam_width; ray->h.frequency = freq / 1000.0; /* GHz */ @@ -842,46 +859,46 @@ RSL_nsig_to_radar 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) + memmove(nsig_twob, &ray_p->range[2*k], 2); + nsig_2byte = NSIG_I2(nsig_twob); + if (nsig_2byte == 0 || nsig_2byte == 65535) ray_data = NSIG_NO_ECHO2; - else ray_data = (float)(nsig_u2byte-32768)/100.; + else ray_data = (float)(nsig_2byte-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) + memmove(nsig_twob, &ray_p->range[2*k], 2); + nsig_2byte = NSIG_I2(nsig_twob); + if (nsig_2byte == 0 || nsig_2byte == 65535) ray_data = NSIG_NO_ECHO2; - else ray_data = (float)nsig_u2byte/100.; + else ray_data = (float)nsig_2byte/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) + memmove(nsig_twob, &ray_p->range[2*k], 2); + nsig_2byte = NSIG_I2(nsig_twob); + if (nsig_2byte == 0 || nsig_2byte == 65535) ray_data = NSIG_NO_ECHO; else - ray_data = 360.*(nsig_u2byte-1)/65534.; + ray_data = 360.*(nsig_2byte-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) + memmove(nsig_twob, &ray_p->range[2*k], 2); + nsig_2byte = NSIG_I2(nsig_twob); + if (nsig_2byte == 0 || nsig_2byte == 65535) ray_data = NSIG_NO_ECHO2; - else ray_data = (float)(nsig_u2byte-1)/65533.; + else ray_data = (float)(nsig_2byte-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) + memmove(nsig_twob, &ray_p->range[2*k], 2); + nsig_2byte = NSIG_I2(nsig_twob); + if (nsig_2byte == 0 || nsig_2byte == 65535) ray_data = NSIG_NO_ECHO2; else - ray_data = nsig_u2byte; + ray_data = nsig_2byte; } if (ray_data == NSIG_NO_ECHO || ray_data == NSIG_NO_ECHO2) diff --git a/prune.c b/prune.c index 1a23913..6ee0459 100644 --- a/prune.c +++ b/prune.c @@ -36,21 +36,6 @@ #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; @@ -65,20 +50,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; @@ -91,19 +76,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; @@ -115,11 +100,8 @@ Radar *RSL_prune_radar(Radar *radar) int i; /* Volume indexes are fixed so we just prune the substructures. */ if (radar == NULL) return NULL; - 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); + for (i=0; ih.nvolumes; i++) + radar->v[i] = RSL_prune_volume(radar->v[i]); return radar; } diff --git a/radar.c b/radar.c index ff138c8..4078f17 100644 --- a/radar.c +++ b/radar.c @@ -123,6 +123,7 @@ Radar *RSL_new_radar(int nvolumes) r = (Radar *) calloc(1, sizeof(Radar)); r->v = (Volume **) calloc(nvolumes, sizeof(Volume *)); r->h.nvolumes = nvolumes; + r->h.scan_mode = PPI; /* default PPI is enum constant defined in rsl.h */ return r; } diff --git a/radar_to_uf.c b/radar_to_uf.c index 6708396..df237f5 100644 --- a/radar_to_uf.c +++ b/radar_to_uf.c @@ -25,13 +25,14 @@ #include #include +#define USE_RSL_VARS #include "rsl.h" extern int radar_verbose_flag; /* Missing data flag : -32768 when a signed short. */ #define UF_NO_DATA 0X8000 -/* Any convensions may be observed. */ +/* Field names. Any convensions may be observed. */ /* Typically: * DZ = Reflectivity (dBZ). * VR = Radial Velocity. @@ -53,14 +54,9 @@ extern int radar_verbose_flag; * KD = KDP wavelength*deg/km * TI = TIME (units unknown). * These fields may appear in any order in the UF file. + * There are more fields than appear here. See rsl.h. */ -char *UF_field_name[] = {"DZ", "VR", "SW", "CZ", "ZT", "DR", "LR", - "ZD", "DM", "RH", "PH", "XZ", "CD", "MZ", - "MD", "ZE", "VE", "KD", "TI", "DX", "CH", - "AH", "CV", "AV", "SQ" -}; - typedef short UF_buffer[16384]; /* Bigger than documented 4096. */ @@ -110,6 +106,7 @@ void RSL_radar_to_uf_fp(Radar *r, FILE *fp) int rec_len, save_rec_len; int nfield; float vr_az; + int max_field_names; struct tm *tm; time_t the_time; @@ -118,6 +115,7 @@ void RSL_radar_to_uf_fp(Radar *r, FILE *fp) int degree, minute; float second; + int uf_sweep_mode = 1; /* default PPI */ /* Here are the arrays for each field type. Each dimension is the number * of fields in the radar structure. I do this because the radar organization @@ -134,11 +132,10 @@ void RSL_radar_to_uf_fp(Radar *r, FILE *fp) float x; if (r == NULL) { - fprintf(stderr, "radar_to_uf_fp: radar pointer NULL\n"); - return; + fprintf(stderr, "radar_to_uf_fp: radar pointer NULL\n"); + return; } - /* Do all the headers first time around. Then, prune OP and LU. */ q_op = q_lu = q_dh = q_fh = 1; @@ -147,6 +144,10 @@ void RSL_radar_to_uf_fp(Radar *r, FILE *fp) sweep_num = ray_num = rec_num = 0; true_nvolumes = nvolumes = maxsweeps = nrays = 0; + /* PPI and RHI are enum constants defined in rsl.h */ + if (r->h.scan_mode == PPI) uf_sweep_mode = 1; + else if (r->h.scan_mode == RHI) uf_sweep_mode = 3; + /* * The organization of the Radar structure is by volumes, then sweeps, then * rays, then gates. This is different from the UF file organization. @@ -171,325 +172,339 @@ void RSL_radar_to_uf_fp(Radar *r, FILE *fp) * the main controlling loop variable. */ for (i=0; iv[i]; - if(volume[i]) { - nsweeps[i] = volume[i]->h.nsweeps; - if (nsweeps[i] > maxsweeps) maxsweeps = nsweeps[i]; - true_nvolumes++; - } + volume[i] = r->v[i]; + if(volume[i]) { + nsweeps[i] = volume[i]->h.nsweeps; + if (nsweeps[i] > maxsweeps) maxsweeps = nsweeps[i]; + true_nvolumes++; + } } if (radar_verbose_flag) { - fprintf(stderr,"True number of volumes for UF is %d\n", true_nvolumes); - fprintf(stderr,"Maximum # of volumes for UF is %d\n", nvolumes); + fprintf(stderr,"True number of volumes for UF is %d\n", true_nvolumes); + fprintf(stderr,"Maximum # of volumes for UF is %d\n", nvolumes); } + + max_field_names = sizeof(RSL_ftype) / 4; + /*-------- * LOOP for all sweeps (typically 11 or 16 for wsr88d data. * */ for (i=0; isweep[i]; - - /* Check if we really can access this sweep. Paul discovered that - * if the actual number of sweeps is less than the maximum that we - * could be chasing a bad pointer (a NON-NULL garbage pointer). - */ - if (i >= nsweeps[k]) sweep[k] = NULL; - - if (sweep[k]) if (sweep[k]->h.nrays > nrays) nrays = sweep[k]->h.nrays; - } - - sweep_num++; /* I guess it will be ok to count NULL sweeps. */ - ray_num = 0; + nrays = 0; + for (k=0; ksweep[i]; + + /* Check if we really can access this sweep. Paul discovered that + * if the actual number of sweeps is less than the maximum that we + * could be chasing a bad pointer (a NON-NULL garbage pointer). + */ + if (i >= nsweeps[k]) sweep[k] = NULL; + + if (sweep[k]) if (sweep[k]->h.nrays > nrays) nrays = sweep[k]->h.nrays; + } + + sweep_num++; /* I guess it will be ok to count NULL sweeps. */ + ray_num = 0; if (radar_verbose_flag) - fprintf(stderr,"Processing sweep %d for %d rays.", i, nrays); + fprintf(stderr,"Processing sweep %d for %d rays.", i, nrays); if (radar_verbose_flag) - if (little_endian()) fprintf(stderr," ... On Little endian.\n"); - else fprintf(stderr,"\n"); + if (little_endian()) fprintf(stderr," ... On Little endian.\n"); + else fprintf(stderr,"\n"); /* Now LOOP for all rays within this particular sweep (i). * Get all the field types together for the ray, see ray[k], and * fill the UF data buffer appropriately. */ - for (j=0; jh.nrays) - if (sweep[k]->ray) - if ((ray = sweep[k]->ray[j])) break; - } - - /* If there is no such ray, then continue on to the next ray. */ - if (ray) { + for (j=0; jh.nrays) + if (sweep[k]->ray) + if ((ray = sweep[k]->ray[j])) break; + } + + /* If there is no such ray, then continue on to the next ray. */ + if (ray) { /* - fprintf(stderr,"Ray: %.4d, Time: %2.2d:%2.2d:%f %.2d/%.2d/%.4d\n", ray_num, ray->h.hour, ray->h.minute, ray->h.sec, ray->h.month, ray->h.day, ray->h.year); + fprintf(stderr,"Ray: %.4d, Time: %2.2d:%2.2d:%f %.2d/%.2d/%.4d\n", ray_num, ray->h.hour, ray->h.minute, ray->h.sec, ray->h.month, ray->h.day, ray->h.year); */ - /* - * ---- Begining of MANDITORY HEADER BLOCK. - */ - uf_ma = uf; - memcpy(&uf_ma[0], "UF", 2); - if (little_endian()) memcpy(&uf_ma[0], "FU", 2); - uf_ma[1] = 0; /* Not known yet. */ - uf_ma[2] = 0; /* Not known yet. Really, I do. */ - uf_ma[3] = 0; /* Not known yet. */ - uf_ma[4] = 0; /* Not known yet. */ - - uf_ma[6] = 1; - uf_ma[7] = ray_num; - uf_ma[8 ] = 1; - uf_ma[9 ] = sweep_num; - memcpy(&uf_ma[10], r->h.radar_name, 8); - if (little_endian()) swap2(&uf_ma[10], 8/2); - memcpy(&uf_ma[14], r->h.name, 8); - if (little_endian()) swap2(&uf_ma[14], 8/2); - /* Convert decimal lat/lon to d:m:s */ - - if (ray->h.lat != 0.0) { - degree = (int)ray->h.lat; - minute = (int)((ray->h.lat - degree) * 60); - second = (ray->h.lat - degree - minute/60.0) * 3600.0; - } else { - degree = r->h.latd; - minute = r->h.latm; - second = r->h.lats; - } - uf_ma[18] = degree; - uf_ma[19] = minute; - if (second > 0.0) uf_ma[20] = second*64 + 0.5; - else uf_ma[20] = second*64 - 0.5; - - if (ray->h.lon != 0.0) { - degree = (int)ray->h.lon; - minute = (int)((ray->h.lon - degree) * 60); - second = (ray->h.lon - degree - minute/60.0) * 3600.0; - } else { - degree = r->h.lond; - minute = r->h.lonm; - second = r->h.lons; - } - uf_ma[21] = degree; - uf_ma[22] = minute; - if (second > 0.0) uf_ma[23] = second*64 + 0.5; - else uf_ma[23] = second*64 - 0.5; - if (ray->h.alt != 0) - uf_ma[24] = ray->h.alt; - else - uf_ma[24] = r->h.height; - - uf_ma[25] = ray->h.year % 100; /* By definition: not year 2000 compliant. */ - uf_ma[26] = ray->h.month; - uf_ma[27] = ray->h.day; - uf_ma[28] = ray->h.hour; - uf_ma[29] = ray->h.minute; - uf_ma[30] = ray->h.sec; - memcpy(&uf_ma[31], "UT", 2); - if (little_endian()) memcpy(&uf_ma[31], "TU", 2); - if (ray->h.azimuth > 0) uf_ma[32] = ray->h.azimuth*64 + 0.5; - else uf_ma[32] = ray->h.azimuth*64 - 0.5; - uf_ma[33] = ray->h.elev*64 + 0.5; - uf_ma[34] = 1; /* Sweep mode: PPI = 1 */ - if (ray->h.fix_angle != 0.) - uf_ma[35] = ray->h.fix_angle*64.0 + 0.5; - else uf_ma[35] = sweep[k]->h.elev*64.0 + 0.5; - uf_ma[36] = ray->h.sweep_rate*(360.0/60.0)*64.0 + 0.5; - - the_time = time(NULL); - tm = gmtime(&the_time); - - uf_ma[37] = tm->tm_year % 100; /* Same format as data year */ - uf_ma[38] = tm->tm_mon+1; - uf_ma[39] = tm->tm_mday; - memcpy(&uf_ma[40], "RSL" RSL_VERSION_STR, 8); - if (little_endian()) swap2(&uf_ma[40], 8/2); - uf_ma[44] = (signed short)UF_NO_DATA; - len_ma = 45; - uf_ma[2] = len_ma+1; - /* - * ---- End of MANDITORY HEADER BLOCK. - */ - - /* ---- Begining of OPTIONAL HEADER BLOCK. */ - len_op = 0; - if (q_op) { - q_op = 0; /* Only once. */ - uf_op = uf+len_ma; - memcpy(&uf_op[0], "TRMMGVUF", 8); - if (little_endian()) swap2(&uf_op[0], 8/2); - uf_op[4] = (signed short)UF_NO_DATA; - uf_op[5] = (signed short)UF_NO_DATA; - uf_op[6] = ray->h.hour; - uf_op[7] = ray->h.minute; - uf_op[8] = ray->h.sec; - memcpy(&uf_op[9], "RADAR_UF", 8); - if (little_endian()) swap2(&uf_op[9], 8/2); - uf_op[13] = 2; - len_op = 14; - } - /* ---- End of OPTIONAL HEADER BLOCK. */ - - /* ---- Begining of LOCAL USE HEADER BLOCK. */ - /* If we have DZ and VR, check to see if their azimuths are - * different. If they are, we store VR azimuth in Local Use - * Header. These differences occur with WSR-88D radars, which - * run separate sweeps for DZ and VR at low elevations. - */ - q_lu = 0; - if (sweep[DZ_INDEX] && sweep[VR_INDEX]) { - if (sweep[DZ_INDEX]->ray[j] && sweep[VR_INDEX]->ray[j]) { - vr_az = sweep[VR_INDEX]->ray[j]->h.azimuth; - if (sweep[DZ_INDEX]->ray[j]->h.azimuth != vr_az) - q_lu = 1; /* Set to use Local Use Header block. */ - } - } - len_lu = 0; - if (q_lu) { - /* Store azimuth for WSR-88D VR ray in Local Use Header. */ - uf_lu = uf+len_ma+len_op; - memcpy(&uf_lu[0], "AZ", 2); - if (little_endian()) memcpy(&uf_lu[0], "ZA", 2); - if (vr_az > 0) uf_lu[1] = vr_az*64 + 0.5; - else uf_lu[1] = vr_az*64 - 0.5; - len_lu = 2; - } - /* ---- End of LOCAL USE HEADER BLOCK. */ - - - /* Here is where we loop on each field type. We need to keep - * track of how many FIELD HEADER and FIELD DATA sections, one - * for each field type, we fill. The variable that tracks this - * index into 'uf' is 'current_fh_index'. It is bumped by - * the length of the FIELD HEADER and FIELD DATA for each field - * type encountered. Field types expected are: Reflectivity, - * Velocity, and Spectrum width; this is a typicial list but it - * is not restricted to it. - */ - - for (k=0; kray) - ray = sweep[k]->ray[j]; - else - ray = NULL; - else ray = NULL; - - if (k >= sizeof(UF_field_name)) { - ray = NULL; - fprintf(stderr,"RSL_uf_to_radar: Unknown field type encountered.\n"); - fprintf(stderr,"The field type index in Radar exceeds the number of known UF field types.\n"); - } - - if (ray) { - /* ---- Begining of DATA HEADER. */ - nfield++; - if (q_dh) { - len_dh = 2*true_nvolumes + 3; - uf_dh = uf+len_ma+len_op+len_lu; - uf_dh[0] = nfield; - uf_dh[1] = 1; - uf_dh[2] = nfield; - /* 'nfield' indexes the field number. - * 'k' indexes the particular field from the volume. - */ - memcpy(&uf_dh[3+2*(nfield-1)], UF_field_name[k], 2); - if (little_endian()) swap2(&uf_dh[3+2*(nfield-1)], 2/2); - if (current_fh_index == 0) current_fh_index = len_ma+len_op+len_lu+len_dh; - uf_dh[4+2*(nfield-1)] = current_fh_index + 1; - } - /* ---- End of DATA HEADER. */ - - /* ---- Begining of FIELD HEADER. */ - if (q_fh) { - uf_fh = uf+current_fh_index; - uf_fh[1] = scale_factor = 100; - uf_fh[2] = ray->h.range_bin1/1000.0; - uf_fh[3] = ray->h.range_bin1 - (1000*uf_fh[2]); - uf_fh[4] = ray->h.gate_size; - uf_fh[5] = ray->h.nbins; - uf_fh[6] = ray->h.pulse_width*(RSL_SPEED_OF_LIGHT/1.0e6); - uf_fh[7] = sweep[k]->h.beam_width*64.0 + 0.5; - uf_fh[8] = sweep[k]->h.beam_width*64.0 + 0.5; - uf_fh[9] = ray->h.frequency*64.0 + 0.5; /* Bandwidth (mHz). */ - uf_fh[10] = 0; /* Horizontal polarization. */ - uf_fh[11] = ray->h.wavelength*64.0*100.0; /* m to cm. */ - uf_fh[12] = ray->h.pulse_count; - memcpy(&uf_fh[13], " ", 2); - uf_fh[14] = (signed short)UF_NO_DATA; - uf_fh[15] = (signed short)UF_NO_DATA; - if (k == DZ_INDEX || k == ZT_INDEX) { - uf_fh[16] = volume[k]->h.calibr_const*100.0 + 0.5; - } - else { - memcpy(&uf_fh[16], " ", 2); - } - if (ray->h.prf != 0) - uf_fh[17] = 1.0/ray->h.prf*1000000.0; /* Pulse repetition time(msec) = 1/prf */ - else - uf_fh[17] = (signed short)UF_NO_DATA; /* Pulse repetition time = 1/prf */ - uf_fh[18] = 16; - if (VR_INDEX == k || VE_INDEX == k) { - uf_fh[19] = scale_factor*ray->h.nyq_vel; - uf_fh[20] = 1; - len_fh = 21; - } else { - len_fh = 19; - } - - uf_fh[0] = current_fh_index + len_fh + 1; - /* ---- End of FIELD HEADER. */ - - /* ---- Begining of FIELD DATA. */ - uf_data = uf+len_fh+current_fh_index; - len_data = ray->h.nbins; - for (m=0; mh.f(ray->range[m]); - if (x == BADVAL || x == RFVAL || x == APFLAG || x == NOECHO) - uf_data[m] = (signed short)UF_NO_DATA; - else - uf_data[m] = scale_factor * x; - } - - current_fh_index += (len_fh+len_data); - } - } - /* ---- End of FIELD DATA. */ - } - /* Fill in some infomation we didn't know. Like, buffer length, - * record number, etc. - */ - rec_num++; - uf_ma[1] = current_fh_index; - uf_ma[3] = len_ma + len_op + 1; - uf_ma[4] = len_ma + len_op + len_lu + 1; - uf_ma[5] = rec_num; - - /* WRITE the UF buffer. */ - rec_len =(int)uf_ma[1]*2; - save_rec_len = rec_len; /* We destroy 'rec_len' when making it - big endian on a little endian machine. */ - if (little_endian()) swap_4_bytes(&rec_len); - (void)fwrite(&rec_len, sizeof(int), 1, fp); - if (little_endian()) swap_uf_buffer(uf); - (void)fwrite(uf, sizeof(char), save_rec_len, fp); - (void)fwrite(&rec_len, sizeof(int), 1, fp); - } /* if (ray) */ - } + /* + * ---- Begining of MANDITORY HEADER BLOCK. + */ + uf_ma = uf; + memcpy(&uf_ma[0], "UF", 2); + if (little_endian()) memcpy(&uf_ma[0], "FU", 2); + uf_ma[1] = 0; /* Not known yet. */ + uf_ma[2] = 0; /* Not known yet. Really, I do. */ + uf_ma[3] = 0; /* Not known yet. */ + uf_ma[4] = 0; /* Not known yet. */ + + uf_ma[6] = 1; + uf_ma[7] = ray_num; + uf_ma[8 ] = 1; + uf_ma[9 ] = sweep_num; + memcpy(&uf_ma[10], r->h.radar_name, 8); + if (little_endian()) swap2(&uf_ma[10], 8/2); + memcpy(&uf_ma[14], r->h.name, 8); + if (little_endian()) swap2(&uf_ma[14], 8/2); + /* Convert decimal lat/lon to d:m:s */ + + if (ray->h.lat != 0.0) { + degree = (int)ray->h.lat; + minute = (int)((ray->h.lat - degree) * 60); + second = (ray->h.lat - degree - minute/60.0) * 3600.0; + } else { + degree = r->h.latd; + minute = r->h.latm; + second = r->h.lats; + } + uf_ma[18] = degree; + uf_ma[19] = minute; + if (second > 0.0) uf_ma[20] = second*64 + 0.5; + else uf_ma[20] = second*64 - 0.5; + + if (ray->h.lon != 0.0) { + degree = (int)ray->h.lon; + minute = (int)((ray->h.lon - degree) * 60); + second = (ray->h.lon - degree - minute/60.0) * 3600.0; + } else { + degree = r->h.lond; + minute = r->h.lonm; + second = r->h.lons; + } + uf_ma[21] = degree; + uf_ma[22] = minute; + if (second > 0.0) uf_ma[23] = second*64 + 0.5; + else uf_ma[23] = second*64 - 0.5; + if (ray->h.alt != 0) + uf_ma[24] = ray->h.alt; + else + uf_ma[24] = r->h.height; + + uf_ma[25] = ray->h.year % 100; /* By definition: not year 2000 compliant. */ + uf_ma[26] = ray->h.month; + uf_ma[27] = ray->h.day; + uf_ma[28] = ray->h.hour; + uf_ma[29] = ray->h.minute; + uf_ma[30] = ray->h.sec; + memcpy(&uf_ma[31], "UT", 2); + if (little_endian()) memcpy(&uf_ma[31], "TU", 2); + if (ray->h.azimuth > 0) uf_ma[32] = ray->h.azimuth*64 + 0.5; + else uf_ma[32] = ray->h.azimuth*64 - 0.5; + uf_ma[33] = ray->h.elev*64 + 0.5; + uf_ma[34] = uf_sweep_mode; + if (ray->h.fix_angle != 0.) + uf_ma[35] = ray->h.fix_angle*64.0 + 0.5; + else uf_ma[35] = sweep[k]->h.elev*64.0 + 0.5; + uf_ma[36] = ray->h.sweep_rate*(360.0/60.0)*64.0 + 0.5; + + the_time = time(NULL); + tm = gmtime(&the_time); + + uf_ma[37] = tm->tm_year % 100; /* Same format as data year */ + uf_ma[38] = tm->tm_mon+1; + uf_ma[39] = tm->tm_mday; + memcpy(&uf_ma[40], "RSL" RSL_VERSION_STR, 8); + if (little_endian()) swap2(&uf_ma[40], 8/2); + uf_ma[44] = (signed short)UF_NO_DATA; + len_ma = 45; + uf_ma[2] = len_ma+1; + /* + * ---- End of MANDITORY HEADER BLOCK. + */ + + /* ---- Begining of OPTIONAL HEADER BLOCK. */ + len_op = 0; + if (q_op) { + q_op = 0; /* Only once. */ + uf_op = uf+len_ma; + memcpy(&uf_op[0], "TRMMGVUF", 8); + if (little_endian()) swap2(&uf_op[0], 8/2); + uf_op[4] = (signed short)UF_NO_DATA; + uf_op[5] = (signed short)UF_NO_DATA; + uf_op[6] = ray->h.hour; + uf_op[7] = ray->h.minute; + uf_op[8] = ray->h.sec; + memcpy(&uf_op[9], "RADAR_UF", 8); + if (little_endian()) swap2(&uf_op[9], 8/2); + uf_op[13] = 2; + len_op = 14; + } + /* ---- End of OPTIONAL HEADER BLOCK. */ + + /* ---- Begining of LOCAL USE HEADER BLOCK. */ + q_lu = 0; + + /* Note: Code within "#ifdef LUHDR_VR_AZ" below is retained for testing + * and is not normally compiled. It was used to deal with azimuth + * differences between DZ and VR in WSR-88D split-cuts, now handled in + * ingest routine. + */ +/* 5/18/2010 Temporarily define LUHDR_VR_AZ until merge_split_cuts is + completed. */ +#define LUHDR_VR_AZ +#ifdef LUHDR_VR_AZ + /* If DZ and VR azimuths are different, store VR azimuth in Local Use + * Header. This is done for WSR-88D split cuts. + */ + if (sweep[DZ_INDEX] && sweep[VR_INDEX]) { + if (sweep[DZ_INDEX]->ray[j] && sweep[VR_INDEX]->ray[j]) { + vr_az = sweep[VR_INDEX]->ray[j]->h.azimuth; + if (sweep[DZ_INDEX]->ray[j]->h.azimuth != vr_az) + q_lu = 1; /* Set to use Local Use Header block. */ + } + } +#endif + len_lu = 0; + if (q_lu) { + /* Store azimuth for WSR-88D VR ray in Local Use Header. */ + uf_lu = uf+len_ma+len_op; + memcpy(&uf_lu[0], "AZ", 2); + if (little_endian()) memcpy(&uf_lu[0], "ZA", 2); + if (vr_az > 0) uf_lu[1] = vr_az*64 + 0.5; + else uf_lu[1] = vr_az*64 - 0.5; + len_lu = 2; + } + /* ---- End of LOCAL USE HEADER BLOCK. */ + + + /* Here is where we loop on each field type. We need to keep + * track of how many FIELD HEADER and FIELD DATA sections, one + * for each field type, we fill. The variable that tracks this + * index into 'uf' is 'current_fh_index'. It is bumped by + * the length of the FIELD HEADER and FIELD DATA for each field + * type encountered. Field types expected are: Reflectivity, + * Velocity, and Spectrum width; this is a typicial list but it + * is not restricted to it. + */ + + for (k=0; kh.nrays && sweep[k]->ray[j]) + ray = sweep[k]->ray[j]; + else + ray = NULL; + else ray = NULL; + + if (ray) { + /* ---- Begining of DATA HEADER. */ + nfield++; + if (q_dh) { + len_dh = 2*true_nvolumes + 3; + uf_dh = uf+len_ma+len_op+len_lu; + uf_dh[0] = nfield; + uf_dh[1] = 1; + uf_dh[2] = nfield; + /* 'nfield' indexes the field number. + * 'k' indexes the particular field from the volume. + * RSL_ftype contains field names and is defined in rsl.h. + */ + if (k > max_field_names-1) { + fprintf(stderr, + "RSL_uf_to_radar: No field name for volume index %d\n", k); + fprintf(stderr,"RSL_ftype must be updated in rsl.h for new field.\n"); + fprintf(stderr,"Quitting now.\n"); + return; + } + memcpy(&uf_dh[3+2*(nfield-1)], RSL_ftype[k], 2); + if (little_endian()) swap2(&uf_dh[3+2*(nfield-1)], 2/2); + if (current_fh_index == 0) current_fh_index = len_ma+len_op+len_lu+len_dh; + uf_dh[4+2*(nfield-1)] = current_fh_index + 1; + } + /* ---- End of DATA HEADER. */ + + /* ---- Begining of FIELD HEADER. */ + if (q_fh) { + uf_fh = uf+current_fh_index; + uf_fh[1] = scale_factor = 100; + uf_fh[2] = ray->h.range_bin1/1000.0; + uf_fh[3] = ray->h.range_bin1 - (1000*uf_fh[2]); + uf_fh[4] = ray->h.gate_size; + uf_fh[5] = ray->h.nbins; + uf_fh[6] = ray->h.pulse_width*(RSL_SPEED_OF_LIGHT/1.0e6); + uf_fh[7] = sweep[k]->h.beam_width*64.0 + 0.5; + uf_fh[8] = sweep[k]->h.beam_width*64.0 + 0.5; + uf_fh[9] = ray->h.frequency*64.0 + 0.5; /* Bandwidth (mHz). */ + uf_fh[10] = 0; /* Horizontal polarization. */ + uf_fh[11] = ray->h.wavelength*64.0*100.0; /* m to cm. */ + uf_fh[12] = ray->h.pulse_count; + memcpy(&uf_fh[13], " ", 2); + uf_fh[14] = (signed short)UF_NO_DATA; + uf_fh[15] = (signed short)UF_NO_DATA; + if (k == DZ_INDEX || k == ZT_INDEX) { + uf_fh[16] = volume[k]->h.calibr_const*100.0 + 0.5; + } + else { + memcpy(&uf_fh[16], " ", 2); + } + if (ray->h.prf != 0) + uf_fh[17] = 1.0/ray->h.prf*1000000.0; /* Pulse repetition time(msec) = 1/prf */ + else + uf_fh[17] = (signed short)UF_NO_DATA; /* Pulse repetition time = 1/prf */ + uf_fh[18] = 16; + if (VR_INDEX == k || VE_INDEX == k) { + uf_fh[19] = scale_factor*ray->h.nyq_vel; + uf_fh[20] = 1; + len_fh = 21; + } else { + len_fh = 19; + } + + uf_fh[0] = current_fh_index + len_fh + 1; + /* ---- End of FIELD HEADER. */ + + /* ---- Begining of FIELD DATA. */ + uf_data = uf+len_fh+current_fh_index; + len_data = ray->h.nbins; + for (m=0; mh.f(ray->range[m]); + if (x == BADVAL || x == RFVAL || x == APFLAG || x == NOECHO) + uf_data[m] = (signed short)UF_NO_DATA; + else + uf_data[m] = scale_factor * x; + } + + current_fh_index += (len_fh+len_data); + } + } + /* ---- End of FIELD DATA. */ + } + /* Fill in some infomation we didn't know. Like, buffer length, + * record number, etc. + */ + rec_num++; + uf_ma[1] = current_fh_index; + uf_ma[3] = len_ma + len_op + 1; + uf_ma[4] = len_ma + len_op + len_lu + 1; + uf_ma[5] = rec_num; + + /* WRITE the UF buffer. */ + rec_len =(int)uf_ma[1]*2; + save_rec_len = rec_len; /* We destroy 'rec_len' when making it + big endian on a little endian machine. */ + if (little_endian()) swap_4_bytes(&rec_len); + (void)fwrite(&rec_len, sizeof(int), 1, fp); + if (little_endian()) swap_uf_buffer(uf); + (void)fwrite(uf, sizeof(char), save_rec_len, fp); + (void)fwrite(&rec_len, sizeof(int), 1, fp); + } /* if (ray) */ + } } } @@ -502,13 +517,13 @@ void RSL_radar_to_uf(Radar *r, char *outfile) { FILE *fp; if (r == NULL) { - fprintf(stderr, "radar_to_uf: radar pointer NULL\n"); - return; + fprintf(stderr, "radar_to_uf: radar pointer NULL\n"); + return; } if ((fp = fopen(outfile, "w")) == NULL) { - perror(outfile); - return; + perror(outfile); + return; } RSL_radar_to_uf_fp(r, fp); @@ -524,13 +539,13 @@ void RSL_radar_to_uf_gzip(Radar *r, char *outfile) { FILE *fp; if (r == NULL) { - fprintf(stderr, "radar_to_uf_gzip: radar pointer NULL\n"); - return; + fprintf(stderr, "radar_to_uf_gzip: radar pointer NULL\n"); + return; } if ((fp = fopen(outfile, "w")) == NULL) { - perror(outfile); - return; + perror(outfile); + return; } fp = compress_pipe(fp); diff --git a/rsl.h b/rsl.h index 358072c..2050806 100644 --- a/rsl.h +++ b/rsl.h @@ -27,7 +27,7 @@ #include "config.h" #endif -#define RSL_VERSION_STR "v1.40" +#define RSL_VERSION_STR "v1.40.3" /**********************************************************************/ /* Configure: Define USE_TWO_BYTE_PRECISION to have RSL store internal*/ @@ -134,12 +134,11 @@ typedef struct { int minute;/* Date for this ray; minute (0-59).*/ float sec; /* Date for this ray; second + fraction of second. */ float unam_rng; /* Unambiguous range. (KM). */ - float azimuth; /* Azimuth angle. (degrees). Must be positive - * 0=North, 90=east, -90/270=west. + float azimuth; /* Azimuth angle (degrees). Must be positive. + * 0=North, 90=east, -90/270=west. * This angle is the mean azimuth for the whole ray. - * Eg. for NSIG the beginning and end azimuths are - * averaged. - */ + * E.g. for NSIG the beginning and end azimuths are averaged. + */ int ray_num; /* Ray no. within elevation scan. */ float elev; /* Elevation angle. (degrees). */ int elev_num; /* Elevation no. within volume scan. */ @@ -151,6 +150,7 @@ typedef struct { float sweep_rate; /* Sweep rate. Full sweeps/min. */ int prf; /* Pulse repetition frequency, in Hz. */ + int prf2; /* Second PRF, for Sigmet dual PRF */ float azim_rate; /* Sweep rate in degrees/sec. */ float fix_angle; /* Elevation angle for the sweep. (degrees). */ float pitch; /* Pitch angle. */ @@ -163,9 +163,9 @@ typedef struct { float lon; /* Longitude (degrees) */ int alt; /* Altitude (m) */ float rvc; /* Radial velocity correction (m/sec) */ - float vel_east; /* Platform velocity to the east (m/sec) */ - float vel_north; /* Platform velocity to the north (m/sec) */ - float vel_up; /* Platform velocity toward up (m/sec) */ + float vel_east; /* Platform velocity to the east (negative for west) (m/sec) */ + float vel_north; /* Platform velocity to the north (negative south) (m/sec) */ + float vel_up; /* Platform velocity toward up (negative down) (m/sec) */ int pulse_count; /* Pulses used in a single dwell time. */ float pulse_width; /* Pulse width (micro-sec). */ float beam_width; /* Beamwidth in degrees. */ @@ -202,7 +202,9 @@ typedef struct { typedef struct { int sweep_num; /* Integer sweep number. This may be redundant, since * this will be the same as the Volume.sweep array index.*/ - float elev; /* Elevation angle (mean) for the sweep. */ + float elev; /* Elevation angle (mean) for the sweep. Value is -999 for + * RHI. */ + float azimuth; /* Azimuth for the sweep (RHI). Value is -999 for PPI. */ float beam_width; /* This is in the ray header too. */ float vert_half_bw; /* Vertical beam width divided by 2 */ float horz_half_bw; /* Horizontal beam width divided by 2 */ @@ -321,6 +323,8 @@ typedef struct { int *data; } Histogram; +enum scan_mode {PPI, RHI}; + typedef struct { int month, day, year; int hour, minute; @@ -355,6 +359,7 @@ typedef struct { int height; /* height of site in meters above sea level*/ int spulse; /* length of short pulse (ns)*/ int lpulse; /* length of long pulse (ns) */ + int scan_mode; /* 0 = PPI, 1 = RHI */ int vcp; /* Volume Coverage Pattern (for WSR-88D only) */ } Radar_header; @@ -374,7 +379,7 @@ typedef struct { * 9 = RH_INDEX = RhoHV: Horz-Vert power corr coeff *10 = PH_INDEX = PhiDP: Differential phase angle *11 = XZ_INDEX = X-band reflectivity. - *12 = CR_INDEX = Corrected DR. + *12 = CD_INDEX = Corrected DR. *13 = MZ_INDEX = DZ mask for 1C-51 HDF. *14 = MR_INDEX = DR mask for 1C-51 HDF. *15 = ZE_INDEX = Edited reflectivity. @@ -387,7 +392,18 @@ typedef struct { *22 = CV_INDEX *23 = AV_INDEX *24 = SQ_INDEX = Signal Quality Index (Sigmet) - */ + *25 = VS_INDEX = Radial Velocity Combined (DORADE) + *26 = VL_INDEX = Radial Velocity Combined (DORADE) + *27 = VG_INDEX = Radial Velocity Combined (DORADE) + *28 = VT_INDEX = Radial Velocity Combined (DORADE) + *29 = NP_INDEX = Normalized Coherent Power (DORADE) + *30 = HC_INDEX = HydroClass (Sigmet) + *31 = VC_INDEX = Radial Velocity Corrected (Sigmet) + *32 = V2_INDEX = Radial Velocity for VCP 121 second Doppler cut. + *33 = S2_INDEX = Spectrum Width for VCP 121 second Doppler cut. + *34 = V3_INDEX = Radial Velocity for VCP 121 third Doppler cut. + *35 = S3_INDEX = Spectrum Width for VCP 121 third Doppler cut. + */ } Radar; @@ -447,19 +463,41 @@ typedef struct { * VE Edited Velocity. VE_INDEX * * KD KDP (deg/km) Differencial Phase KD_INDEX - * [Sigmet, Lassen] + * (Sigmet, Lassen) * * TI TIME (unknown) for MCTEX data. TI_INDEX - * SQ SQI: Signal Quality Index. [Sigmet] SQ_INDEX + * + * SQ SQI: Signal Quality Index. (Sigmet) SQ_INDEX * Decimal fraction from 0 to 1, where 0 * is noise, 1 is noiseless. + * + * VS Radial Velocity, Short PRT (m/s) (DORADE) VS_INDEX + * + * VL Radial Velocity, Long PRT (m/s) (DORADE) VL_INDEX + * + * VG Radial Velocity, combined (m/s) (DORADE) VG_INDEX + * + * VT Radial Velocity, combined (m/s) (DORADE) VT_INDEX + * + * NP Normalized Coherent Power. (DORADE) NP_INDEX + * + * HC HydroClass: enumerated class. (Sigmet) HC_INDEX + * + * VC Radial Velocity corrected for (Sigmet) VC_INDEX + * Nyquist unfolding. */ /* * The number of *_INDEX must never exceed MAX_RADAR_VOLUMES. - * Increase MAX_RADAR_VOLUMES appropriately, for new ingest formats. + * Increase MAX_RADAR_VOLUMES appropriately, for new ingest formats. + * + * Also, when adding new *_INDEXes, you must update the following three arrays + * located near the end of this file: RSL_ftype, RSL_f_list, and RSL_invf_list. + * You also need to modify volume.c, updating the initialization of array + * rsl_qfield by adding a '1' for each new volume index. */ -#define MAX_RADAR_VOLUMES 25 + +#define MAX_RADAR_VOLUMES 42 #define DZ_INDEX 0 #define VR_INDEX 1 @@ -486,6 +524,23 @@ typedef struct { #define CV_INDEX 22 #define AV_INDEX 23 #define SQ_INDEX 24 +#define VS_INDEX 25 +#define VL_INDEX 26 +#define VG_INDEX 27 +#define VT_INDEX 28 +#define NP_INDEX 29 +#define HC_INDEX 30 +#define VC_INDEX 31 +#define V2_INDEX 32 +#define S2_INDEX 33 +#define V3_INDEX 34 +#define S3_INDEX 35 +#define CR_INDEX 36 +#define CC_INDEX 37 +#define PR_INDEX 38 +#define SD_INDEX 39 +#define ZZ_INDEX 40 +#define RD_INDEX 41 /* Prototypes for functions. */ @@ -633,6 +688,8 @@ void RSL_load_blue_table(char *infile); void RSL_print_histogram(Histogram *histogram, int min_range, int max_range, char *filename); void RSL_print_version(); +void RSL_prune_radar_on(); +void RSL_prune_radar_off(); void RSL_radar_to_uf(Radar *r, char *outfile); void RSL_radar_to_uf_gzip(Radar *r, char *outfile); void RSL_radar_verbose_off(void); @@ -745,6 +802,14 @@ float AH_F(Range x); float CV_F(Range x); float AV_F(Range x); float SQ_F(Range x); +float VS_F(Range x); +float VL_F(Range x); +float VG_F(Range x); +float VT_F(Range x); +float NP_F(Range x); +float HC_F(Range x); +float VC_F(Range x); +float SD_F(Range x); Range DZ_INVF(float x); Range VR_INVF(float x); @@ -771,6 +836,14 @@ Range AH_INVF(float x); Range CV_INVF(float x); Range AV_INVF(float x); Range SQ_INVF(float x); +Range VS_INVF(float x); +Range VL_INVF(float x); +Range VG_INVF(float x); +Range VT_INVF(float x); +Range NP_INVF(float x); +Range HC_INVF(float x); +Range VC_INVF(float x); +Range SD_INVF(float x); /* If you like these variables, you can use them in your application @@ -778,23 +851,29 @@ Range SQ_INVF(float x); */ #ifdef USE_RSL_VARS static char *RSL_ftype[] = {"DZ", "VR", "SW", "CZ", "ZT", "DR", - "LR", "ZD", "DM", "RH", "PH", "XZ", - "CD", "MZ", "MD", "ZE", "VE", "KD", - "TI", "DX", "CH", "AH", "CV", "AV", - "SQ"}; + "LR", "ZD", "DM", "RH", "PH", "XZ", + "CD", "MZ", "MD", "ZE", "VE", "KD", + "TI", "DX", "CH", "AH", "CV", "AV", + "SQ", "VS", "VL", "VG", "VT", "NP", + "HC", "VC", "V2", "S2", "V3", "S3", + "CR", "CC", "PR", "SD", "ZZ", "RD"}; static float (*RSL_f_list[])(Range x) = {DZ_F, VR_F, SW_F, CZ_F, ZT_F, DR_F, - LR_F, ZD_F, DM_F, RH_F, PH_F, XZ_F, - CD_F, MZ_F, MD_F, ZE_F, VE_F, KD_F, - TI_F, DX_F, CH_F, AH_F, CV_F, AV_F, - SQ_F}; + LR_F, ZD_F, DM_F, RH_F, PH_F, XZ_F, + CD_F, MZ_F, MD_F, ZE_F, VE_F, KD_F, + TI_F, DX_F, CH_F, AH_F, CV_F, AV_F, + SQ_F, VS_F, VL_F, VG_F, VT_F, NP_F, + HC_F, VC_F, VR_F, SW_F, VR_F, SW_F, + DZ_F, CZ_F, PH_F, SD_F, DZ_F, DZ_F}; static Range (*RSL_invf_list[])(float x) - = {DZ_INVF, VR_INVF, SW_INVF, CZ_INVF, ZT_INVF, DR_INVF, - LR_INVF, ZD_INVF, DM_INVF, RH_INVF, PH_INVF, XZ_INVF, - CD_INVF, MZ_INVF, MD_INVF, ZE_INVF, VE_INVF, KD_INVF, - TI_INVF, DX_INVF, CH_INVF, AH_INVF, CV_INVF, AV_INVF, - SQ_INVF}; + = {DZ_INVF, VR_INVF, SW_INVF, CZ_INVF, ZT_INVF, DR_INVF, + LR_INVF, ZD_INVF, DM_INVF, RH_INVF, PH_INVF, XZ_INVF, + CD_INVF, MZ_INVF, MD_INVF, ZE_INVF, VE_INVF, KD_INVF, + TI_INVF, DX_INVF, CH_INVF, AH_INVF, CV_INVF, AV_INVF, + SQ_INVF, VS_INVF, VL_INVF, VG_INVF, VT_INVF, NP_INVF, + HC_INVF, VC_INVF, VR_INVF, SW_INVF, VR_INVF, SW_INVF, + DZ_INVF, CZ_INVF, PH_INVF, SD_INVF, DZ_INVF, DZ_INVF}; #endif /* Secret routines that are quite useful and useful to developers. */ void radar_load_date_time(Radar *radar); diff --git a/uf_to_radar.c b/uf_to_radar.c index ed51f63..341ec33 100644 --- a/uf_to_radar.c +++ b/uf_to_radar.c @@ -25,12 +25,14 @@ #include #include +/* This allows us to use RSL_ftype, RSL_f_list, RSL_invf_list from rsl.h. */ +#define USE_RSL_VARS #include "rsl.h" extern int radar_verbose_flag; typedef short UF_buffer[16384]; /* Some UF files are bigger than 4096 - * that the UF doc's specify. - */ + * that the UF doc's specify. + */ #define UF_MORE 0 #define UF_DONE 1 @@ -43,10 +45,10 @@ Volume *reset_nsweeps_in_volume(Volume *volume) int i; if (volume == NULL) return NULL; for (i=volume->h.nsweeps; i>0; i--) - if (volume->sweep[i-1] != NULL) { - volume->h.nsweeps = i; - break; - } + if (volume->sweep[i-1] != NULL) { + volume->h.nsweeps = i; + break; + } return volume; } Radar *reset_nsweeps_in_all_volumes(Radar *radar) @@ -54,7 +56,7 @@ Radar *reset_nsweeps_in_all_volumes(Radar *radar) int i; if (radar == NULL) return NULL; for (i=0; ih.nvolumes; i++) - radar->v[i] = reset_nsweeps_in_volume(radar->v[i]); + radar->v[i] = reset_nsweeps_in_volume(radar->v[i]); return radar; } @@ -68,7 +70,7 @@ Volume *copy_sweeps_into_volume(Volume *new_volume, Volume *old_volume) new_volume->h = old_volume->h; new_volume->h.nsweeps = nsweeps; for (i=0; ih.nsweeps; i++) - new_volume->sweep[i] = old_volume->sweep[i]; /* Just copy pointers. */ + new_volume->sweep[i] = old_volume->sweep[i]; /* Just copy pointers. */ /* Free the old sweep array, not the pointers to sweeps. */ free(old_volume->sweep); return new_volume; @@ -79,8 +81,8 @@ void swap2(short *buf, int n) short *end_addr; end_addr = buf + n; while (buf < end_addr) { - swap_2_bytes(buf); - buf++; + swap_2_bytes(buf); + buf++; } } @@ -135,10 +137,11 @@ int uf_into_radar(UF_buffer uf, Radar **the_radar) int current_fh_index; float scale_factor; - int nfields, isweep, ifield, iray, i, m; + int nfields, isweep, ifield, iray, i, j, m; static int pulled_time_from_first_ray = 0; + static int need_scan_mode = 1; char *field_type; /* For printing the field type upon error. */ - short proj_name[4]; + short proj_name[4]; Ray *ray; Sweep *sweep; Radar *radar; @@ -168,8 +171,8 @@ int uf_into_radar(UF_buffer uf, Radar **the_radar) isweep = uf_ma[9] - 1; if (rsl_qsweep != NULL) { - if (isweep > rsl_qsweep_max) return UF_DONE; - if (rsl_qsweep[isweep] == 0) return UF_MORE; + if (isweep > rsl_qsweep_max) return UF_DONE; + if (rsl_qsweep[isweep] == 0) return UF_MORE; } @@ -185,15 +188,29 @@ int uf_into_radar(UF_buffer uf, Radar **the_radar) /* Sticky solution here. */ #else if (radar == NULL) { - radar = RSL_new_radar(MAX_RADAR_VOLUMES); - *the_radar = radar; - pulled_time_from_first_ray = 0; - for (i=0; iv[i] = RSL_new_volume(20); + radar = RSL_new_radar(MAX_RADAR_VOLUMES); + *the_radar = radar; + pulled_time_from_first_ray = 0; + for (i=0; iv[i] = RSL_new_volume(20); } #endif + + if (need_scan_mode) { + /* PPI and RHI are enum constants defined in rsl.h */ + if (uf_ma[34] == 1) radar->h.scan_mode = PPI; + else if (uf_ma[34] == 3) radar->h.scan_mode = RHI; + else { + fprintf(stderr,"Warning: UF sweep mode = %d\n", uf_ma[34]); + fprintf(stderr," Expected 1 or 3 (PPI or RHI)\n"); + fprintf(stderr," Setting radar->h.scan_mode to PPI\n"); + radar->h.scan_mode = PPI; + } + need_scan_mode = 0; + } + /* For LITTLE ENDIAN: * WE "UNSWAP" character strings. Because there are so few of them, * it is easier to catch them here. The entire UF buffer is swapped prior @@ -202,264 +219,223 @@ int uf_into_radar(UF_buffer uf, Radar **the_radar) */ for (i=0; iv[ifield] == NULL) continue; /* Nope. */ - - if (isweep >= radar->v[ifield]->h.nsweeps) { /* Exceeded sweep limit. - * Allocate more sweeps. - * Copy all previous sweeps. - */ - if (radar_verbose_flag) - fprintf(stderr,"Exceeded sweep allocation of %d. Adding 20 more.\n", isweep); - new_volume = RSL_new_volume(radar->v[ifield]->h.nsweeps+20); - new_volume = copy_sweeps_into_volume(new_volume, radar->v[ifield]); - radar->v[ifield] = new_volume; - } - - if (radar->v[ifield]->sweep[isweep] == NULL) { - if (radar_verbose_flag) - fprintf(stderr,"Allocating new sweep for field %d, isweep %d\n", ifield, isweep); - radar->v[ifield]->sweep[isweep] = RSL_new_sweep(1000); - radar->v[ifield]->sweep[isweep]->h.nrays = 0; /* Increment this for each - * ray encountered. - */ - radar->v[ifield]->h.f = f; - radar->v[ifield]->h.invf = invf; - radar->v[ifield]->sweep[isweep]->h.f = f; - radar->v[ifield]->sweep[isweep]->h.invf = invf; - radar->v[ifield]->sweep[isweep]->h.sweep_num = uf_ma[9]; - radar->v[ifield]->sweep[isweep]->h.elev = uf_ma[35] / 64.0; - } - - - - current_fh_index = uf_dh[4+2*i]; - uf_fh = uf + current_fh_index - 1; - sweep = radar->v[ifield]->sweep[isweep]; - iray = sweep->h.nrays; - nbins = uf_fh[5]; - radar->v[ifield]->sweep[isweep]->ray[iray] = RSL_new_ray(nbins); - ray = radar->v[ifield]->sweep[isweep]->ray[iray]; - sweep->h.nrays += 1; - - - if (ray) { - /* - * ---- Beginning of MANDATORY HEADER BLOCK. - */ - ray->h.ray_num = uf_ma[7]; - if (little_endian()) swap2(&uf_ma[10], 8); - memcpy(radar->h.radar_name, &uf_ma[10], 8); - if (little_endian()) swap2(&uf_ma[10], 8/2); - memcpy(radar->h.name, &uf_ma[14], 8); - if (little_endian()) swap2(&uf_ma[14], 8/2); - + if (little_endian()) swap_2_bytes(&uf_dh[3+2*i]); /* Unswap. */ + ifield = -1; + field_type = (char *)&uf_dh[3+2*i]; + for (j=0; jv[ifield] == NULL) continue; /* Nope. */ + + if (isweep >= radar->v[ifield]->h.nsweeps) { /* Exceeded sweep limit. + * Allocate more sweeps. + * Copy all previous sweeps. + */ + if (radar_verbose_flag) + fprintf(stderr,"Exceeded sweep allocation of %d. Adding 20 more.\n", isweep); + new_volume = RSL_new_volume(radar->v[ifield]->h.nsweeps+20); + new_volume = copy_sweeps_into_volume(new_volume, radar->v[ifield]); + radar->v[ifield] = new_volume; + } + + if (radar->v[ifield]->sweep[isweep] == NULL) { + if (radar_verbose_flag) + fprintf(stderr,"Allocating new sweep for field %d, isweep %d\n", ifield, isweep); + radar->v[ifield]->sweep[isweep] = RSL_new_sweep(1000); + radar->v[ifield]->sweep[isweep]->h.nrays = 0; /* Increment this for each + * ray encountered. + */ + radar->v[ifield]->h.f = f; + radar->v[ifield]->h.invf = invf; + radar->v[ifield]->sweep[isweep]->h.f = f; + radar->v[ifield]->sweep[isweep]->h.invf = invf; + radar->v[ifield]->sweep[isweep]->h.sweep_num = uf_ma[9]; + radar->v[ifield]->sweep[isweep]->h.elev = uf_ma[35] / 64.0; + } + + + + current_fh_index = uf_dh[4+2*i]; + uf_fh = uf + current_fh_index - 1; + sweep = radar->v[ifield]->sweep[isweep]; + iray = sweep->h.nrays; + nbins = uf_fh[5]; + radar->v[ifield]->sweep[isweep]->ray[iray] = RSL_new_ray(nbins); + ray = radar->v[ifield]->sweep[isweep]->ray[iray]; + sweep->h.nrays += 1; + + + if (ray) { + /* + * ---- Beginning of MANDATORY HEADER BLOCK. + */ + ray->h.ray_num = uf_ma[7]; + if (little_endian()) swap2(&uf_ma[10], 8); + memcpy(radar->h.radar_name, &uf_ma[10], 8); + if (little_endian()) swap2(&uf_ma[10], 8/2); + memcpy(radar->h.name, &uf_ma[14], 8); + if (little_endian()) swap2(&uf_ma[14], 8/2); + /* All components of lat/lon are the same sign. If not, then * what ever wrote the UF was in error. A simple RSL program * can repair the damage, however, not here. */ - ray->h.lat = uf_ma[18] + uf_ma[19]/60.0 + uf_ma[20]/64.0/3600; - ray->h.lon = uf_ma[21] + uf_ma[22]/60.0 + uf_ma[23]/64.0/3600; - ray->h.alt = uf_ma[24]; - ray->h.year = uf_ma[25]; - if (ray->h.year < 1900) { - ray->h.year += 1900; - if (ray->h.year < 1980) ray->h.year += 100; /* Year >= 2000. */ - } - ray->h.month = uf_ma[26]; - ray->h.day = uf_ma[27]; - ray->h.hour = uf_ma[28]; - ray->h.minute = uf_ma[29]; - ray->h.sec = uf_ma[30]; - ray->h.azimuth = uf_ma[32] / 64.0; - - /* If Local Use Header is present and contains azimuth, use that - * azimuth for VR and SW. This is for WSR-88D, which runs separate - * scans for DZ and VR/SW at the lower elevations, which means DZ - * VR/SW and have different azimuths in the "same" ray. - */ - len_lu = uf_ma[4] - uf_ma[3]; - if (len_lu == 2 && (ifield == VR_INDEX || ifield == SW_INDEX)) { - if (strncmp((char *)uf_lu,"ZA",2) == 0 || - strncmp((char *)uf_lu,"AZ",2) == 0) - ray->h.azimuth = uf_lu[1] / 64.0; - } - if (ray->h.azimuth < 0.) ray->h.azimuth += 360.; /* make it 0 to 360. */ - ray->h.elev = uf_ma[33] / 64.0; - ray->h.elev_num = sweep->h.sweep_num; - ray->h.fix_angle = sweep->h.elev = uf_ma[35] / 64.0; - ray->h.azim_rate = uf_ma[36] / 64.0; - ray->h.sweep_rate = ray->h.azim_rate * (60.0/360.0); - missing_data = uf_ma[44]; - - if (pulled_time_from_first_ray == 0) { - radar->h.height = uf_ma[24]; - radar->h.latd = uf_ma[18]; - radar->h.latm = uf_ma[19]; - radar->h.lats = uf_ma[20] / 64.0; - radar->h.lond = uf_ma[21]; - radar->h.lonm = uf_ma[22]; - radar->h.lons = uf_ma[23] / 64.0; - radar->h.year = ray->h.year; - radar->h.month = ray->h.month; - radar->h.day = ray->h.day; - radar->h.hour = ray->h.hour; - radar->h.minute = ray->h.minute; - radar->h.sec = ray->h.sec; - strcpy(radar->h.radar_type, "uf"); - pulled_time_from_first_ray = 1; - } - /* - * ---- End of MANDATORY HEADER BLOCK. - */ - - /* ---- Optional header used for MCTEX files. */ - /* If this is a MCTEX file, the first 4 words following the - mandatory header contain the string 'MCTEX'. */ - memcpy(proj_name, (short *)(uf + uf_ma[2] - 1), 8); - if (little_endian()) swap2(proj_name, 4); - - - /* ---- Local Use Header (if present) was checked during Mandatory - * Header processing above. - */ - - /* ---- Begining of FIELD HEADER. */ - uf_fh = uf+current_fh_index - 1; - scale_factor = uf_fh[1]; - ray->h.range_bin1 = uf_fh[2] * 1000.0 + uf_fh[3]; - ray->h.gate_size = uf_fh[4]; - - ray->h.nbins = uf_fh[5]; - ray->h.pulse_width = uf_fh[6]/(RSL_SPEED_OF_LIGHT/1.0e6); - - if (strncmp((char *)proj_name, "MCTEX", 5) == 0) /* MCTEX? */ - { - /* The beamwidth values are not correct in Mctex UF files. */ - ray->h.beam_width = 1.0; - sweep->h.beam_width = ray->h.beam_width; - sweep->h.horz_half_bw = ray->h.beam_width/2.0; - sweep->h.vert_half_bw = ray->h.beam_width/2.0; - } - else /* Not MCTEX */ - { - ray->h.beam_width = uf_fh[7] / 64.0; - sweep->h.beam_width = uf_fh[7] / 64.0; - sweep->h.horz_half_bw = uf_fh[7] / 128.0; /* DFF 4/4/95 */ - sweep->h.vert_half_bw = uf_fh[8] / 128.0; /* DFF 4/4/95 */ - } - /* fprintf (stderr, "uf_fh[7] = %d, [8] = %d\n", (int)uf_fh[7], (int)uf_fh[8]); */ - if((int)uf_fh[7] == -32768) { - ray->h.beam_width = 1; - sweep->h.beam_width = 1; - sweep->h.horz_half_bw = .5; - sweep->h.vert_half_bw = .5; - } - - ray->h.frequency = uf_fh[9] / 64.0; - ray->h.wavelength = uf_fh[11] / 64.0 / 100.0; /* cm to m. */ - ray->h.pulse_count = uf_fh[12]; - if (ifield == DZ_INDEX || ifield == ZT_INDEX) { - radar->v[ifield]->h.calibr_const = uf_fh[16] / 100.0; - /* uf value scaled by 100 */ - } - else { - radar->v[ifield]->h.calibr_const = 0.0; - } - if (uf_fh[17] == (short)UF_NO_DATA) x = 0; - else x = uf_fh[17] / 1000000.0; /* PRT in seconds. */ - if (x != 0) { - ray->h.prf = 1/x; - ray->h.unam_rng = RSL_SPEED_OF_LIGHT / (2.0 * ray->h.prf * 1000.0); - } - else { - ray->h.prf = 0.0; - ray->h.unam_rng = 0.0; - } - - if (VR_INDEX == ifield || VE_INDEX == ifield) { - ray->h.nyq_vel = uf_fh[19] / scale_factor; - } - - /* ---- End of FIELD HEADER. */ - - ray->h.f = f; - ray->h.invf = invf; - - /* ---- Begining of FIELD DATA. */ - uf_data = uf+uf_fh[0] - 1; - - len_data = ray->h.nbins; /* Known because of RSL_new_ray. */ - for (m=0; mrange[m] = invf(BADVAL); /* BADVAL */ - else { - if(uf_data[m] == missing_data) - ray->range[m] = invf(NOECHO); /* NOECHO */ - else - ray->range[m] = invf((float)uf_data[m]/scale_factor); - } - } - } + ray->h.lat = uf_ma[18] + uf_ma[19]/60.0 + uf_ma[20]/64.0/3600; + ray->h.lon = uf_ma[21] + uf_ma[22]/60.0 + uf_ma[23]/64.0/3600; + ray->h.alt = uf_ma[24]; + ray->h.year = uf_ma[25]; + if (ray->h.year < 1900) { + ray->h.year += 1900; + if (ray->h.year < 1980) ray->h.year += 100; /* Year >= 2000. */ + } + ray->h.month = uf_ma[26]; + ray->h.day = uf_ma[27]; + ray->h.hour = uf_ma[28]; + ray->h.minute = uf_ma[29]; + ray->h.sec = uf_ma[30]; + ray->h.azimuth = uf_ma[32] / 64.0; + + /* If Local Use Header is present and contains azimuth, use that + * azimuth for VR and SW. This is for WSR-88D, which runs separate + * scans for DZ and VR/SW at the lower elevations, which means DZ + * VR/SW and have different azimuths in the "same" ray. + */ + len_lu = uf_ma[4] - uf_ma[3]; + if (len_lu == 2 && (ifield == VR_INDEX || ifield == SW_INDEX)) { + if (strncmp((char *)uf_lu,"ZA",2) == 0 || + strncmp((char *)uf_lu,"AZ",2) == 0) + ray->h.azimuth = uf_lu[1] / 64.0; + } + if (ray->h.azimuth < 0.) ray->h.azimuth += 360.; /* make it 0 to 360. */ + ray->h.elev = uf_ma[33] / 64.0; + ray->h.elev_num = sweep->h.sweep_num; + ray->h.fix_angle = sweep->h.elev = uf_ma[35] / 64.0; + ray->h.azim_rate = uf_ma[36] / 64.0; + ray->h.sweep_rate = ray->h.azim_rate * (60.0/360.0); + missing_data = uf_ma[44]; + + if (pulled_time_from_first_ray == 0) { + radar->h.height = uf_ma[24]; + radar->h.latd = uf_ma[18]; + radar->h.latm = uf_ma[19]; + radar->h.lats = uf_ma[20] / 64.0; + radar->h.lond = uf_ma[21]; + radar->h.lonm = uf_ma[22]; + radar->h.lons = uf_ma[23] / 64.0; + radar->h.year = ray->h.year; + radar->h.month = ray->h.month; + radar->h.day = ray->h.day; + radar->h.hour = ray->h.hour; + radar->h.minute = ray->h.minute; + radar->h.sec = ray->h.sec; + strcpy(radar->h.radar_type, "uf"); + pulled_time_from_first_ray = 1; + } + /* + * ---- End of MANDATORY HEADER BLOCK. + */ + + /* ---- Optional header used for MCTEX files. */ + /* If this is a MCTEX file, the first 4 words following the + mandatory header contain the string 'MCTEX'. */ + memcpy(proj_name, (short *)(uf + uf_ma[2] - 1), 8); + if (little_endian()) swap2(proj_name, 4); + + + /* ---- Local Use Header (if present) was checked during Mandatory + * Header processing above. + */ + + /* ---- Begining of FIELD HEADER. */ + uf_fh = uf+current_fh_index - 1; + scale_factor = uf_fh[1]; + ray->h.range_bin1 = uf_fh[2] * 1000.0 + uf_fh[3]; + ray->h.gate_size = uf_fh[4]; + + ray->h.nbins = uf_fh[5]; + ray->h.pulse_width = uf_fh[6]/(RSL_SPEED_OF_LIGHT/1.0e6); + + if (strncmp((char *)proj_name, "MCTEX", 5) == 0) /* MCTEX? */ + { + /* The beamwidth values are not correct in Mctex UF files. */ + ray->h.beam_width = 1.0; + sweep->h.beam_width = ray->h.beam_width; + sweep->h.horz_half_bw = ray->h.beam_width/2.0; + sweep->h.vert_half_bw = ray->h.beam_width/2.0; + } + else /* Not MCTEX */ + { + ray->h.beam_width = uf_fh[7] / 64.0; + sweep->h.beam_width = uf_fh[7] / 64.0; + sweep->h.horz_half_bw = uf_fh[7] / 128.0; /* DFF 4/4/95 */ + sweep->h.vert_half_bw = uf_fh[8] / 128.0; /* DFF 4/4/95 */ + } + /* fprintf (stderr, "uf_fh[7] = %d, [8] = %d\n", (int)uf_fh[7], (int)uf_fh[8]); */ + if((int)uf_fh[7] == -32768) { + ray->h.beam_width = 1; + sweep->h.beam_width = 1; + sweep->h.horz_half_bw = .5; + sweep->h.vert_half_bw = .5; + } + + ray->h.frequency = uf_fh[9] / 64.0; + ray->h.wavelength = uf_fh[11] / 64.0 / 100.0; /* cm to m. */ + ray->h.pulse_count = uf_fh[12]; + if (ifield == DZ_INDEX || ifield == ZT_INDEX) { + radar->v[ifield]->h.calibr_const = uf_fh[16] / 100.0; + /* uf value scaled by 100 */ + } + else { + radar->v[ifield]->h.calibr_const = 0.0; + } + if (uf_fh[17] == (short)UF_NO_DATA) x = 0; + else x = uf_fh[17] / 1000000.0; /* PRT in seconds. */ + if (x != 0) { + ray->h.prf = 1/x; + ray->h.unam_rng = RSL_SPEED_OF_LIGHT / (2.0 * ray->h.prf * 1000.0); + } + else { + ray->h.prf = 0.0; + ray->h.unam_rng = 0.0; + } + + if (VR_INDEX == ifield || VE_INDEX == ifield) { + ray->h.nyq_vel = uf_fh[19] / scale_factor; + } + + /* ---- End of FIELD HEADER. */ + + ray->h.f = f; + ray->h.invf = invf; + + /* ---- Begining of FIELD DATA. */ + uf_data = uf+uf_fh[0] - 1; + + len_data = ray->h.nbins; /* Known because of RSL_new_ray. */ + for (m=0; mrange[m] = invf(BADVAL); /* BADVAL */ + else { + if(uf_data[m] == missing_data) + ray->range[m] = invf(NOECHO); /* NOECHO */ + else + ray->range[m] = invf((float)uf_data[m]/scale_factor); + } + } + } } return UF_MORE; } @@ -479,7 +455,7 @@ void swap_uf_buffer(UF_buffer uf) addr_end = uf + sizeof(UF_buffer)/sizeof(short); while (uf < addr_end) - swap_2_bytes(uf++); + swap_2_bytes(uf++); } enum UF_type {NOT_UF, TRUE_UF, TWO_BYTE_UF, FOUR_BYTE_UF}; @@ -496,9 +472,9 @@ enum UF_type {NOT_UF, TRUE_UF, TWO_BYTE_UF, FOUR_BYTE_UF}; Radar *RSL_uf_to_radar_fp(FILE *fp) { union { - char buf[6]; - short sword; - int word; + char buf[6]; + short sword; + int word; } magic; Radar *radar; int nbytes; @@ -521,73 +497,73 @@ Radar *RSL_uf_to_radar_fp(FILE *fp) switch (uf_type) { case FOUR_BYTE_UF: - if (radar_verbose_flag) fprintf(stderr,"UF file with 4 byte FORTRAN record delimeters.\n"); - /* Handle first record specially, since we needed magic information. */ - nbytes = magic.word; - if (little_endian()) swap_4_bytes(&nbytes); - memcpy(uf, &magic.buf[4], 2); - (void)fread(&uf[1], sizeof(char), nbytes-2, fp); - if (little_endian()) swap_uf_buffer(uf); - (void)fread(&nbytes, sizeof(int), 1, fp); - if (uf_into_radar(uf, &radar) == UF_DONE) break; - /* Now the rest of the file. */ - while(fread(&nbytes, sizeof(int), 1, fp) > 0) { - if (little_endian()) swap_4_bytes(&nbytes); - - (void)fread(uf, sizeof(char), nbytes, fp); - if (little_endian()) swap_uf_buffer(uf); - - (void)fread(&nbytes, sizeof(int), 1, fp); - - if (uf_into_radar(uf, &radar) == UF_DONE) break; - } - break; + if (radar_verbose_flag) fprintf(stderr,"UF file with 4 byte FORTRAN record delimeters.\n"); + /* Handle first record specially, since we needed magic information. */ + nbytes = magic.word; + if (little_endian()) swap_4_bytes(&nbytes); + memcpy(uf, &magic.buf[4], 2); + (void)fread(&uf[1], sizeof(char), nbytes-2, fp); + if (little_endian()) swap_uf_buffer(uf); + (void)fread(&nbytes, sizeof(int), 1, fp); + if (uf_into_radar(uf, &radar) == UF_DONE) break; + /* Now the rest of the file. */ + while(fread(&nbytes, sizeof(int), 1, fp) > 0) { + if (little_endian()) swap_4_bytes(&nbytes); + + (void)fread(uf, sizeof(char), nbytes, fp); + if (little_endian()) swap_uf_buffer(uf); + + (void)fread(&nbytes, sizeof(int), 1, fp); + + if (uf_into_radar(uf, &radar) == UF_DONE) break; + } + break; case TWO_BYTE_UF: - if (radar_verbose_flag) fprintf(stderr,"UF file with 2 byte FORTRAN record delimeters.\n"); - /* Handle first record specially, since we needed magic information. */ - sbytes = magic.sword; - if (little_endian()) swap_2_bytes(&sbytes); - memcpy(uf, &magic.buf[2], 4); - (void)fread(&uf[2], sizeof(char), sbytes-4, fp); - if (little_endian()) swap_uf_buffer(uf); - (void)fread(&sbytes, sizeof(short), 1, fp); - uf_into_radar(uf, &radar); - /* Now the rest of the file. */ - while(fread(&sbytes, sizeof(short), 1, fp) > 0) { - if (little_endian()) swap_2_bytes(&sbytes); - - (void)fread(uf, sizeof(char), sbytes, fp); - if (little_endian()) swap_uf_buffer(uf); - - (void)fread(&sbytes, sizeof(short), 1, fp); - - if (uf_into_radar(uf, &radar) == UF_DONE) break; - } - break; + if (radar_verbose_flag) fprintf(stderr,"UF file with 2 byte FORTRAN record delimeters.\n"); + /* Handle first record specially, since we needed magic information. */ + sbytes = magic.sword; + if (little_endian()) swap_2_bytes(&sbytes); + memcpy(uf, &magic.buf[2], 4); + (void)fread(&uf[2], sizeof(char), sbytes-4, fp); + if (little_endian()) swap_uf_buffer(uf); + (void)fread(&sbytes, sizeof(short), 1, fp); + uf_into_radar(uf, &radar); + /* Now the rest of the file. */ + while(fread(&sbytes, sizeof(short), 1, fp) > 0) { + if (little_endian()) swap_2_bytes(&sbytes); + + (void)fread(uf, sizeof(char), sbytes, fp); + if (little_endian()) swap_uf_buffer(uf); + + (void)fread(&sbytes, sizeof(short), 1, fp); + + if (uf_into_radar(uf, &radar) == UF_DONE) break; + } + break; case TRUE_UF: - if (radar_verbose_flag) fprintf(stderr,"UF file with no FORTRAN record delimeters. Good.\n"); - /* Handle first record specially, since we needed magic information. */ - memcpy(&sbytes, &magic.buf[2], 2); /* Record length is in word #2. */ - if (little_endian()) swap_2_bytes(&sbytes); /* # of 2 byte words. */ - - memcpy(uf, &magic.buf[0], 6); - (void)fread(&uf[3], sizeof(short), sbytes-3, fp); - if (little_endian()) swap_uf_buffer(uf); - uf_into_radar(uf, &radar); - /* Now the rest of the file. */ - while(fread(uf, sizeof(short), 2, fp) > 0) { - memcpy(&sbytes, &uf[1], 2); /* Record length is in word #2. */ - if (little_endian()) swap_2_bytes(&sbytes); - - (void)fread(&uf[2], sizeof(short), sbytes-2, fp); /* Have words 1,2. */ - if (little_endian()) swap_uf_buffer(uf); - - if (uf_into_radar(uf, &radar) == UF_DONE) break; - } - break; - + if (radar_verbose_flag) fprintf(stderr,"UF file with no FORTRAN record delimeters. Good.\n"); + /* Handle first record specially, since we needed magic information. */ + memcpy(&sbytes, &magic.buf[2], 2); /* Record length is in word #2. */ + if (little_endian()) swap_2_bytes(&sbytes); /* # of 2 byte words. */ + + memcpy(uf, &magic.buf[0], 6); + (void)fread(&uf[3], sizeof(short), sbytes-3, fp); + if (little_endian()) swap_uf_buffer(uf); + uf_into_radar(uf, &radar); + /* Now the rest of the file. */ + while(fread(uf, sizeof(short), 2, fp) > 0) { + memcpy(&sbytes, &uf[1], 2); /* Record length is in word #2. */ + if (little_endian()) swap_2_bytes(&sbytes); + + (void)fread(&uf[2], sizeof(short), sbytes-2, fp); /* Have words 1,2. */ + if (little_endian()) swap_uf_buffer(uf); + + if (uf_into_radar(uf, &radar) == UF_DONE) break; + } + break; + case NOT_UF: return NULL; break; } radar = reset_nsweeps_in_all_volumes(radar); @@ -618,16 +594,16 @@ Radar *RSL_uf_to_radar(char *infile) radar = NULL; if (infile == NULL) { - int save_fd; - save_fd = dup(0); - fp = fdopen(save_fd, "r"); + int save_fd; + save_fd = dup(0); + fp = fdopen(save_fd, "r"); } else if ((fp = fopen(infile, "r")) == NULL) { - perror(infile); - return radar; + perror(infile); + return radar; } fp = uncompress_pipe(fp); /* Transparently gunzip. */ radar = RSL_uf_to_radar_fp(fp); rsl_pclose(fp); - + return radar; } diff --git a/volume.c b/volume.c index 2ff6aa9..a30368b 100644 --- a/volume.c +++ b/volume.c @@ -114,7 +114,7 @@ extern int radar_verbose_flag; float DZ_F(Range x) { if (x >= F_OFFSET) /* This test works when Range is unsigned. */ - return (((float)x-F_OFFSET)/F_FACTOR - F_DZ_RANGE_OFFSET); /* Default wsr88d. */ + return (((float)x-F_OFFSET)/F_FACTOR - F_DZ_RANGE_OFFSET); /* Default wsr88d. */ if (x == 0) return BADVAL; if (x == 1) return RFVAL; if (x == 2) return APFLAG; @@ -125,9 +125,9 @@ 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 - F_VR_OFFSET); /* Default wsr88d coding. */ - /* fprintf(stderr, "x=%d, val=%f\n", x, val); */ - return val; + val = (((float)x-F_OFFSET)/F_FACTOR - F_VR_OFFSET); /* Default wsr88d coding. */ + /* fprintf(stderr, "x=%d, val=%f\n", x, val); */ + return val; } if (x == 0) return BADVAL; if (x == 1) return RFVAL; @@ -139,8 +139,8 @@ 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 - F_DR_OFFSET); - return val; + val = (((float)x-F_OFFSET)/F_DR_FACTOR - F_DR_OFFSET); + return val; } if (x == 0) return BADVAL; if (x == 1) return RFVAL; @@ -151,7 +151,7 @@ float DR_F(Range x) { /* Differential reflectivity */ float LR_F(Range x) {/* From MCTEX */ if (x >= F_OFFSET) /* This test works when Range is unsigned. */ - return (float) (x - 250.)/6.; + return (float) (x - 250.)/6.; if (x == 0) return BADVAL; if (x == 1) return RFVAL; if (x == 2) return APFLAG; @@ -228,13 +228,23 @@ float KD_F(Range x) return (x-32768.)/100.; } -float NP_F(Range x) { /* Normalized Coherent Power (DORADE) */ +/* Normalized Coherent Power (DORADE) */ +float NP_F(Range x) +{ + if (x == 0) return BADVAL; + return (float)(x - 1) / 100.; +} + +/* Standard Deviation (for Dual-pole QC testing.) */ +float SD_F(Range x) +{ if (x == 0) return BADVAL; return (float)x / 100.; } /* Signal Quality Index */ -float SQ_F(Range x) { +float SQ_F(Range x) +{ if (x == 0) return BADVAL; return (float)(x-1) / 65533.; } @@ -360,28 +370,34 @@ Range PH_INVF(float x) { /* 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. + return (Range)(x * 100. + 32768. + 0.5); +/****** Old code for 1-byte Sigmet native data format commented-out: if (x == RFVAL) return (Range)1; if (x == APFLAG) return (Range)2; if (x == NOECHO) return (Range)3; if (rsl_kdp_wavelen == 0.0) return (Range)0; if (x < 0) { - x = 127 - - 126 * (log((double)-x) - log((double)(0.25/rsl_kdp_wavelen))) / - log((double)600.0) + - 0.5; + x = 127 - + 126 * (log((double)-x) - log((double)(0.25/rsl_kdp_wavelen))) / + log((double)600.0) + + 0.5; } else if (x > 0) { - x = 129 + - 126 * (log((double)x) - log((double)0.25/rsl_kdp_wavelen)) / - log((double)600.0) + - 0.5; + x = 129 + + 126 * (log((double)x) - log((double)0.25/rsl_kdp_wavelen)) / + log((double)600.0) + + 0.5; } else { - x = 128; + x = 128; } x += F_OFFSET; ******/ - return (Range)(x * 100. + 32768. + 0.5); - +} + +/* Standard Deviation (for Dual-pole QC testing.) */ +Range SD_INVF(float x) +{ + if (x == BADVAL) return (Range)0; + return (Range)(x * 100.); } /* Signal Quality Index */ @@ -391,11 +407,11 @@ Range SQ_INVF(float x) return (Range)(x * 65533. + 1. +.5); } - -Range NP_INVF(float x) /* Normalized Coherent Power (DORADE) */ +/* Normalized Coherent Power (DORADE) */ +Range NP_INVF(float x) { if (x == BADVAL) return (0); - return (Range)(x * 100.); + return (Range)(x * 100. + 1.); } Range TI_INVF(float x) /* MCTEX */ @@ -497,9 +513,9 @@ typedef struct { #define STATIC STATIC int RSL_max_sweeps = 0; /* Initial allocation for sweep_list. - * RSL_new_sweep will allocate the space first - * time around. - */ + * RSL_new_sweep will allocate the space first + * time around. + */ STATIC int RSL_nsweep_addr = 0; /* A count of sweeps in the table. */ STATIC Sweep_list *RSL_sweep_list = NULL; STATIC int RSL_nextents = 0; @@ -516,7 +532,7 @@ void FREE_HASH_TABLE(Hash_table *table) int i; if (table == NULL) return; for (i=0; inindexes; i++) - FREE_HASH_NODE(table->indexes[i]); /* A possible linked list of Rays. */ + FREE_HASH_NODE(table->indexes[i]); /* A possible linked list of Rays. */ free(table->indexes); free(table); } @@ -527,7 +543,7 @@ void REMOVE_SWEEP(Sweep *s) int j; /* Find where it goes, split the list and slide the tail down one. */ for (i=0; i= RSL_max_sweeps) { /* Current list is too small. */ - RSL_nextents++; - new_list = (Sweep_list *) calloc(100*RSL_nextents, sizeof(Sweep_list)); - if (new_list == NULL) { - perror("INSERT_SWEEP"); - exit(2); - } + RSL_nextents++; + new_list = (Sweep_list *) calloc(100*RSL_nextents, sizeof(Sweep_list)); + if (new_list == NULL) { + perror("INSERT_SWEEP"); + exit(2); + } /* Copy the old list to the new one. */ - for (i=0; ii; j--) - RSL_sweep_list[j] = RSL_sweep_list[j-1]; + RSL_sweep_list[j] = RSL_sweep_list[j-1]; RSL_sweep_list[i].s_addr = s; RSL_sweep_list[i].hash = NULL; @@ -581,7 +597,7 @@ int SWEEP_INDEX(Sweep *s) /* Simple linear search; but this will be a binary search. */ int i; for (i=0; iray = (Ray **) calloc(max_rays, sizeof(Ray*)); if (s->ray == NULL) perror("RSL_new_sweep, Ray*"); s->h.nrays = max_rays; /* A default setting. */ + s->h.elev = -999.; + s->h.azimuth = -999.; return s; } @@ -633,7 +651,7 @@ Sweep *RSL_clear_sweep(Sweep *s) int i; if (s == NULL) return s; for (i=0; ih.nrays; i++) { - RSL_clear_ray(s->ray[i]); + RSL_clear_ray(s->ray[i]); } return s; } @@ -642,7 +660,7 @@ Volume *RSL_clear_volume(Volume *v) int i; if (v == NULL) return v; for (i=0; ih.nsweeps; i++) { - RSL_clear_sweep(v->sweep[i]); + RSL_clear_sweep(v->sweep[i]); } return v; } @@ -664,7 +682,7 @@ void RSL_free_sweep(Sweep *s) int i; if (s == NULL) return; for (i=0; ih.nrays; i++) { - RSL_free_ray(s->ray[i]); + RSL_free_ray(s->ray[i]); } if (s->ray) free(s->ray); REMOVE_SWEEP(s); /* Remove from internal Sweep list. */ @@ -676,9 +694,9 @@ void RSL_free_volume(Volume *v) if (v == NULL) return; for (i=0; ih.nsweeps; i++) - { - RSL_free_sweep(v->sweep[i]); - } + { + RSL_free_sweep(v->sweep[i]); + } if (v->sweep) free(v->sweep); free(v); } @@ -711,7 +729,7 @@ Sweep *RSL_copy_sweep(Sweep *s) n_sweep->h = s->h; for (i=0; ih.nrays; i++) { - n_sweep->ray[i] = RSL_copy_ray(s->ray[i]); + n_sweep->ray[i] = RSL_copy_ray(s->ray[i]); } return n_sweep; } @@ -728,7 +746,7 @@ Volume *RSL_copy_volume(Volume *v) new_vol->h = v->h; for (i=0; ih.nsweeps; i++) { - new_vol->sweep[i] = RSL_copy_sweep(v->sweep[i]); + new_vol->sweep[i] = RSL_copy_sweep(v->sweep[i]); } return new_vol; } @@ -817,9 +835,9 @@ Ray *RSL_get_next_ccwise_ray(Sweep *s, Ray *ray) double cwise_angle_diff(float x,float y) { /* Returns the clockwise angle difference of x to y. - * If x = 345 and y = 355 return 10. - * If x = 345 and y = 335 return 350 - */ + * If x = 345 and y = 355 return 10. + * If x = 345 and y = 335 return 350 + */ double d; d = (double)(y - x); @@ -836,9 +854,9 @@ double cwise_angle_diff(float x,float y) double ccwise_angle_diff(float x,float y) { /* Returns the counterclockwise angle differnce of x to y. - * If x = 345 and y = 355 return 350. - * If x = 345 and y = 335 return 10 - */ + * If x = 345 and y = 355 return 350. + * If x = 345 and y = 335 return 10 + */ double d; d = (double)(x - y); @@ -862,48 +880,48 @@ Azimuth_hash *the_closest_hash(Azimuth_hash *hash, float ray_angle) if (hash == NULL) return NULL; /* Set low pointer to hash index with ray angle just below - * requested angle and high pointer to just above requested - * angle. - */ + * requested angle and high pointer to just above requested + * angle. + */ /* set low and high pointers to initial search locations*/ low = hash; high = hash->ray_high; /* Search until clockwise angle to high is less then clockwise - * angle to low. - */ + * angle to low. + */ clow = cwise_angle_diff(ray_angle,low->ray->h.azimuth); chigh = cwise_angle_diff(ray_angle,high->ray->h.azimuth); cclow = ccwise_angle_diff(ray_angle,low->ray->h.azimuth); while((chigh > clow) && (clow != 0)) - { - if (clow < cclow) - { - low = low->ray_low; - high = low->ray_high; /* Not the same low as line before ! */ - } - else - { - low = low->ray_high; - high = low->ray_high; /* Not the same low as line before ! */ - } - - clow = cwise_angle_diff(ray_angle,low->ray->h.azimuth); - chigh = cwise_angle_diff(ray_angle,high->ray->h.azimuth); - cclow = ccwise_angle_diff(ray_angle,low->ray->h.azimuth); - } + { + if (clow < cclow) + { + low = low->ray_low; + high = low->ray_high; /* Not the same low as line before ! */ + } + else + { + low = low->ray_high; + high = low->ray_high; /* Not the same low as line before ! */ + } + + clow = cwise_angle_diff(ray_angle,low->ray->h.azimuth); + chigh = cwise_angle_diff(ray_angle,high->ray->h.azimuth); + cclow = ccwise_angle_diff(ray_angle,low->ray->h.azimuth); + } if(chigh <= cclow) - { - return high; - } + { + return high; + } else - { - return low; - } + { + return low; + } } @@ -1025,8 +1043,8 @@ int hash_bin(Hash_table *table,float angle) * why bother? */ while(table->indexes[hash] == NULL) { - hash++; - if(hash >= table->nindexes) hash = 0; + hash++; + if(hash >= table->nindexes) hash = 0; } return hash; @@ -1038,13 +1056,13 @@ Hash_table *hash_table_for_sweep(Sweep *s) i = SWEEP_INDEX(s); if (i==-1) { /* Obviously, an unregistered sweep. Most likely the - * result of pointer assignments. - */ - i = INSERT_SWEEP(s); + * result of pointer assignments. + */ + i = INSERT_SWEEP(s); } if (RSL_sweep_list[i].hash == NULL) { /* First time. Construct the table. */ - RSL_sweep_list[i].hash = construct_sweep_hash_table(s); + RSL_sweep_list[i].hash = construct_sweep_hash_table(s); } return RSL_sweep_list[i].hash; @@ -1061,7 +1079,7 @@ Ray *RSL_get_closest_ray_from_sweep(Sweep *s,float ray_angle, float limit) /* * Return closest Ray in Sweep within limit (angle) specified * in parameter list. Assume PPI mode. - */ + */ int hindex; Hash_table *hash_table; Azimuth_hash *closest; @@ -1078,8 +1096,8 @@ Ray *RSL_get_closest_ray_from_sweep(Sweep *s,float ray_angle, float limit) closest = the_closest_hash(hash_table->indexes[hindex],ray_angle); /* Is closest ray within limit parameter ? If - * so return ray, else return NULL. - */ + * so return ray, else return NULL. + */ close_diff = angle_diff(ray_angle,closest->ray->h.azimuth); @@ -1134,13 +1152,13 @@ float RSL_get_value_from_ray(Ray *ray, float r) if (ray == NULL) return BADVAL; if(ray->h.gate_size == 0) - { - if(radar_verbose_flag) - { - fprintf(stderr,"RSL_get_value_from_ray: ray->h.gate_size == 0\n"); - } - return BADVAL; - } + { + if(radar_verbose_flag) + { + fprintf(stderr,"RSL_get_value_from_ray: ray->h.gate_size == 0\n"); + } + return BADVAL; + } /* range_bin1 is range to center of first bin */ bin_index = (int)(((rm - ray->h.range_bin1)/ray->h.gate_size) + 0.5); @@ -1254,10 +1272,10 @@ Ray *RSL_get_ray_above(Volume *v, Ray *current_ray) i++; while( i < v->h.nsweeps) - { - if(v->sweep[i] != NULL) break; - i++; - } + { + if(v->sweep[i] != NULL) break; + i++; + } if(i >= v->h.nsweeps) return NULL; @@ -1283,10 +1301,10 @@ Ray *RSL_get_ray_below(Volume *v, Ray *current_ray) i--; while( i >= 0) - { - if(v->sweep[i] != NULL) break; - i--; - } + { + if(v->sweep[i] != NULL) break; + i--; + } if(i < 0) return NULL; @@ -1332,13 +1350,13 @@ Ray *RSL_get_first_ray_of_sweep(Sweep *s) smallest_ray_num = 9999999; if (s == NULL) return r; for (i=0; ih.nrays; i++) - if (s->ray[i]) { - if (s->ray[i]->h.ray_num <= 1) return s->ray[i]; - if (s->ray[i]->h.ray_num < smallest_ray_num) { - r = s->ray[i]; - smallest_ray_num = r->h.ray_num; - } - } + if (s->ray[i]) { + if (s->ray[i]->h.ray_num <= 1) return s->ray[i]; + if (s->ray[i]->h.ray_num < smallest_ray_num) { + r = s->ray[i]; + smallest_ray_num = r->h.ray_num; + } + } return r; } @@ -1361,7 +1379,7 @@ Sweep *RSL_get_first_sweep_of_volume(Volume *v) int i; if (v == NULL) return NULL; for (i=0; ih.nsweeps; i++) - if (RSL_get_first_ray_of_sweep(v->sweep[i])) return v->sweep[i]; + if (RSL_get_first_ray_of_sweep(v->sweep[i])) return v->sweep[i]; return NULL; } @@ -1384,6 +1402,8 @@ 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 }; @@ -1438,26 +1458,26 @@ void RSL_select_fields(char *field_type, ...) if (radar_verbose_flag) fprintf(stderr,"Selected fields for ingest:"); while (c_field) { - /* CHECK EACH FIELD. This is a fancier case statement than C provides. */ - if (radar_verbose_flag) fprintf(stderr," %s", c_field); - if (strcasecmp(c_field, "all") == 0) { - for (i=0; i> specified.\n", c_field); - } - } - c_field = va_arg(ap, char *); + /* CHECK EACH FIELD. This is a fancier case statement than C provides. */ + if (radar_verbose_flag) fprintf(stderr," %s", c_field); + if (strcasecmp(c_field, "all") == 0) { + for (i=0; i> specified.\n", c_field); + } + } + c_field = va_arg(ap, char *); } if (radar_verbose_flag) fprintf(stderr,"\n"); @@ -1497,15 +1517,15 @@ int rsl_query_field(char *c_field) RSL_invf_list[0] = RSL_invf_list[0]; for (i=0; i> specified.\n", c_field); + * properly written ingest code. + */ + fprintf(stderr, "rsl_query_field: Invalid field name <<%s>> specified.\n", c_field); } /* 'i' is the index. Is it set? */ @@ -1515,8 +1535,8 @@ int rsl_query_field(char *c_field) /* Could be static and force use of 'rsl_query_sweep' */ int *rsl_qsweep = NULL; /* If NULL, then read all sweeps. Otherwise, - * read what is on the list. - */ + * read what is on the list. + */ #define RSL_MAX_QSWEEP 500 /* It'll be rediculious to have more. :-) */ int rsl_qsweep_max = RSL_MAX_QSWEEP; @@ -1555,35 +1575,35 @@ void RSL_read_these_sweeps(char *csweep, ...) if (radar_verbose_flag) fprintf(stderr,"Selected sweeps for ingest:"); - for (;c_sweep; c_sweep = va_arg(ap, char *)) - { - /* CHECK EACH FIELD. This is a fancier case statement than C provides. */ - if (radar_verbose_flag) fprintf(stderr," %s", c_sweep); - if (strcasecmp(c_sweep, "all") == 0) { - for (i=0; i RSL_MAX_QSWEEP) { - if (radar_verbose_flag) fprintf(stderr,"\nRSL_read_these_sweeps: parameter %s not in [0,%d). Ignoring.\n", c_sweep, RSL_MAX_QSWEEP); - continue; - } - - if (isweep > rsl_qsweep_max) rsl_qsweep_max = isweep; - rsl_qsweep[isweep] = 1; - } + for (;c_sweep; c_sweep = va_arg(ap, char *)) + { + /* CHECK EACH FIELD. This is a fancier case statement than C provides. */ + if (radar_verbose_flag) fprintf(stderr," %s", c_sweep); + if (strcasecmp(c_sweep, "all") == 0) { + for (i=0; i RSL_MAX_QSWEEP) { + if (radar_verbose_flag) fprintf(stderr,"\nRSL_read_these_sweeps: parameter %s not in [0,%d). Ignoring.\n", c_sweep, RSL_MAX_QSWEEP); + continue; + } + + if (isweep > rsl_qsweep_max) rsl_qsweep_max = isweep; + rsl_qsweep[isweep] = 1; + } } if (radar_verbose_flag) fprintf(stderr,"\n"); @@ -1639,9 +1659,9 @@ void RSL_add_dbz_offset_to_ray(Ray *r, float dbz_offset) if (r == NULL) return; for (ibin=0; ibinh.nbins; ibin++) { - val = r->h.f(r->range[ibin]); - if ( val >= (float)NOECHO ) continue; /* Invalid value */ - r->range[ibin] = r->h.invf(val + dbz_offset); + val = r->h.f(r->range[ibin]); + if ( val >= (float)NOECHO ) continue; /* Invalid value */ + r->range[ibin] = r->h.invf(val + dbz_offset); } } @@ -1655,7 +1675,7 @@ void RSL_add_dbz_offset_to_sweep(Sweep *s, float dbz_offset) int iray; if (s == NULL) return; for (iray=0; irayh.nrays; iray++) - RSL_add_dbz_offset_to_ray(s->ray[iray], dbz_offset); + RSL_add_dbz_offset_to_ray(s->ray[iray], dbz_offset); } /*********************************************************************/ @@ -1668,5 +1688,5 @@ void RSL_add_dbz_offset_to_volume(Volume *v, float dbz_offset) int isweep; if (v == NULL) return; for (isweep=0; isweeph.nsweeps; isweep++) - RSL_add_dbz_offset_to_sweep(v->sweep[isweep], dbz_offset); + RSL_add_dbz_offset_to_sweep(v->sweep[isweep], dbz_offset); } diff --git a/wsr88d.c b/wsr88d.c index 3339851..eeadeb8 100644 --- a/wsr88d.c +++ b/wsr88d.c @@ -700,6 +700,9 @@ int wsr88d_get_volume_coverage(Wsr88d_ray *ray) if (ray->vol_cpat == 31) return 31; if (ray->vol_cpat == 32) return 32; if (ray->vol_cpat == 121) return 121; + if (ray->vol_cpat == 211) return 211; + if (ray->vol_cpat == 212) return 212; + if (ray->vol_cpat == 221) return 221; return 0; } diff --git a/wsr88d_locations.dat b/wsr88d_locations.dat index d044206..0df1702 100644 --- a/wsr88d_locations.dat +++ b/wsr88d_locations.dat @@ -108,6 +108,7 @@ 53833 KOHX NASHVILLE TN 36 14 50 -86 33 45 176 94703 KOKX NEW_YORK_CITY NY 40 51 56 -72 51 50 26 4106 KOTX SPOKANE WA 47 40 49 -117 37 36 727 +3948 KOUN NORMAN OK 35 14 10 -97 27 44 390 3816 KPAH PADUCAH KY 37 4 6 -88 46 19 119 4832 KPBZ PITTSBURGH PA 40 31 54 -80 13 6 361 24155 KPDT PENDLETON OR 45 41 26 -118 51 10 462 diff --git a/wsr88d_m31.c b/wsr88d_m31.c index 3c260a7..4373df3 100644 --- a/wsr88d_m31.c +++ b/wsr88d_m31.c @@ -28,6 +28,7 @@ #include "rsl.h" #include "wsr88d.h" +#include /* Data descriptions in the following data structures are from the "Interface * Control Document for the RDA/RPG", Build 10.0 Draft, WSR-88D Radar @@ -35,17 +36,15 @@ */ typedef struct { - short rpg[6]; /* Unused. Not really part of message header, but is - * inserted by RPG Communications Manager for each message. - */ - unsigned short msg_size; /* for this segment, in halfwords */ + short rpg[6]; /* 12 bytes inserted by RPG Communications Mgr. Ignored. */ + unsigned short msg_size; /* Message size for this segment, in halfwords */ unsigned char channel; /* RDA Redundant Channel */ - unsigned char msg_type; /* Message Type */ - unsigned short id_seq; /* I.d. Seq = 0 to 7FFF, then roll over to 0 */ - unsigned short msg_date; /* modified Julian date from 1/1/70 */ - unsigned int msg_time; /* packet generation time in ms past midnite */ - unsigned short num_segs; - unsigned short seg_num; + unsigned char msg_type; /* Message type. For example, 31 */ + unsigned short id_seq; /* Msg seq num = 0 to 7FFF, then roll over to 0 */ + unsigned short msg_date; /* Modified Julian date from 1/1/70 */ + unsigned int msg_time; /* Packet generation time in ms past midnight */ + unsigned short num_segs; /* Number of segments for this message */ + unsigned short seg_num; /* Number of this segment */ } Wsr88d_msg_hdr; typedef struct { @@ -91,32 +90,18 @@ typedef struct { float offset; } Data_moment_hdr; -typedef struct { - Data_moment_hdr data_hdr; - /* TODO: data type will need to be changed to "void" to handle Dual Pole: - * some data are in 2-byte words. - */ - unsigned char *data; -} Data_moment; +#define MAX_RADIAL_LENGTH 14288 typedef struct { Ray_header_m31 ray_hdr; - short unamb_rng; - short nyq_vel; - /* Pointers to data */ - /* TODO: Replace data moment pointers with single pointer to radial. - * Can simply make unsigned char array as in MSG1 version. - * Maximum radial length is 14288 bytes. - */ - Data_moment *ref; - Data_moment *vel; - Data_moment *sw; + float unamb_rng; + float nyq_vel; + unsigned char data[MAX_RADIAL_LENGTH]; } Wsr88d_ray_m31; -#define MAXRAYS_M31 800 -#define MAXSWEEPS 20 -enum {START_OF_ELEV, INTERMED_RADIAL, END_OF_ELEV, BEGIN_VOS, END_VOS}; +enum radial_status {START_OF_ELEV, INTERMED_RADIAL, END_OF_ELEV, BEGIN_VOS, + END_VOS}; void wsr88d_swap_m31_hdr(Wsr88d_msg_hdr *msghdr) @@ -132,7 +117,7 @@ void wsr88d_swap_m31_hdr(Wsr88d_msg_hdr *msghdr) void wsr88d_swap_m31_ray(Ray_header_m31 *wsr88d_ray) { - int *fullword; + int *data_ptr; swap_4_bytes(&wsr88d_ray->ray_time); swap_2_bytes(&wsr88d_ray->ray_date); @@ -141,13 +126,13 @@ void wsr88d_swap_m31_ray(Ray_header_m31 *wsr88d_ray) swap_2_bytes(&wsr88d_ray->radial_len); swap_4_bytes(&wsr88d_ray->elev); swap_2_bytes(&wsr88d_ray->data_block_count); - fullword = (int *) &wsr88d_ray->dbptr_vol_const; - for (; fullword <= (int *) &wsr88d_ray->dbptr_rho; fullword++) - swap_4_bytes(fullword); + data_ptr = (int *) &wsr88d_ray->dbptr_vol_const; + for (; data_ptr <= (int *) &wsr88d_ray->dbptr_rho; data_ptr++) + swap_4_bytes(data_ptr); } -void wsr88d_swap_data_moment(Data_moment_hdr *this_field) +void wsr88d_swap_data_hdr(Data_moment_hdr *this_field) { short *halfword; halfword = (short *) &this_field->ngates; @@ -158,25 +143,6 @@ void wsr88d_swap_data_moment(Data_moment_hdr *this_field) } -void testprt(Ray_header_m31 ray_hdr) -{ - /* For testing: Print some values from data block header. */ - printf("\nray_time: %d\n",ray_hdr.ray_time); - printf("ray_date: %d\n",ray_hdr.ray_date); - printf("azm_num: %d\n",ray_hdr.azm_num); - printf("azm: %f\n",ray_hdr.azm); - printf("radial_len: %d\n",ray_hdr.radial_len); - printf("elev: %f\n",ray_hdr.elev); - printf("data block count: %d\n",ray_hdr.data_block_count); - printf("dbptr_vol_const: %d\n",ray_hdr.dbptr_vol_const); - printf("dbptr_elev_const: %d\n",ray_hdr.dbptr_elev_const); - printf("dbptr_radial_const: %d\n",ray_hdr.dbptr_radial_const); - printf("dbptr_ref: %d\n",ray_hdr.dbptr_ref); - printf("dbptr_vel: %d\n",ray_hdr.dbptr_vel); - printf("dbptr_sw: %d\n",ray_hdr.dbptr_sw); -} - - float wsr88d_get_angle(short bitfield) { short mask = 1; @@ -273,31 +239,39 @@ void wsr88d_get_vcp_data(short *msgtype5) } -Data_moment *read_data_moment(unsigned char *radial, int begin_data_block) +void get_wsr88d_unamb_and_nyq_vel(Wsr88d_ray_m31 *wsr88d_ray, float *unamb_rng, + float *nyq_vel) { - Data_moment *this_field; - unsigned char* data; - unsigned char* dptr; - int hdr_size, n; - - /* - printf("read_data_moment: "); - printf("sizeof(Data_moment) = %d\n", sizeof(Data_moment)); - printf("sizeof(Data_moment_hdr) = %d\n", sizeof(Data_moment_hdr)); - */ - this_field = malloc(sizeof(Data_moment)); - hdr_size = sizeof(Data_moment_hdr); - /* printf("hdr_size = %d\n", hdr_size); */ - - dptr = radial + begin_data_block; - memcpy(&this_field->data_hdr, dptr, hdr_size); - dptr += hdr_size; - - if (little_endian()) wsr88d_swap_data_moment(&this_field->data_hdr); - data = (unsigned char *) malloc(this_field->data_hdr.ngates); - memcpy(data, dptr, this_field->data_hdr.ngates); - this_field->data = data; - return this_field; + int dbptr, found, ray_hdr_len; + short nyq_vel_sh, unamb_rng_sh; + + found = 0; + ray_hdr_len = sizeof(wsr88d_ray->ray_hdr); + dbptr = wsr88d_ray->ray_hdr.dbptr_radial_const - ray_hdr_len; + if (strncmp((char *) &wsr88d_ray->data[dbptr], "RRAD", 4) == 0) found = 1; + else { + dbptr = wsr88d_ray->ray_hdr.dbptr_elev_const - ray_hdr_len; + if (strncmp((char *) &wsr88d_ray->data[dbptr], "RRAD", 4) == 0) + found = 1; + else { + dbptr = wsr88d_ray->ray_hdr.dbptr_vol_const - ray_hdr_len; + if (strncmp((char *) &wsr88d_ray->data[dbptr], "RRAD", 4) == 0) + found = 1; + } + } + if (found) { + memcpy(&unamb_rng_sh, &wsr88d_ray->data[dbptr+6], 2); + memcpy(&nyq_vel_sh, &wsr88d_ray->data[dbptr+16], 2); + if (little_endian()) { + swap_2_bytes(&unamb_rng_sh); + swap_2_bytes(&nyq_vel_sh); + } + *unamb_rng = unamb_rng_sh / 10.; + *nyq_vel = nyq_vel_sh / 100.; + } else { + *unamb_rng = 0.; + *nyq_vel = 0.; + } } @@ -305,8 +279,8 @@ int read_wsr88d_ray_m31(Wsr88d_file *wf, int msg_size, Wsr88d_ray_m31 *wsr88d_ray) { Ray_header_m31 ray_hdr; - int n, begin_data_block, bytes_to_read, ray_hdr_len; - unsigned char *radial; + int n, bytes_to_read, ray_hdr_len; + float nyq_vel, unamb_rng; /* Read ray data header block */ n = fread(&wsr88d_ray->ray_hdr, sizeof(Ray_header_m31), 1, wf->fptr); @@ -314,170 +288,33 @@ int read_wsr88d_ray_m31(Wsr88d_file *wf, int msg_size, fprintf(stderr,"read_wsr88d_ray_m31: Read failed.\n"); return 0; } + + /* Byte swap header if needed. */ if (little_endian()) wsr88d_swap_m31_ray(&wsr88d_ray->ray_hdr); ray_hdr = wsr88d_ray->ray_hdr; - - /* - if (ray_hdr.azm_num == 1) testprt(wsr88d_ray->ray_hdr); - */ - ray_hdr_len = sizeof(ray_hdr); - /* Read in radial - * Use header_size offset with data pointers - * Copy data into structures - */ + /* Read data portion of radial. */ bytes_to_read = msg_size - ray_hdr_len; - radial = (unsigned char *) malloc(bytes_to_read); - n = fread(radial, bytes_to_read, 1, wf->fptr); + n = fread(wsr88d_ray->data, bytes_to_read, 1, wf->fptr); if (n < 1) { fprintf(stderr,"read_wsr88d_ray_m31: Read failed.\n"); return 0; } - - if (ray_hdr.dbptr_radial_const != 0) { - begin_data_block = ray_hdr.dbptr_radial_const - ray_hdr_len; - memcpy(&wsr88d_ray->unamb_rng, &radial[begin_data_block+6], 2); - memcpy(&wsr88d_ray->nyq_vel, &radial[begin_data_block+16], 2); - if (little_endian()) { - swap_2_bytes(&wsr88d_ray->unamb_rng); - swap_2_bytes(&wsr88d_ray->nyq_vel); - } - } - else { - wsr88d_ray->unamb_rng = 0; - wsr88d_ray->nyq_vel = 0; - } - - if (ray_hdr.dbptr_ref != 0) { - begin_data_block = ray_hdr.dbptr_ref - ray_hdr_len; - wsr88d_ray->ref = read_data_moment(radial, begin_data_block); - } - if (ray_hdr.dbptr_vel != 0) { - begin_data_block = ray_hdr.dbptr_vel - ray_hdr_len; - wsr88d_ray->vel = read_data_moment(radial, begin_data_block); - } - if (ray_hdr.dbptr_sw != 0) { - begin_data_block = ray_hdr.dbptr_sw - ray_hdr_len; - wsr88d_ray->sw = read_data_moment(radial, begin_data_block); - } - - - /* For testing: print reflectivity data. */ - int prtdata = 0; - { int i; - if (prtdata) { - for (i=0; i < wsr88d_ray->ref->data_hdr.ngates; i++) { - if (i % 10 == 0) printf("\n"); - printf(" %d", wsr88d_ray->ref->data[i]); - } - printf("\n"); - } - } - - free(radial); - return 1; -} - - -void wsr88d_load_ray_data(Data_moment *data_block, Ray *ray) -{ - Data_moment_hdr data_hdr; - int ngates; - int i; - float value, scale, offset; - unsigned char *data; - Range (*invf)(float x); - float (*f)(Range x); - - data_hdr = data_block->data_hdr; - data = data_block->data; - - /* - printf("wsr88d_load_ray_data: "); - printf(" DataName: %s\n", data_hdr.dataname); - */ - - ngates = data_hdr.ngates; - /* - printf(" ngates = %d\n", ngates); - printf(" scale = %f\n", data_hdr.scale); - printf(" offset = %f\n", data_hdr.offset); - */ - offset = data_hdr.offset; - scale = data_hdr.scale; - - /* Note: data range is 2-255. 0 means signal is below threshold, and 1 - * means range folded. + /* We retrieve unambiguous range and Nyquist velocity here so that we don't + * have to do it repeatedly for each data moment later. */ + get_wsr88d_unamb_and_nyq_vel(wsr88d_ray, &unamb_rng, &nyq_vel); + wsr88d_ray->unamb_rng = unamb_rng; + wsr88d_ray->nyq_vel = nyq_vel; - /* Test print - for (i=0; i < 10; i++) { - value = (data[i] - offset) / scale; - printf(" bytevalue = %d, value = %f\n", data[i], value); - } - */ - - /* Reflectivity */ - if (strncmp(data_hdr.dataname, "DREF", 4) == 0) { - /* convert data to float - * convert float to range and put in ray.range - */ - f = DZ_F; - invf = DZ_INVF; - for (i = 0; i < ngates; i++) { - if (data[i] > 1) - value = (data[i] - offset) / scale; - else value = (data[i] == 0) ? BADVAL : RFVAL; - ray->range[i] = invf(value); - ray->h.f = f; - ray->h.invf = invf; - } - } - - /* Velocity */ - if (strncmp(data_hdr.dataname, "DVEL", 4) == 0) { - /* convert data to float - * convert float to range and put in ray.range - */ - f = VR_F; - invf = VR_INVF; - for (i = 0; i < ngates; i++) { - if (data[i] > 1) - value = (data[i] - offset) / scale; - else value = (data[i] == 0) ? BADVAL : RFVAL; - ray->range[i] = invf(value); - ray->h.f = f; - ray->h.invf = invf; - } - } - - /* Spectrum Width */ - if (strncmp(data_hdr.dataname, "DSW", 3) == 0) { - /* convert data to float - * convert float to range and put in ray.range - */ - f = SW_F; - invf = SW_INVF; - for (i = 0; i < ngates; i++) { - if (data[i] > 1) - value = (data[i] - offset) / scale; - else value = (data[i] == 0) ? BADVAL : RFVAL; - ray->range[i] = invf(value); - ray->h.f = f; - ray->h.invf = invf; - } - } - ray->h.range_bin1 = data_hdr.range_first_gate; - ray->h.gate_size = data_hdr.range_samp_interval; - ray->h.nbins = ngates; + return 1; } - void wsr88d_load_ray_hdr(Wsr88d_ray_m31 wsr88d_ray, Ray *ray) { int month, day, year, hour, minute, sec; @@ -501,118 +338,100 @@ void wsr88d_load_ray_hdr(Wsr88d_ray_m31 wsr88d_ray, Ray *ray) ray->h.ray_num = ray_hdr.azm_num; ray->h.elev = ray_hdr.elev; ray->h.elev_num = ray_hdr.elev_num; - ray->h.unam_rng = wsr88d_ray.unamb_rng / 10.; - ray->h.nyq_vel = wsr88d_ray.nyq_vel / 100.; + ray->h.unam_rng = wsr88d_ray.unamb_rng; + ray->h.nyq_vel = wsr88d_ray.nyq_vel; int elev_index; elev_index = ray_hdr.elev_num - 1; ray->h.azim_rate = vcp_data.azim_rate[elev_index]; ray->h.fix_angle = vcp_data.fixed_angle[elev_index]; ray->h.vel_res = vcp_data.vel_res; + if (ray_hdr.azm_res != 1) + ray->h.beam_width = 1.0; + else ray->h.beam_width = 0.5; - /* Get some values using message type 1 routines. - * First load VCP and elevation numbers into msg 1 ray. + /* For convenience, use message type 1 routines to get some values. + * First load VCP and elevation numbers into a msg 1 ray. */ m1_ray.vol_cpat = vcp_data.vcp; m1_ray.elev_num = ray_hdr.elev_num; - m1_ray.unam_rng = wsr88d_ray.unamb_rng; - if (ray_hdr.azm_res != 1) - ray->h.beam_width = 1.0; - else ray->h.beam_width = 0.5; + m1_ray.unam_rng = (short) (wsr88d_ray.unamb_rng * 10.); + /* Get values from message type 1 routines. */ ray->h.frequency = wsr88d_get_frequency(&m1_ray); ray->h.pulse_width = wsr88d_get_pulse_width(&m1_ray); ray->h.pulse_count = wsr88d_get_pulse_count(&m1_ray); - ray->h.prf = wsr88d_get_prf(&m1_ray); + ray->h.prf = (int) wsr88d_get_prf(&m1_ray); ray->h.wavelength = 0.1071; } int wsr88d_get_vol_index(char* dataname) { - int vol_index; + int vol_index = -1; if (strncmp(dataname, "DREF", 4) == 0) vol_index = DZ_INDEX; if (strncmp(dataname, "DVEL", 4) == 0) vol_index = VR_INDEX; - if (strncmp(dataname, "DSW", 3) == 0) vol_index = SW_INDEX; - /* TODO: Add the other data moments. */ + if (strncmp(dataname, "DSW", 3) == 0) vol_index = SW_INDEX; + if (strncmp(dataname, "DZDR", 3) == 0) vol_index = DR_INDEX; + if (strncmp(dataname, "DPHI", 3) == 0) vol_index = PH_INDEX; + if (strncmp(dataname, "DRHO", 3) == 0) vol_index = RH_INDEX; return vol_index; } -void wsr88d_load_ray_into_radar(Wsr88d_ray_m31 wsr88d_ray, int isweep, int iray, - Radar *radar) +#define MAXRAYS_M31 800 +#define MAXSWEEPS 20 + +void wsr88d_load_ray(Wsr88d_ray_m31 wsr88d_ray, int data_ptr, + int isweep, int iray, Radar *radar) { - Ray *ray; + /* Load data into ray structure for this field or data moment. */ + + Data_moment_hdr data_hdr; + int ngates; + int i, hdr_size; + float value, scale, offset; + unsigned char *data; Range (*invf)(float x); float (*f)(Range x); + Ray *ray; int vol_index, waveform; - enum waveforms {surveillance=1, doppler_ambres, doppler_no_ambres, batch}; + enum waveforms {surveillance=1, doppler_w_amb_res, doppler_no_amb_res, + batch}; -/* for each data moment do - get the name of the data moment - get new volume if needed - get new sweep if needed - get new ray - put this data moment in ray - endfor each data moment -*/ - - /* Note: The data moment type can only be determined by the name within the - * individual data moment block. The name of the pointer is not reliable. - * For example, in legacy resolution mode, dbptr_ref points to the velocity - * data moment in velocity split cuts. Apparently the location of the first - * data moment is always stored in the reflectivity pointer (dbptr_ref), and - * sometimes this is velocity. When this occurs, the velocity data pointer - * (dbptr_vel) then points to spectrum width. - * With super resolution, there actually is reflectivity in the velocity - * split cut. It's there for use with the recombination algorithm. - */ + /* Get data moment header. */ + hdr_size = sizeof(data_hdr); + memcpy(&data_hdr, &wsr88d_ray.data[data_ptr], hdr_size); + data_ptr += hdr_size; + if (little_endian()) wsr88d_swap_data_hdr(&data_hdr); - if (wsr88d_ray.ray_hdr.dbptr_ref > 0) { - vol_index = wsr88d_get_vol_index(wsr88d_ray.ref->data_hdr.dataname); - switch (vol_index) { - case DZ_INDEX: f = DZ_F; invf = DZ_INVF; break; - case VR_INDEX: f = VR_F; invf = VR_INVF; break; - case SW_INDEX: f = SW_F; invf = SW_INVF; break; - default: f = DZ_F; invf = DZ_INVF; break; - } - /* If this is reflectivity, check the waveform type to make sure - * it isn't from a Doppler split cut. - * We only keep reflectivity if the waveform type is surveillance or - * batch, or the elevation is above the split cut elevations. - */ - waveform = vcp_data.waveform[isweep]; - if (vol_index != DZ_INDEX || - (waveform == surveillance || waveform == batch || - vcp_data.fixed_angle[isweep] >= 6.0)) - { - if (radar->v[vol_index] == NULL) { - radar->v[vol_index] = RSL_new_volume(MAXSWEEPS); - radar->v[vol_index]->h.f = f; - radar->v[vol_index]->h.invf = invf; - } - if (radar->v[vol_index]->sweep[isweep] == NULL) { - radar->v[vol_index]->sweep[isweep] = RSL_new_sweep(MAXRAYS_M31); - radar->v[vol_index]->sweep[isweep]->h.f = f; - radar->v[vol_index]->sweep[isweep]->h.invf = invf; - } - ray = RSL_new_ray(wsr88d_ray.ref->data_hdr.ngates); - wsr88d_load_ray_data(wsr88d_ray.ref, ray); - wsr88d_load_ray_hdr(wsr88d_ray, ray); - radar->v[vol_index]->sweep[isweep]->ray[iray] = ray; - radar->v[vol_index]->sweep[isweep]->h.nrays = iray+1; - } + vol_index = wsr88d_get_vol_index(data_hdr.dataname); + if (vol_index < 0) { + fprintf(stderr,"wsr88d_load_ray: Unknown dataname %s. isweep = %d, " + "iray = %d.\n", data_hdr.dataname, isweep, iray); + return; + } + switch (vol_index) { + case DZ_INDEX: f = DZ_F; invf = DZ_INVF; break; + case VR_INDEX: f = VR_F; invf = VR_INVF; break; + case SW_INDEX: f = SW_F; invf = SW_INVF; break; + case DR_INDEX: f = DR_F; invf = DR_INVF; break; + case PH_INDEX: f = PH_F; invf = PH_INVF; break; + case RH_INDEX: f = RH_F; invf = RH_INVF; break; } - if (wsr88d_ray.ray_hdr.dbptr_vel > 0) { - vol_index = wsr88d_get_vol_index(wsr88d_ray.vel->data_hdr.dataname); - switch (vol_index) { - case DZ_INDEX: f = DZ_F; invf = DZ_INVF; break; - case VR_INDEX: f = VR_F; invf = VR_INVF; break; - case SW_INDEX: f = SW_F; invf = SW_INVF; break; - default: f = DZ_F; invf = DZ_INVF; break; - } + waveform = vcp_data.waveform[isweep]; + /* The following somewhat complicated if-statement says we'll do the + * body of the statement if: + * a) the data moment is not reflectivity, or + * b) the data moment is reflectivity and one of the following is true: + * - waveform is surveillance + * - waveform is batch + * - elevation is greater than highest split cut (i.e., 6 deg.) + */ + if (vol_index != DZ_INDEX || (waveform == surveillance || + waveform == batch || vcp_data.fixed_angle[isweep] >= 6.0)) { if (radar->v[vol_index] == NULL) { radar->v[vol_index] = RSL_new_volume(MAXSWEEPS); radar->v[vol_index]->h.f = f; @@ -623,61 +442,69 @@ void wsr88d_load_ray_into_radar(Wsr88d_ray_m31 wsr88d_ray, int isweep, int iray, radar->v[vol_index]->sweep[isweep]->h.f = f; radar->v[vol_index]->sweep[isweep]->h.invf = invf; } - ray = RSL_new_ray(wsr88d_ray.vel->data_hdr.ngates); - wsr88d_load_ray_data(wsr88d_ray.vel, ray); - wsr88d_load_ray_hdr(wsr88d_ray, ray); - radar->v[vol_index]->sweep[isweep]->ray[iray] = ray; - radar->v[vol_index]->sweep[isweep]->h.nrays = iray+1; - } + ngates = data_hdr.ngates; + ray = RSL_new_ray(ngates); - if (wsr88d_ray.ray_hdr.dbptr_sw > 0) { - vol_index = wsr88d_get_vol_index(wsr88d_ray.sw->data_hdr.dataname); - switch (vol_index) { - case DZ_INDEX: f = DZ_F; invf = DZ_INVF; break; - case VR_INDEX: f = VR_F; invf = VR_INVF; break; - case SW_INDEX: f = SW_F; invf = SW_INVF; break; - default: f = DZ_F; invf = DZ_INVF; break; - } - if (radar->v[vol_index] == NULL) { - radar->v[vol_index] = RSL_new_volume(MAXSWEEPS); - radar->v[vol_index]->h.f = f; - radar->v[vol_index]->h.invf = invf; - } - if (radar->v[vol_index]->sweep[isweep] == NULL) { - radar->v[vol_index]->sweep[isweep] = RSL_new_sweep(MAXRAYS_M31); - radar->v[vol_index]->sweep[isweep]->h.f = f; - radar->v[vol_index]->sweep[isweep]->h.invf = invf; + /* Convert data to float, then use range function to store in ray. + * Note: data range is 2-255. 0 means signal is below threshold, and 1 + * means range folded. + */ + + offset = data_hdr.offset; + scale = data_hdr.scale; + if (data_hdr.scale == 0) scale = 1.0; + data = &wsr88d_ray.data[data_ptr]; + for (i = 0; i < ngates; i++) { + if (data[i] > 1) + value = (data[i] - offset) / scale; + else value = (data[i] == 0) ? BADVAL : RFVAL; + ray->range[i] = invf(value); + ray->h.f = f; + ray->h.invf = invf; } - ray = RSL_new_ray(wsr88d_ray.sw->data_hdr.ngates); - wsr88d_load_ray_data(wsr88d_ray.sw, ray); wsr88d_load_ray_hdr(wsr88d_ray, ray); + ray->h.range_bin1 = data_hdr.range_first_gate; + ray->h.gate_size = data_hdr.range_samp_interval; + ray->h.nbins = ngates; radar->v[vol_index]->sweep[isweep]->ray[iray] = ray; radar->v[vol_index]->sweep[isweep]->h.nrays = iray+1; } - } +void wsr88d_load_ray_into_radar(Wsr88d_ray_m31 wsr88d_ray, int isweep, int iray, + Radar *radar) +{ + int *data_ptr, hdr_size; + int i, ndatablocks, nconstblocks = 3; + + hdr_size = sizeof(wsr88d_ray.ray_hdr); + + ndatablocks = wsr88d_ray.ray_hdr.data_block_count; + data_ptr = (int *) &wsr88d_ray.ray_hdr.dbptr_ref; + for (i=0; i < ndatablocks-nconstblocks; i++) { + wsr88d_load_ray(wsr88d_ray, *data_ptr-hdr_size, isweep, iray, radar); + data_ptr++; + } +} + -void wsr88d_load_sweep_header(Radar *radar, int isweep, - Wsr88d_ray_m31 wsr88d_ray) +void wsr88d_load_sweep_header(Radar *radar, int isweep) { - int ivolume; - Ray_header_m31 ray_hdr; + int ivolume, nrays; Sweep *sweep; - int vcp; - - ray_hdr = wsr88d_ray.ray_hdr; + Ray *last_ray; for (ivolume=0; ivolume < MAX_RADAR_VOLUMES; ivolume++) { - if (radar->v[ivolume] != NULL && radar->v[ivolume]->sweep[isweep] != NULL) { + if (radar->v[ivolume] != NULL && + radar->v[ivolume]->sweep[isweep] != NULL) { sweep = radar->v[ivolume]->sweep[isweep]; - radar->v[ivolume]->sweep[isweep]->h.sweep_num = ray_hdr.elev_num; - radar->v[ivolume]->sweep[isweep]->h.elev = - vcp_data.fixed_angle[isweep]; - if (ray_hdr.azm_res != 1) - sweep->h.beam_width = 1.0; - else sweep->h.beam_width = 0.5; + nrays = sweep->h.nrays; + if (nrays == 0) continue; + last_ray = sweep->ray[nrays-1]; + sweep->h.sweep_num = last_ray->h.elev_num; + sweep->h.elev = vcp_data.fixed_angle[isweep]; + sweep->h.beam_width = last_ray->h.beam_width; sweep->h.vert_half_bw = sweep->h.beam_width / 2.; sweep->h.horz_half_bw = sweep->h.beam_width / 2.; } @@ -694,41 +521,64 @@ Radar *load_wsr88d_m31_into_radar(Wsr88d_file *wf) int msg_hdr_size, msg_size, n; Radar *radar = NULL; -/* Message type 31 is a variable length message, whereas all other message - * types are made up of 2432 byte segments. To handle this, we read the - * message header and check the message type. If it is not 31, then we read - * the remainder of the constant size segment. If message type is 31, we read - * the remainder of the message by the size given in message header. - * When reading the message header, we must include 12 bytes inserted - * by RPG, which we ignore, followed by the 8 halfwords (short integers) which - * make up the actual message header. - * For more information, see the "Interface Control Document for the RDA/RPG" - * at the WSR-88D Radar Operations Center web site. - */ + /* Message type 31 is a variable length message. All other message types + * are made up of 2432-byte segments. To handle all types, we read the + * message header and check the message type. If it is not 31, we simply + * read the remainder of the 2432-byte segment. If message type is 31, we + * use the size given in message header to determine how many bytes to + * read. For more information, see the "Interface Control Document for the + * RDA/RPG" at the WSR-88D Radar Operations Center web site. + */ n = fread(&msghdr, sizeof(Wsr88d_msg_hdr), 1, wf->fptr); /* printf("msgtype = %d\n", msghdr.msg_type); */ msg_hdr_size = sizeof(Wsr88d_msg_hdr) - sizeof(msghdr.rpg); - radar = RSL_new_radar(MAX_RADAR_VOLUMES); - + while (! end_of_vos) { if (msghdr.msg_type == 31) { if (little_endian()) wsr88d_swap_m31_hdr(&msghdr); - /* Get size in bytes of message following message header. - * The size given in message header is in halfwords, so double it. + /* Get size of the remainder of message. The given size is in + * halfwords, but we want it in bytes, so double it. */ msg_size = (int) msghdr.msg_size * 2 - msg_hdr_size; n = read_wsr88d_ray_m31(wf, msg_size, &wsr88d_ray); + /* Assume error message was issued from read routine */ if (n <= 0) return NULL; - /* Load this ray into radar structure ray. */ + /* We need to check radial status for start of new elevation. + * Sometimes this occurs without an end-of-elevation flag for the + * previous sweep, which is the trigger for loading the sweep + * header. "iray" is set to zero after end-of-elev is received, + * so that's why we check it. We issue a warning because when this + * occurs, the number of rays in the previous sweep is usually less + * than expected. + */ + if (wsr88d_ray.ray_hdr.radial_status == START_OF_ELEV && + iray != 0) { + fprintf(stderr,"Warning: Radial status is Start-of-Elevation, " + "but End-of-Elevation was not\n" + "issued for elevation number %d. Number of rays = %d" + "\n", isweep+1, iray); + wsr88d_load_sweep_header(radar, isweep); + isweep++; + iray = 0; + } + + /* Load ray into radar structure. */ wsr88d_load_ray_into_radar(wsr88d_ray, isweep, iray, radar); iray++; + if (iray >= MAXRAYS_M31) { + fprintf(stderr,"Error: iray = %d, equals or exceeds MAXRAYS_M31" + " (%d)\n", iray, MAXRAYS_M31); + fprintf(stderr,"isweep = %d\n", isweep); + RSL_free_radar(radar); + return NULL; + } } else { /* msg_type not 31 */ n = fread(&non31_seg_remainder, sizeof(non31_seg_remainder), 1, @@ -737,10 +587,11 @@ Radar *load_wsr88d_m31_into_radar(Wsr88d_file *wf) fprintf(stderr,"Warning: load_wsr88d_m31_into_radar: "); if (feof(wf->fptr) != 0) fprintf(stderr, "Unexpected end of file.\n"); - else fprintf(stderr,"Read failed.\n"); - fprintf(stderr,"Current sweep number: %d\n" - "Last ray read: %d\n", isweep+1, iray); - wsr88d_load_sweep_header(radar, isweep, wsr88d_ray); + else + fprintf(stderr,"Read failed.\n"); + fprintf(stderr,"Current sweep index: %d\n" + "Last ray read: %d\n", isweep, iray); + wsr88d_load_sweep_header(radar, isweep); return radar; } if (msghdr.msg_type == 5) { @@ -752,11 +603,12 @@ Radar *load_wsr88d_m31_into_radar(Wsr88d_file *wf) /* Check for end of sweep */ if (wsr88d_ray.ray_hdr.radial_status == END_OF_ELEV) { - wsr88d_load_sweep_header(radar, isweep, wsr88d_ray); + wsr88d_load_sweep_header(radar, isweep); isweep++; iray = 0; } + /* If not at end of volume scan, read next message header. */ if (wsr88d_ray.ray_hdr.radial_status != END_VOS) { n = fread(&msghdr, sizeof(Wsr88d_msg_hdr), 1, wf->fptr); if (n < 1) { @@ -764,18 +616,18 @@ Radar *load_wsr88d_m31_into_radar(Wsr88d_file *wf) if (feof(wf->fptr) != 0) fprintf(stderr, "Unexpected end of file.\n"); else fprintf(stderr,"Failed reading msghdr.\n"); - fprintf(stderr,"Current sweep number: %d\n" - "Last ray read: %d\n", isweep+1, iray); - wsr88d_load_sweep_header(radar, isweep, wsr88d_ray); + fprintf(stderr,"Current sweep index: %d\n" + "Last ray read: %d\n", isweep, iray); + wsr88d_load_sweep_header(radar, isweep); return radar; } } else { end_of_vos = 1; - wsr88d_load_sweep_header(radar, isweep, wsr88d_ray); + wsr88d_load_sweep_header(radar, isweep); } if (feof(wf->fptr) != 0) end_of_vos = 1; - } + } /* while not end of vos */ return radar; } diff --git a/wsr88d_to_radar.c b/wsr88d_to_radar.c index 9ef0597..0ea0672 100644 --- a/wsr88d_to_radar.c +++ b/wsr88d_to_radar.c @@ -223,6 +223,7 @@ Radar *RSL_wsr88d_to_radar(char *infile, char *call_or_first_tape_file) Wsr88d_tape_header wsr88d_tape_header; int n; int nsweep; + int i; int iv; int nvolumes; int volume_mask[] = {WSR88D_DZ, WSR88D_VR, WSR88D_SW}; @@ -300,14 +301,15 @@ Radar *RSL_wsr88d_to_radar(char *infile, char *call_or_first_tape_file) */ if (n > 0) { strncpy(version, wsr88d_file_header.title.filename, 8); - if (strncmp(version,"AR2V0004",8) == 0 || - strncmp(version,"AR2V0003",8) ==0 || - strncmp(version,"AR2V0002",8) == 0) { - expected_msgtype = 31; + if (strncmp(version,"AR2V0006",8) == 0 || + strncmp(version,"AR2V0004",8) == 0 || + strncmp(version,"AR2V0003",8) == 0 || + strncmp(version,"AR2V0002",8) == 0) { + expected_msgtype = 31; } else if (strncmp(version,"ARCHIVE2",8) == 0 || - strncmp(version,"AR2V0001",8) == 0) { - expected_msgtype = 1; + strncmp(version,"AR2V0001",8) == 0) { + expected_msgtype = 1; } } @@ -386,6 +388,14 @@ Radar *RSL_wsr88d_to_radar(char *infile, char *call_or_first_tape_file) } } } + if (nsweep == 0) { + /* Get Volume Coverage Pattern number for radar header. */ + i=0; + while (i < MAX_RAYS_IN_SWEEP && wsr88d_sweep.ray[i] == NULL) i++; + if (i < MAX_RAYS_IN_SWEEP) radar->h.vcp = wsr88d_get_volume_coverage( + wsr88d_sweep.ray[i]); + } + free_and_clear_sweep(&wsr88d_sweep, 0, MAX_RAYS_IN_SWEEP); } -- 2.43.2