From: Andy Spencer Date: Sun, 27 Apr 2014 23:05:57 +0000 (+0000) Subject: RSL v1.44 X-Git-Url: http://pileus.org/git/?p=~andy%2Frsl;a=commitdiff_plain;h=01560cc947c59aa600473158dfe26ea0dfaf0823 RSL v1.44 --- diff --git a/CHANGES b/CHANGES index 1d956c1..d9af1cf 100644 --- a/CHANGES +++ b/CHANGES @@ -1,7 +1,57 @@ /* Changes for RSL * *--------------------------------------------------------------------- - * v1.43 Released 4/30/2011 + * v1.44 Released 11/21/2013 + * + * 1. wsr88d_m31.c: Changed the criteria for discarding reflectivity belonging + * to a split-cut velocity sweep. Previously checked for elev < 6 degrees + * to identify split cuts; now check waveform for "Contiguous Doppler". The + * elev test was inadequate for VCP 31, which has non-split cuts below 6 + * degrees. Thanks to Chris Theisen for the bug report. + * Also added code to handle field selection. + * 2. New routines: + * RSL_wsr88d_keep_short_refl: + * Calling this function instructs the WSR-88D Level II ingest to keep + * short-range reflectivity collected in the velocity portion of a split + * cut. Normally this short-range reflectivity is discarded. + * Calling this function also changes how split cuts are handled in the + * Radar structure and when writing to UF. Split cut surveillance and + * Dopper sweeps are normally placed in the same sweep index in the Radar + * structure, as though they were taken in a single sweep, so that + * reflectivity (long-range), radial velocity, and spectrum width are at + * the same sweep index. If RSL_wsr88d_keep_short_refl() is called, each + * individual sweep of a split cut is stored at a separate index, so for + * example, the first sweep contains only long-range reflectivity (and + * dual-pol fields if present), and the second sweep contains short-range + * reflectivity, radial velocity, and spectrum width. This pattern is + * repeated with the remaining split cuts. If a UF is produced, the + * split cuts are written in like manner. + * wsr88d_align_splitcut_rays: + * This is called by RSL_radar_to_uf when writing WSR-88D data to UF. It + * reorders rays in the velocity sweep of a split-cut to match those in + * the reflectivity sweep by azimuth. + * wsr88d_merge_split_cuts: + * Moves the surveillance (long-range reflectivity) and Doppler (radial + * velocity) sweeps in a WSR-88D split cut to a single sweep index in the + * Radar structure. It also takes care of the multiple velocity sweeps in + * each VCP 121 split cut by assigning them unique field names. + * 3. radar_to_uf.c: Added code to handle WSR-88D split cuts. Also Changed + * optional header to use time from radar header instead of first ray. + * This was done because some radar data is apparently sorted by azimuth, + * which results in the first ray not being the chronological first ray. + * 4. gzip.c (uncompress_pipe): Implemented the fix for the "too many open + * files" problem, sent independently by Yu Zhang and Lee Burchett. + * 5. Restored EDGE routines for backward compatibility. + * 6. nsig_to_radar.c: Fixed bug where string exceeded size of variable into + * which it was written (radar->h.state). This caused a buffer overflow + * error on Ubuntu systems. Thanks to Jonathan Helmus for finding this one. + * 7. volume.c (DM_F, DM_INVF): Thanks to David Kingsmill for providing improved + * conversion functions for Returned Power (DM). + * volume.c (RH_F, RH_INVF, PH_F, PH_INVF): These conversion functions were + * specific to Sigmet raw data. Made them more general. + * 8. Thanks to Mark Parncutt for updated wsr88d_locations.dat. + *--------------------------------------------------------------------- + * v1.43 Released 4/30/2012 * * 1. nsig_to_radar.c: Added antenna scan mode to radar header. * Added azimuth to sweep header for RHI. diff --git a/Makefile.am b/Makefile.am index b34824a..895521d 100644 --- a/Makefile.am +++ b/Makefile.am @@ -8,17 +8,20 @@ colordir = $(libdir)/colors lib_LTLIBRARIES = librsl.la -librsl_la_LDFLAGS = -version-info 1:43 +librsl_la_LDFLAGS = -version-info 1:44 librsl_la_SOURCES = \ $(rapic_c) $(radtec_c)\ dorade.c dorade_print.c dorade_to_radar.c\ lassen.c lassen_to_radar.c \ +edge_to_radar.c \ radar.c volume.c image_gen.c cappi.c fraction.c read_write.c farea.c \ range.c radar_to_uf.c uf_to_radar.c wsr88d_to_radar.c \ carpi.c cube.c sort_rays.c toga_to_radar.c gts.c histogram.c \ ray_indexes.c anyformat_to_radar.c get_win.c endian.c mcgill_to_radar.c \ mcgill.c interp.c toga.c wsr88d.c wsr88d_get_site.c wsr88d_m31.c \ gzip.c prune.c reverse.c fix_headers.c \ + wsr88d_align_split_cut_rays.c \ + wsr88d_merge_split_cuts.c \ nsig_to_radar.c nsig.c nsig2_to_radar.c \ africa_to_radar.c africa.c \ radar_to_hdf_2.c hdf_to_radar.c toolkit_memory_mgt.c \ diff --git a/README b/README index 4ab945a..b43a860 100644 --- a/README +++ b/README @@ -1,4 +1,4 @@ -v1.43 (Released 4/30/2012) +v1.44 (Released 11/21/2013) This is the README file for the Radar Software Library (RSL). diff --git a/THANKS b/THANKS index fe07807..bbde807 100644 --- a/THANKS +++ b/THANKS @@ -1,15 +1,20 @@ Thanks to the following folks who helped make RSL a better product. -Ther contributions include reporting bugs, offering suggestions, and sending -code for bug fixes. Thanks again! +Their contributions include bug reports, suggestions, and/or code. +Mark Askelson Cesar Beneti Thiago Biscaro Stacy Brodzik Eric Bruning +Lee Burchett Scott Collis Patrick Gatlin John Hall +Jonathan Helmus David Kingsmill +Mark Parncutt Fabio Sato Daniel Sheldon Andy Spencer +Chris Theisen +Yu Zhang diff --git a/anyformat_to_radar.c b/anyformat_to_radar.c index 6f0d8fb..45e3667 100644 --- a/anyformat_to_radar.c +++ b/anyformat_to_radar.c @@ -157,7 +157,6 @@ Radar *RSL_anyformat_to_radar(char *infile, ...) case TOGA_FILE: radar = RSL_toga_to_radar(infile); break; case NSIG_FILE_V1: radar = RSL_nsig_to_radar(infile); break; case NSIG_FILE_V2: radar = RSL_nsig2_to_radar(infile); break; - case RAPIC_FILE: radar = RSL_rapic_to_radar(infile); break; case RADTEC_FILE: radar = RSL_radtec_to_radar(infile); break; case RSL_FILE: radar = RSL_read_radar(infile); break; #ifdef HAVE_LIBTSDISTK diff --git a/configure.in b/configure.in index 7f7d8b6..8dbb3ac 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.43) +AC_INIT(rsl, v1.44) AC_CONFIG_SRCDIR(volume.c) AM_INIT_AUTOMAKE diff --git a/doc/RSL_anyformat_to_radar.html b/doc/RSL_anyformat_to_radar.html index 53074b4..e3adb1b 100644 --- a/doc/RSL_anyformat_to_radar.html +++ b/doc/RSL_anyformat_to_radar.html @@ -38,7 +38,6 @@ thus, are recognized by RSL_anyformat_to_radar:

NSIG (SIGMET version 1 and version 2 raw product files)

MCGILL

RSL output format -

RAPIC (Berrimah Austrailia)

RADTEC (SPANDAR, requires PKWARE compression tools)

DORADE

For TOGA and KWAJALEIN files, there is no standard magic information @@ -76,8 +75,6 @@ See also RSL_nsig2_to_radar, RSL_mcgill_to_radar RSL_kwaj_to_radar, -RSL_rapic_to_radar -, RSL_radtec_to_radar.

diff --git a/doc/RSL_radar_intro.html b/doc/RSL_radar_intro.html index 08f49b3..e30d371 100644 --- a/doc/RSL_radar_intro.html +++ b/doc/RSL_radar_intro.html @@ -19,13 +19,13 @@ Introduction

By John H. Merritt and David B. Wolff; NASA/TRMM Office
-Software Verson 1.42 (released 7/12/2011)

+Software Verson 1.44 (released 11/21/2013)
This library is an object oriented programming environment to keep application programming simple, for the casual C programmer, as well as for analysis software applicable to all RADAR data related to the TRMM GV effort. This library reads the wsr88d, lassen, kwajalein, mcgill, toga, -UF, RAPIC and native RSL file formats. The most important functions provided +UF, and native RSL file formats. The most important functions provided are those which load any one of the radar file formats into memory; see RSL_anyformat_to_radar. Additional functions are provide to mainpulate the RSL objects. Nearly all of the diff --git a/doc/RSL_ray_header_struct.html b/doc/RSL_ray_header_struct.html index 500d6e5..b9b35fb 100644 --- a/doc/RSL_ray_header_struct.html +++ b/doc/RSL_ray_header_struct.html @@ -13,7 +13,7 @@   
-
typedef struct {
  int   month; /* Date for this ray; month (1-12). */
  int   day;   /* Date for this ray; day (1-31).   */
  int   year;  /* Date for this ray; year (eg. 1993). */
  int   hour;  /* Time for this ray; hour (0-23). */
  int   minute;/* Time for this ray; minute (0-59).*/
  float sec;   /* Time 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.
                    * This angle is the mean azimuth for the whole ray.
                    * Eg. 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. */
  
  int   range_bin1; /* Range to first gate.(meters) */
  int   gate_size;  /* Data gate size (meters)*/
  
  float  vel_res;    /* Doppler velocity resolution */
  float sweep_rate;   /* Sweep rate. Full sweeps/min. */
  
  int prf;          /* Pulse repitition frequency, in Hz. */
  int prf2;         /* Second PRF, for Sigmet dual PRF. */
  float azim_rate;  /* Sweep rate in degrees/second.*/
  float fix_angle;  /* Elevation angle for the sweep. (degrees). */
  float pitch;      /* Pitch angle. */
  float roll;       /* Roll  angle. */
  float heading;    /* Heading. */
  float pitch_rate; /* (angle/sec) */
  float roll_rate;  /* (angle/sec) */
  float heading_rate; /* (angle/sec) */
  float lat;          /* Latitude (degrees) */
  float lon;          /* Longitude (degrees) */
  int   alt;          /* Altitude (m) */
  float rvc;          /* Radial velocity correction (m/sec) */
  float vel_east;     /* Platform velocity to the east (negative for west)   (m/sec) */
  float vel_north;    /* Platform velocity to the north (negative for south) (m/sec) */
  float vel_up;       /* Platform velocity toward up (negative for down)     (m/sec) */
  float pulse_count;  /* Pulses used in a single dwell time. */
  float pulse_width; /* Pulse width (micro-sec). */
  float beam_width;  /* Beamwidth in degrees. */
  float frequency;   /* Bandwidth MHz. */
  float wavelength;  /* Wavelength. Meters. */
  float nyq_vel;    /* Nyquist velocity. m/s */
  float (*f)(Range x);       /* Data conversion function. f(x). */
  Range (*invf)(float x);    /* Data conversion function. invf(x). */
  int   nbins;               /* Number of array elements for 'Range'. */
} Ray_header;
+
typedef struct {
  int   month; /* Date for this ray; month (1-12). */
  int   day;   /* Date for this ray; day (1-31).   */
  int   year;  /* Date for this ray; year (eg. 1993). */
  int   hour;  /* Time for this ray; hour (0-23). */
  int   minute;/* Time for this ray; minute (0-59).*/
  float sec;   /* Time 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.
                    * This angle is the mean azimuth for the whole ray.
                    * Eg. 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. */
  
  int   range_bin1; /* Range to first gate.(meters) */
  int   gate_size;  /* Data gate size (meters)*/
  
  float  vel_res;    /* Doppler velocity resolution */
  float sweep_rate;   /* Sweep rate. Full sweeps/min. */
  
  int prf;          /* Pulse repitition frequency, in Hz. */
  int prf2;         /* Second PRF, for Sigmet dual PRF. */
  float azim_rate;  /* Sweep rate in degrees/second.*/
  float fix_angle;  /* Elevation angle for the sweep. (degrees). */
  float pitch;      /* Pitch angle. */
  float roll;       /* Roll  angle. */
  float heading;    /* Heading. */
  float pitch_rate; /* (angle/sec) */
  float roll_rate;  /* (angle/sec) */
  float heading_rate; /* (angle/sec) */
  float lat;          /* Latitude (degrees) */
  float lon;          /* Longitude (degrees) */
  int   alt;          /* Altitude (m) */
  float rvc;          /* Radial velocity correction (m/sec) */
  float vel_east;     /* Platform velocity to the east (negative for west)   (m/sec) */
  float vel_north;    /* Platform velocity to the north (negative for south) (m/sec) */
  float vel_up;       /* Platform velocity toward up (negative for down)     (m/sec) */
  float pulse_count;  /* Pulses used in a single dwell time. */
  float pulse_width; /* Pulse width (micro-sec). */
  float beam_width;  /* Beamwidth in degrees. */
  float frequency;   /* Bandwidth GHz. */
  float wavelength;  /* Wavelength. Meters. */
  float nyq_vel;    /* Nyquist velocity. m/s */
  float (*f)(Range x);       /* Data conversion function. f(x). */
  Range (*invf)(float x);    /* Data conversion function. invf(x). */
  int   nbins;               /* Number of array elements for 'Range'. */
} Ray_header;

diff --git a/doc/RSL_wsr88d_keep_short_refl.html b/doc/RSL_wsr88d_keep_short_refl.html new file mode 100644 index 0000000..36925c7 --- /dev/null +++ b/doc/RSL_wsr88d_keep_short_refl.html @@ -0,0 +1,41 @@ + + + + + +
+ + +

RSL_wsr88d_keep_short_refl

+ +
+ +

Synopsis

+#include "rsl.h"
+void RSL_wsr88d_keep_short_refl(void); + +

+
Description

+Calling this function instructs the WSR-88D Level II ingest to keep +short-range reflectivity collected in the velocity portion of a split +cut. Normally this short-range reflectivity is discarded. +Calling this function also changes how split cuts are handled in the +Radar structure and when writing to UF. Split cut surveillance and +Dopper sweeps are normally placed in the same sweep index in the Radar +structure, as though they were taken in a single sweep, so that +reflectivity (long-range), radial velocity, and spectrum width are at +the same sweep index. If RSL_wsr88d_keep_short_refl is called, each +individual sweep of a split cut is stored at a separate index, so for +example, the first sweep contains only long-range reflectivity (and +dual-pol fields if present), and the second sweep contains short-range +reflectivity, radial velocity, and spectrum width. This pattern is +repeated with the remaining split cuts. If a UF is produced, the +split cuts are written in like manner. +
+
+ +

Return value

+None. +
+ + diff --git a/doc/functionality_index.html b/doc/functionality_index.html index 5df2230..bdbee77 100644 --- a/doc/functionality_index.html +++ b/doc/functionality_index.html @@ -23,7 +23,6 @@ Input
Radar *RSL_nsig2_to_radar(char *infile);
Radar *RSL_radtec_to_radar(char *infile); -
Radar *RSL_rapic_to_radar(char *infile);
Radar *RSL_read_radar(char *infile);
Radar *RSL_toga_to_radar(char *infile);
Radar *RSL_uf_to_radar(char *infile); @@ -41,6 +40,7 @@ Input mds, float calibr_slope, float calibr_intercept);
void RSL_radar_verbose_off(void);
void RSL_radar_verbose_on(void); +
void RSL_wsr88d_keep_short_refl(void);

Output

void RSL_radar_to_uf(Radar *r, char *outfile); diff --git a/doc/index.html b/doc/index.html index 5ecf853..33d5f5d 100644 --- a/doc/index.html +++ b/doc/index.html @@ -21,7 +21,7 @@

-

Current RSL Version 1.43, released 4/30/2012

+

Current RSL Version 1.44, released 11/21/2013

In support of the Tropical Rainfall Measuring Mission's (TRMM) Global @@ -31,7 +31,7 @@ library written in C.

This library is an object oriented programming environment for writing software applicable to all RADAR data related to the TRMM GV effort. This library reads the WSR88D, Lassen, Sigmet, McGill, UF, -HDF, RAPIC, RADTEC and native RSL file formats. Additional functions are +HDF, RADTEC and native RSL file formats. Additional functions are provided to manipulate the RSL objects. Nearly all of the functions return objects. When they don't, they usually perform actions like output, making images, etc. The most general object in RSL is Radar. The structure @@ -49,7 +49,7 @@ etc. There are approximately 20 field types. See the Users Guide for more information.  Also, check out What's New.
-Send questions or comments to help@radar.gsfc.nasa.gov. +Send questions or comments to gsfc-rsl-help@lists.nasa.gov.

Indexed Reference
 @@ -82,7 +82,7 @@ Supported Radar Data Formats

-Lassen +Sigmet Raw Product Yes @@ -90,7 +90,7 @@ Supported Radar Data Formats -McGill +WSR-88D Level II Yes @@ -98,15 +98,15 @@ Supported Radar Data Formats -Sigmet +UF - Universal Format Yes -No +Yes -WSR-88D +Lassen Yes @@ -114,7 +114,7 @@ Supported Radar Data Formats -RAPIC +McGill Yes @@ -122,7 +122,7 @@ Supported Radar Data Formats -UF - Universal Format +HDF 1B-51, 1C-51 Yes @@ -130,15 +130,37 @@ Supported Radar Data Formats -HDF - Hierarchical Data Format +RADTEC Yes -Yes +No +
+ +
+ - + + + + + + + + + + + + + + + + diff --git a/doc/quick_ref.html b/doc/quick_ref.html index 314d8da..e63cca4 100644 --- a/doc/quick_ref.html +++ b/doc/quick_ref.html @@ -37,6 +37,7 @@ Alphabetical index to RSL routines (alphabetical and grouped by object returned)Radar *RSL_anyformat_to_radar(char *infile [, char *callid_or_first_file]); +
Radar *RSL_edge_to_radar(char *infile);  (Not supported)
Radar *RSL_get_window_from_radar(Radar *r, float min_range, float max_range, float low_azim, float hi_azim);
Radar *RSL_hdf_to_radar(char *infile); @@ -50,7 +51,7 @@ float min_range, float max_range, float low_azim, float hi_azim);
Radar *RSL_prune_radar(Radar *radar);
Radar *RSL_radtec_to_radar(char *infile); -
Radar *RSL_rapic_to_radar(char *infile); +
Radar *RSL_rapic_to_radar(char *infile);  (Not supported)
Radar *RSL_read_radar(char *infile);
Radar *RSL_sort_radar(Radar *r);
Radar *RSL_toga_to_radar(char *infile); @@ -257,6 +258,7 @@ char *image, int xdim, int ydim, char c_ table[256][3]); char *image, int xdim, int ydim);
void RSL_write_ppm(char *outfile, unsigned char *image, int xdim, int ydim, char c_t able[256][3]); +
void RSL_wsr88d_keep_short_refl(void);
Cappi *RSL_cappi_at_h(Volume *v, float h, float grnd_range);
Cappi *RSL_new_cappi(Sweep *sweep, float diff --git a/doc/users_guide.html b/doc/users_guide.html index fb575c0..9f582a0 100644 --- a/doc/users_guide.html +++ b/doc/users_guide.html @@ -126,6 +126,12 @@ Output routine + + + + + +
+

+Not Supported but Retained for Backward Compatibility

+
RADTECFormat TypeInputOutput
EDGEYesNo
RAPIC Yes None
EDGERSL_edge_to_radarNone
RSL is designed to provide you with a uniform data structure so that you diff --git a/doc/whats_new.html b/doc/whats_new.html index e517a3e..89c5173 100644 --- a/doc/whats_new.html +++ b/doc/whats_new.html @@ -12,6 +12,8 @@

What's new?

+

Version 1.44:

+Mainly internal improvements. See CHANGES.

Version 1.43:

Mostly bug fixes. See CHANGES.

diff --git a/examples/qlook.c b/examples/qlook.c index f304da0..cf85a63 100644 --- a/examples/qlook.c +++ b/examples/qlook.c @@ -268,6 +268,7 @@ main(int argc, char **argv) { radar->h.year, radar->h.month, radar->h.day, radar->h.hour, radar->h.minute); + if (radar->h.vcp > 0) printf("VCP %d\n", radar->h.vcp); /* Print the radar/volume info. */ diff --git a/gzip.c b/gzip.c index 5b8f902..8308af4 100644 --- a/gzip.c +++ b/gzip.c @@ -84,6 +84,7 @@ FILE *uncompress_pipe (FILE *fp) close(0); dup(save_fd); close(save_fd); + fclose(fp); return fpipe; } diff --git a/nsig.h b/nsig.h index b742ae7..a7e246a 100644 --- a/nsig.h +++ b/nsig.h @@ -42,6 +42,8 @@ #define NSIG_DTB_HCLASS2 56 #define NSIG_DTB_ZDRC 57 #define NSIG_DTB_ZDRC2 58 +#define NSIG_DTB_DBTE8 71 +#define NSIG_DTB_DBZE8 73 /* Product type code ,value for byte 12 in product configuration * struct, III-35 diff --git a/nsig_to_radar.c b/nsig_to_radar.c index d134b55..da9a8ef 100644 --- a/nsig_to_radar.c +++ b/nsig_to_radar.c @@ -174,9 +174,10 @@ RSL_nsig_to_radar float max_vel, sweep_rate, azim_rate; float ray_data; float az1, az2; + float elev1, elev2; double tmp; float sqi, log, csr, sig, cal_dbz; - char radar_type[50], state[2], city[15]; + char radar_type[50], state[4], city[15]; char site_name[16]; NSIG_Product_file *prod_file; short id; @@ -190,6 +191,7 @@ RSL_nsig_to_radar twob nsig_twob; Sweep *sweep; int msec; + const unsigned short low10bits = 0x3ff; float azm, elev, pitch, roll, heading, azm_rate, elev_rate, pitch_rate, roll_rate, heading_rate, lat, lon; @@ -527,7 +529,7 @@ RSL_nsig_to_radar sweep_day = NSIG_I2(nsig_sweep[itype]->idh.time.day); sweep_sec = NSIG_I4(nsig_sweep[itype]->idh.time.sec); #ifdef NSIG_VER2 - msec = NSIG_I2(nsig_sweep[itype]->idh.time.msec); + msec = NSIG_I2(nsig_sweep[itype]->idh.time.msec) & low10bits; /* printf("....... msec == %d\n", msec); */ #endif /* converting seconds since mid to time of day */ @@ -623,15 +625,27 @@ RSL_nsig_to_radar ifield = HC_INDEX; f = HC_F; invf = HC_INVF; + break; case NSIG_DTB_DBZ2: ifield = CZ_INDEX; f = CZ_F; invf = CZ_INVF; + break; case NSIG_DTB_ZDRC2: ifield = ZD_INDEX; f = ZD_F; invf = ZD_INVF; break; + case NSIG_DTB_DBTE8: + ifield = ET_INDEX; + f = ZT_F; + invf = ZT_INVF; + break; + case NSIG_DTB_DBZE8: + ifield = EZ_INDEX; + f = DZ_F; + invf = DZ_INVF; + break; default: fprintf(stderr,"Unknown field type: %d Skipping it.\n", data_type); continue; @@ -736,7 +750,7 @@ RSL_nsig_to_radar ray->h.wavelength = wave/100.0; /* meters */ ray->h.nyq_vel = max_vel; /* m/s */ if (elev == 0.) elev = sweep->h.elev; - ray->h.elev = (nsig_from_bang(ray_p->h.end_elev)+nsig_from_bang(ray_p->h.beg_elev))/2.0; + /* Compute mean azimuth angle for ray. */ az1 = nsig_from_bang(ray_p->h.beg_azm); az2 = nsig_from_bang(ray_p->h.end_azm); @@ -752,6 +766,21 @@ RSL_nsig_to_radar if (az1 > 360) az1 -= 360; ray->h.azimuth = az1; + /* Compute mean elevation angle for ray. */ + elev1 = nsig_from_bang(ray_p->h.beg_elev); + elev2 = nsig_from_bang(ray_p->h.end_elev); + /* printf("elev1, %f, elev2 %f\n", elev1, elev2); */ + if(elev1 > elev2) + if((elev1 - elev2) > 180.0) elev2 += 360.0; + else + ; + else + if((elev2 - elev1) > 180.0) elev2 -= 360.0; + + elev1 = (elev1 + elev2) / 2.0; + if (elev1 > 360) elev1 -= 360; + ray->h.elev = elev1; + /* From the extended header information, we learn the following. */ ray->h.pitch = pitch; ray->h.roll = roll; @@ -782,6 +811,8 @@ RSL_nsig_to_radar switch(data_type) { case NSIG_DTB_UCR: case NSIG_DTB_CR: + case NSIG_DTB_DBTE8: + case NSIG_DTB_DBZE8: if (ray_p->range[k] == 0) ray_data = NSIG_NO_ECHO; else ray_data = (float)((ray_p->range[k]-64.0)/2.0); break; diff --git a/radar_to_uf.c b/radar_to_uf.c index 0c9ee3d..dcd61a0 100644 --- a/radar_to_uf.c +++ b/radar_to_uf.c @@ -63,6 +63,8 @@ typedef short UF_buffer[20000]; /* Bigger than documented 4096. */ void swap_uf_buffer(UF_buffer uf); void swap2(short *buf, int n); +Radar *wsr88d_align_split_cut_rays(Radar *radar); + /**********************************************************************/ /* */ @@ -132,11 +134,24 @@ void RSL_radar_to_uf_fp(Radar *r, FILE *fp) int sweep_num, ray_num, rec_num; float x; + Radar *r_save = NULL; + if (r == NULL) { fprintf(stderr, "radar_to_uf_fp: radar pointer NULL\n"); return; } + /* If this is WSR-88D data, reorder VR and SW rays in split cuts so that + * their azimuths agree with DZ rays. + */ + if (r->h.vcp > 0) { + if (wsr88d_merge_split_cuts_is_set()) { + /* Save a copy of the input radar; we'll restore it later. */ + r_save = r; + r = wsr88d_align_split_cut_rays(r); + } + } + /* Do all the headers first time around. Then, prune OP and LU. */ q_op = q_lu = q_dh = q_fh = 1; @@ -339,9 +354,9 @@ void RSL_radar_to_uf_fp(Radar *r, FILE *fp) 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; + uf_op[6] = r->h.hour; + uf_op[7] = r->h.minute; + uf_op[8] = r->h.sec; memcpy(&uf_op[9], "RADAR_UF", 8); if (little_endian()) swap2(&uf_op[9], 8/2); uf_op[13] = 2; @@ -351,38 +366,11 @@ void RSL_radar_to_uf_fp(Radar *r, FILE *fp) /* ---- Begining of LOCAL USE HEADER BLOCK. */ q_lu = 0; - - /* TODO: Code within "#ifdef LUHDR_VR_AZ" below should be removed - * once testing of merge_split_cuts is completed. - */ -/* 7/25/2011 There's really no need to remove this code after merge_split_cuts. - * We can continue using it for exact VR azimuth in RSL. Non-RSL programs can - * use the mandatory header azimuth, which is exact for DZ and a close - * approximation for VR when merge-split-cuts is applied. - */ -/* Preprocessor directives are for testing. */ -#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; + q_lu = 0; /* Only once. */ + len_lu = 0; /* I don't have anything yet. */ } /* ---- End of LOCAL USE HEADER BLOCK. */ @@ -445,7 +433,7 @@ void RSL_radar_to_uf_fp(Radar *r, FILE *fp) 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[9] = ray->h.frequency * 1000.; /* 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; @@ -511,6 +499,9 @@ void RSL_radar_to_uf_fp(Radar *r, FILE *fp) } /* if (ray) */ } } + + /* If Radar argument "r" was modified (WSR-88D), restore the saved copy. */ + if (r_save != NULL) r = r_save; } /**********************************************************************/ diff --git a/rsl.h b/rsl.h index 2e32654..41cb299 100644 --- a/rsl.h +++ b/rsl.h @@ -27,7 +27,7 @@ #include "config.h" #endif -#define RSL_VERSION_STR "v1.43" +#define RSL_VERSION_STR "v1.44" /**********************************************************************/ /* Configure: Define USE_TWO_BYTE_PRECISION to have RSL store internal*/ @@ -46,7 +46,7 @@ /* so you shouldn't have to modify anything here. */ /**********************************************************************/ #ifndef COLORDIR -#define COLORDIR "/usr/local/trmm/lib/colors" +#define COLORDIR "/home/kelley/trmm/lib/colors" #endif /* These are the color table indexes. See RSL_set/get_color_table. */ @@ -402,6 +402,8 @@ typedef struct { *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. + *42 = ET_INDEX = Total Power Enhanced (Sigmet) + *43 = EZ_INDEX = Clutter Corr. Reflectivity Enhanced (Sigmet) */ } Radar; @@ -496,7 +498,7 @@ typedef struct { * rsl_qfield by adding a '1' for each new volume index. */ -#define MAX_RADAR_VOLUMES 42 +#define MAX_RADAR_VOLUMES 44 #define DZ_INDEX 0 #define VR_INDEX 1 @@ -540,6 +542,8 @@ typedef struct { #define SD_INDEX 39 #define ZZ_INDEX 40 #define RD_INDEX 41 +#define ET_INDEX 42 +#define EZ_INDEX 43 /* Prototypes for functions. */ @@ -851,7 +855,8 @@ static char *RSL_ftype[] = {"DZ", "VR", "SW", "CZ", "ZT", "DR", "TI", "DX", "CH", "AH", "CV", "AV", "SQ", "VS", "VL", "VG", "VT", "NP", "HC", "VC", "V2", "S2", "V3", "S3", - "CR", "CC", "PR", "SD", "ZZ", "RD"}; + "CR", "CC", "PR", "SD", "ZZ", "RD", + "ET", "EZ"}; 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, @@ -859,7 +864,8 @@ static float (*RSL_f_list[])(Range x) = {DZ_F, VR_F, SW_F, CZ_F, ZT_F, DR_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}; + DZ_F, CZ_F, PH_F, SD_F, DZ_F, DZ_F, + ZT_F, DZ_F}; static Range (*RSL_invf_list[])(float x) = {DZ_INVF, VR_INVF, SW_INVF, CZ_INVF, ZT_INVF, DR_INVF, @@ -868,7 +874,8 @@ static Range (*RSL_invf_list[])(float x) 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}; + DZ_INVF, CZ_INVF, PH_INVF, SD_INVF, DZ_INVF, DZ_INVF, + ZT_INVF, DZ_INVF}; #endif /* Secret routines that are quite useful and useful to developers. */ void radar_load_date_time(Radar *radar); @@ -883,6 +890,13 @@ Hash_table *construct_sweep_hash_table(Sweep *s); double angle_diff(float x, float y); int rsl_query_field(char *c_field); +/* Functions for controlling handling of WSR-88D split cuts. */ +void RSL_wsr88d_merge_split_cuts_on(); +void RSL_wsr88d_merge_split_cuts_off(); +void RSL_wsr88d_keep_short_refl(); +int wsr88d_merge_split_cuts_is_set(); + +Radar *wsr88d_merge_split_cuts(Radar *radar); /* Debugging prototypes. */ void poke_around_volume(Volume *v); diff --git a/uf_to_radar.c b/uf_to_radar.c index e87b68d..9835a19 100644 --- a/uf_to_radar.c +++ b/uf_to_radar.c @@ -24,6 +24,7 @@ #include #include #include +#include /* This allows us to use RSL_ftype, RSL_f_list, RSL_invf_list from rsl.h. */ #define USE_RSL_VARS @@ -87,6 +88,63 @@ void swap2(short *buf, int n) } } +static void put_start_time_in_radar_header(Radar *radar) +{ + /* Get the earliest ray time and store it in radar header. + * The search is necessary because rays are not always in chronological order. + * For example, we have received data in which rays were apparently sorted by + * azimuth in some upstream software. This results in the ray times being out + * of order, because a sweep rarely actually begins at zero degrees. + * + * Written by Bart Kelley, SSAI, June 19, 2013 + */ + + int i = 0; + Sweep *sweep; + Ray *ray; + + int prevdate, thisdate; + float prevtime, thistime; + + /* Get first sweep of first available field. */ + for (i=0; i < MAX_RADAR_VOLUMES; i++) { + if ((sweep = radar->v[i]->sweep[0]) != NULL) break; + } + /* This shouldn't happen. */ + if (i >= MAX_RADAR_VOLUMES) { + printf("put_start_time_in_radar_header: No radar volumes contained " + "sweep at index 0.\n"); + return; + } + + /* Get first ray and its time. */ + i = 0; + while (!sweep->ray[i] && i < sweep->h.nrays) i++; + ray = sweep->ray[i]; + prevdate = ray->h.year * 10000 + ray->h.month * 100 + ray->h.day; + prevtime = ray->h.hour * 10000 + ray->h.minute * 100 + ray->h.sec; + + /* Compare times of remaining rays for earliest time. */ + for (i=0; ih.nrays; i++) { + ray = sweep->ray[i]; + thisdate = ray->h.year * 10000 + ray->h.month * 100 + ray->h.day; + thistime = ray->h.hour * 10000 + ray->h.minute * 100 + ray->h.sec; + if (thisdate == prevdate) { + if (thistime < prevtime) prevtime = thistime; + } + else if (thisdate < prevdate) { + prevdate = thisdate; + prevtime = thistime; + } + } + + radar->h.year = prevdate / 10000; + radar->h.month = prevdate / 100 % 100; + radar->h.day = prevdate % 100; + radar->h.hour = (int) prevtime / 10000; + radar->h.minute = (int) prevtime / 100 % 100; + radar->h.sec = fmod(prevtime,100.); +} /********************************************************************/ /*********************************************************************/ @@ -150,6 +208,7 @@ int uf_into_radar(UF_buffer uf, Radar **the_radar) short missing_data; Volume *new_volume; int nbins; + float frequency; extern int rsl_qfield[]; extern int *rsl_qsweep; /* See RSL_read_these_sweeps in volume.c */ extern int rsl_qsweep_max; @@ -336,6 +395,10 @@ int uf_into_radar(UF_buffer uf, Radar **the_radar) radar->h.lond = uf_ma[21]; radar->h.lonm = uf_ma[22]; radar->h.lons = uf_ma[23] / 64.0; + /* Note that radar header time is now handled at end of ingest by + * function put_start_time_in_radar_header(). The values below are + * replaced. --BLK, 6/19/13 + */ radar->h.year = ray->h.year; radar->h.month = ray->h.month; radar->h.day = ray->h.day; @@ -392,7 +455,14 @@ int uf_into_radar(UF_buffer uf, Radar **the_radar) sweep->h.vert_half_bw = .5; } - ray->h.frequency = uf_fh[9] / 64.0; + frequency = uf_fh[9]; + /* This corrects an error in v1.43 and earlier where frequency was + * multiplied by 64. Correct units for UF are MHz; radar structure + * uses GHz. + */ + if (frequency < 1000.) frequency = frequency/64.; + else frequency = frequency/1000.; + ray->h.frequency = frequency; 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) { @@ -568,6 +638,7 @@ Radar *RSL_uf_to_radar_fp(FILE *fp) case NOT_UF: return NULL; break; } radar = reset_nsweeps_in_all_volumes(radar); + put_start_time_in_radar_header(radar); radar = RSL_prune_radar(radar); return radar; diff --git a/volume.c b/volume.c index a30368b..c2072ba 100644 --- a/volume.c +++ b/volume.c @@ -89,10 +89,13 @@ extern int radar_verbose_flag; * The conversion functions may NOT be macros. */ +/* Modification to properly convert DM field (DEK 21 Nov 2012) */ + #ifdef USE_TWO_BYTE_PRECISION #define F_FACTOR 100.0 #define F_DR_FACTOR 1000.0 #define F_DZ_RANGE_OFFSET 50 +#define F_DM_RANGE_OFFSET 130 /* DEK */ #else #define F_FACTOR 2.0 #define F_DR_FACTOR 10.0 @@ -122,6 +125,18 @@ float DZ_F(Range x) { return BADVAL; /* Can't get here, but quiets the compiler. */ } + +float DM_F(Range x) { + if (x >= F_OFFSET) /* This test works when Range is unsigned. */ + return (((float)x-F_OFFSET)/F_FACTOR - F_DM_RANGE_OFFSET); /* DEK */ + if (x == 0) return BADVAL; + if (x == 1) return RFVAL; + if (x == 2) return APFLAG; + if (x == 3) return NOECHO; + return BADVAL; /* Can't get here, but quiets the compiler. */ +} + + float VR_F(Range x) { float val; if (x >= F_OFFSET) { /* This test works when Range is unsigned. */ @@ -177,7 +192,7 @@ float HC_F(Range x) { /* HydroClass (Sigmet) */ float RH_F(Range x) { if (x == 0) return BADVAL; /* return (float)(sqrt((double)((x-1.0)/253.0))); */ - return (float)(x-1) / 65533.; + return (float)x / 100.; } /***************************** @@ -193,7 +208,7 @@ float RH_F(Range x) { float PH_F(Range x) { if (x == 0) return BADVAL; /*return (float)(180.0*((x-1.0)/254.0));*/ - return (360.*(x-1.))/65534.; + return x/10.; } /* TODO: Should this be 5. cm. instead of 0.5? Or maybe 10. cm.? */ @@ -203,29 +218,8 @@ float rsl_kdp_wavelen = 0.5; /* Default radar wavelen = .5 cm. See /* KD_F for 1 or 2 byte. */ float KD_F(Range x) { -/****** Commented-out code for 1-byte Sigmet native data format. - if (x >= F_OFFSET) { - x -= F_OFFSET; - if (rsl_kdp_wavelen == 0.0) return BADVAL; - if (x < 128) - return (float)( - -0.25 * pow((double)600.0,(double)((127-x)/126.0)) - )/rsl_kdp_wavelen; - else if (x > 128) - return (float)( - 0.25 * pow((double)600.0,(double)((x-129)/126.0)) - )/rsl_kdp_wavelen; - else - return 0.0; - } - if (x == 1) return RFVAL; - if (x == 2) return APFLAG; - if (x == 3) return NOECHO; - return BADVAL; -******/ - if (x == 0) return BADVAL; - return (x-32768.)/100.; + return x/100.; } /* Normalized Coherent Power (DORADE) */ @@ -269,7 +263,7 @@ float MZ_F(Range x) { return (float)x; } /* DZ Mask */ float MD_F(Range x) { return MZ_F(x); } /* ZD Mask */ float ZE_F(Range x) { return DZ_F(x); } float VE_F(Range x) { return VR_F(x); } -float DM_F(Range x) { return DZ_F(x); } +/* float DM_F(Range x) { return DZ_F(x); } DEK */ float DX_F(Range x) { return DZ_F(x); } float CH_F(Range x) { return DZ_F(x); } float AH_F(Range x) { return DZ_F(x); } @@ -298,6 +292,16 @@ Range DZ_INVF(float x) return (Range)(F_FACTOR*(x+F_DZ_RANGE_OFFSET)+.5 + F_OFFSET); /* Default wsr88d. */ } +Range DM_INVF(float x) +{ + if (x == BADVAL) return (Range)0; + if (x == RFVAL) return (Range)1; + if (x == APFLAG) return (Range)2; + if (x == NOECHO) return (Range)3; + return (Range)(F_FACTOR*(x+F_DM_RANGE_OFFSET)+.5 + F_OFFSET); /* DEK */ +} + + Range VR_INVF(float x) { if (x == BADVAL) return (Range)0; @@ -347,7 +351,7 @@ Range LR_INVF(float x) /* MCTEX */ Range RH_INVF(float x) { if (x == BADVAL) return (Range)0; /* return (Range)(x * x * 253.0 + 1.0 + 0.5); */ - return (Range)(x * 65533. + 1. +.5); + return (Range)(x * 100. +.5); } /****************************** @@ -363,34 +367,14 @@ Range RH_INVF(float x) { Range PH_INVF(float x) { if (x == BADVAL) return (Range)0; /* return (Range)((x / 180.0) * 254.0 + 1.0 + 0.5); */ - return (Range)(x*65534./360. + 1.0 + 0.5); + return (Range)(x*10.+ 0.5); } /* KD_INVF for 1 or 2 byte data. */ Range KD_INVF(float x) { if (x == BADVAL) return (Range)0; - 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; - } else if (x > 0) { - x = 129 + - 126 * (log((double)x) - log((double)0.25/rsl_kdp_wavelen)) / - log((double)600.0) + - 0.5; - } else { - x = 128; - } - x += F_OFFSET; -******/ + return (Range)(x * 100. + 0.5); } /* Standard Deviation (for Dual-pole QC testing.) */ @@ -434,7 +418,7 @@ Range MZ_INVF(float x) { return (Range)x; } /* DZ Mask */ Range MD_INVF(float x) { return MZ_INVF(x); } /* ZD Mask */ Range ZE_INVF(float x) { return DZ_INVF(x); } Range VE_INVF(float x) { return VR_INVF(x); } -Range DM_INVF(float x) { return DZ_INVF(x); } +/* Range DM_INVF(float x) { return DZ_INVF(x); } DEK */ Range DX_INVF(float x) { return DZ_INVF(x); } Range CH_INVF(float x) { return DZ_INVF(x); } Range AH_INVF(float x) { return DZ_INVF(x); } @@ -1404,7 +1388,7 @@ 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 }; diff --git a/wsr88d.h b/wsr88d.h index 62490fc..de996fc 100644 --- a/wsr88d.h +++ b/wsr88d.h @@ -7,7 +7,7 @@ * The directory should be the same as the LIBDIR in the makefile. */ #ifndef WSR88D_SITE_INFO_FILE -#define WSR88D_SITE_INFO_FILE "/usr/local/trmm/lib/wsr88d_locations.dat" +#define WSR88D_SITE_INFO_FILE "/home/kelley/trmm/lib/wsr88d_locations.dat" #endif /*===============================================================*/ typedef struct { @@ -139,8 +139,7 @@ typedef struct radar_site { } Wsr88d_site_info; typedef struct { - FILE *fptr; /* this usually points to the gzip pipe */ - FILE *orig; /* save the original file pointer for cleanup */ + FILE *fptr; } Wsr88d_file; #define PACKET_SIZE 2432 @@ -211,8 +210,6 @@ float wsr88d_get_prf(Wsr88d_ray *ray); float wsr88d_get_prt(Wsr88d_ray *ray); float wsr88d_get_wavelength(Wsr88d_ray *ray); float wsr88d_get_frequency(Wsr88d_ray *ray); -void wsr88d_keep_hi_prf_dz(); -void wsr88d_no_hi_prf_dz(); int no_command (char *cmd); FILE *uncompress_pipe (FILE *fp); diff --git a/wsr88d_align_split_cut_rays.c b/wsr88d_align_split_cut_rays.c new file mode 100644 index 0000000..12fdd00 --- /dev/null +++ b/wsr88d_align_split_cut_rays.c @@ -0,0 +1,428 @@ +#include +#include +#include +#ifndef _rsl_h +#include "rsl.h" +#endif + +void replace_ray_gaps_with_null(Sweep *sweep) +{ + /* Check for ray gaps and replace any with null. A ray gap means that + * the difference in ray numbers in a pair of consecutive rays is greater + * than 1. When such a case is found, nulls are inserted for missing rays. + */ + + int i, j, maxrays, nrays; + Sweep *newsweep; + + if (sweep == NULL) { + fprintf(stderr,"replace_ray_gaps_with_null: Null sweep.\n"); + return; + } + + nrays = sweep->h.nrays; + /* If last ray of sweep is non-null, and its ray number is same as nrays, + * then there are no gaps. + */ + if (sweep->ray[nrays-1] && sweep->ray[nrays-1]->h.ray_num == nrays) return; + + /* Determine expected maximum rays from beamwidth. */ + maxrays = 720; + if (sweep->h.beam_width > 0.9) maxrays = 360; + if (sweep->h.nrays > maxrays) maxrays = sweep->h.nrays; + newsweep = RSL_new_sweep(maxrays); + + /* Copy rays to new sweep based on ray number. New sweep will contain + * NULL where rays are missing. + */ + for (i=0, j=0; iray[j] && sweep->ray[j]->h.ray_num == i+1) + newsweep->ray[i] = sweep->ray[j++]; + + /* Copy rays back to original sweep. */ + for (i=0; iray[i] = newsweep->ray[i]; + + sweep->h.nrays = maxrays; + free(newsweep); +} + + +int get_ray_index(Sweep *sweep, Ray *ray) +{ + /* Find the array index for the given ray in the sweep structure. */ + + float target_azimuth, this_azimuth; + int iray, incr, maxrays, nchk; + + if (ray == NULL) return -1; + + target_azimuth = ray->h.azimuth; + + /* Start with ray_num - 1. That will normally be a match. */ + iray = ray->h.ray_num - 1; + if (sweep->ray[iray]->h.azimuth == target_azimuth) + return iray; + + /* Search for ray with matching azimuth by incrementing ray index forward + * or backward, depending on azimuth. + */ + this_azimuth = sweep->ray[iray]->h.azimuth; + if (this_azimuth < target_azimuth) incr = 1; + else incr = -1; + iray+=incr; + maxrays = sweep->h.nrays; + for (nchk = 1; nchk <= maxrays; nchk++) { + if (iray < 0) iray = maxrays - 1; + if (iray >= maxrays) iray = 0; + if (sweep->ray[iray]->h.azimuth == target_azimuth) return iray; + iray+=incr; + } + /* Control shouldn't end up here. Ray came from this sweep, so there + * has to be a match. + */ + fprintf(stderr,"get_ray_index: Didn't find matching azimuth, but there " + "should be one.\n"); + fprintf(stderr,"Current values:\n"); + fprintf(stderr, "target_azimuth: %.2f\n" + "this_azimuth: %.2f\n" + "iray: %d\n" + "incr: %d\n" + "maxrays: %d\n" + "nchk: %d\n", + target_azimuth, this_azimuth, iray, incr, maxrays, nchk); + return -1; +} + + +int get_azimuth_match(Sweep *sweep1, Sweep *sweep2, int *iray1, int *iray2) +{ + /* Using the azimuth of sweep1->ray[iray1], find the ray in sweep2 with + * matching azimuth. The index of matching ray is returned through + * argument iray2. + * + * Return value: function returns 1 for successful match, 0 otherwise. + */ + + Ray *ray2; + float matchlimit; + const notfound = 0, found = 1; + + matchlimit = sweep2->h.horz_half_bw; + /* Assure that matchlimit is wide enough--a problem with pre-Build10. */ + if (sweep2->h.beam_width > .9) matchlimit = 0.5; + + ray2 = RSL_get_closest_ray_from_sweep(sweep2, + sweep1->ray[*iray1]->h.azimuth, matchlimit); + + /* If no match, try again with a slightly wider limit. */ + if (ray2 == NULL) ray2 = RSL_get_closest_ray_from_sweep(sweep2, + sweep1->ray[*iray1]->h.azimuth, matchlimit + .04 * matchlimit); + + if (ray2 == NULL) { + *iray2 = -1; + return notfound; + } + + *iray2 = get_ray_index(sweep2, ray2); + return found; +} + + +int get_first_azimuth_match(Sweep *dzsweep, Sweep *vrsweep, int *dz_ray_index, + int *vr_ray_index) +{ + /* Find the DZ ray whose azimuth matches the first ray of VR sweep. + * Indexes for the matching rays are returned through arguments. + * Return value: 1 for successful match, 0 otherwise. + */ + + Ray *dz_ray, *vr_ray; + int iray_dz = 0, iray_vr = 0, match_found = 0; + const notfound = 0, found = 1; + + /* Get first non-null VR ray. */ + vr_ray = vrsweep->ray[iray_vr]; + while (vr_ray == NULL && iray_vr < vrsweep->h.nrays) + iray_vr++; + if (vr_ray == NULL) { + fprintf(stderr,"get_first_azimuth_match: All rays in VR sweep number %d" + " are NULL.\nMerge split cut not done for this sweep.\n", + vrsweep->h.sweep_num); + return notfound; + } + + /* Get DZ ray with azimuth closest to first VR ray. */ + + while (!match_found) { + if (get_azimuth_match(vrsweep, dzsweep, &iray_vr, &iray_dz)) + match_found = 1; + else { + /* If still no match, get next non-null VR ray. */ + iray_vr++; + while ((vr_ray = vrsweep->ray[iray_vr]) == NULL) { + iray_vr++; + if (iray_vr >= vrsweep->h.nrays) { + fprintf(stderr,"get_first_azimuth_match: Couldn't find a " + "matching azimuth.\nSweeps compared:\n"); + fprintf(stderr,"DZ: sweep number %d, elev %f\n" + "VR: sweep number %d, elev %f\n", + dzsweep->h.sweep_num, dzsweep->h.elev, + vrsweep->h.sweep_num, vrsweep->h.elev); + return notfound; + } + } + } + } + + *vr_ray_index = iray_vr; + *dz_ray_index = iray_dz; + return found; +} + + +void reorder_rays(Sweep *dzsweep, Sweep *vrsweep, Sweep *swsweep) +{ + /* Reorder rays of the VR and SW sweeps so that they match by azimuth + * the rays in the DZ sweep. + */ + int i, iray_dz, iray_vr, iray_dz_start, ray_num; + int maxrays; + int more_rays; + Sweep *new_vrsweep, *new_swsweep; + + /* Allocate new VR and SW sweeps. Allocate maximum number of rays based + * on beam width. + */ + + if (vrsweep != NULL) { + maxrays = (vrsweep->h.beam_width < 0.8) ? 720 : 360; + if (vrsweep->h.nrays > maxrays) maxrays = vrsweep->h.nrays; + } + else { + fprintf(stderr,"reorder_rays: Velocity sweep number %d is NULL.\n" + "Merge split cut not done for this sweep.\n", + vrsweep->h.sweep_num); + return; + } + + replace_ray_gaps_with_null(dzsweep); + replace_ray_gaps_with_null(vrsweep); + if (swsweep) replace_ray_gaps_with_null(swsweep); + + /* Get ray indices for first match of DZ and VR azimuths. */ + if (get_first_azimuth_match(dzsweep, vrsweep, &iray_dz, &iray_vr) == 0) { + fprintf(stderr,"reorder_rays: Merge split cut not done for DZ sweep " + "number %d and VR sweep %d.\n", + dzsweep->h.sweep_num, vrsweep->h.sweep_num); + return; + } + + /* If DZ and VR azimuths already match, no need to do ray alignment. */ + if (iray_dz == iray_vr) return; + + new_vrsweep = RSL_new_sweep(maxrays); + memcpy(&new_vrsweep->h, &vrsweep->h, sizeof(Sweep_header)); + if (swsweep) { + new_swsweep = RSL_new_sweep(maxrays); + memcpy(&new_swsweep->h, &swsweep->h, sizeof(Sweep_header)); + } + + /* Copy VR rays to new VR sweep such that their azimuths match those of + * DZ rays at the same indices. Do the same for SW. + */ + + iray_dz_start = iray_dz; + ray_num = dzsweep->ray[iray_dz]->h.ray_num; + more_rays = 1; + while (more_rays) { + if (vrsweep->ray[iray_vr]) { + new_vrsweep->ray[iray_dz] = vrsweep->ray[iray_vr]; + new_vrsweep->ray[iray_dz]->h.ray_num = ray_num; + } + if (swsweep && swsweep->ray[iray_vr]) { + new_swsweep->ray[iray_dz] = swsweep->ray[iray_vr]; + new_swsweep->ray[iray_dz]->h.ray_num = ray_num; + } + iray_dz++; + iray_vr++; + ray_num++; + + /* Handle sweep wraparound. */ + if (iray_dz == dzsweep->h.nrays || iray_dz == maxrays) { + iray_dz = 0; + ray_num = 1; + /* Check if azimuth rematch is needed. */ + if (dzsweep->ray[iray_dz] && vrsweep->ray[iray_vr] && + fabsf(dzsweep->ray[iray_dz]->h.azimuth - + vrsweep->ray[iray_vr]->h.azimuth) + > dzsweep->h.horz_half_bw) { + get_azimuth_match(dzsweep, vrsweep, &iray_dz, &iray_vr); + } + } + if (iray_vr == vrsweep->h.nrays) iray_vr = 0; + + /* If we're back at first match, we're done. */ + if (iray_dz == iray_dz_start) more_rays = 0; + } + + /* Copy ray pointers in new order back to original sweeps. */ + /* TODO: Can we simply replace the sweeps themselves with the new sweeps? */ + + for (i=0; i < maxrays; i++) { + vrsweep->ray[i] = new_vrsweep->ray[i]; + if (swsweep) swsweep->ray[i] = new_swsweep->ray[i]; + } + + /* TODO: + * Update nrays in sweep header with ray number of last ray + * (ray count includes NULL rays). + */ + vrsweep->h.nrays = maxrays; + if (swsweep) swsweep->h.nrays = maxrays; +} + + +void align_this_split_cut(Radar *radar, int iswp) +{ + int vrindex[] = {VR_INDEX, V2_INDEX, V3_INDEX}; + int swindex[] = {SW_INDEX, S2_INDEX, S3_INDEX}; + int ivr, nvrfields; + Sweep *dzsweep, *vrsweep, *swsweep; + + dzsweep = radar->v[DZ_INDEX]->sweep[iswp]; + if (dzsweep == NULL) { + fprintf(stderr,"align_this_split_cut: DZ sweep at index %d is NULL.\n" + "Merge split cut not done for this sweep.\n", iswp); + return; + } + + if (radar->h.vcp != 121) nvrfields = 1; + else { + if (iswp < 4) nvrfields = 3; + else nvrfields = 2; + } + + for (ivr = 0; ivr < nvrfields; ivr++) { + if (radar->v[vrindex[ivr]]) { + vrsweep = radar->v[vrindex[ivr]]->sweep[iswp]; + } + else { + fprintf(stderr,"align_this_split_cut: No VR at INDEX %d.\n" + "No ray alignment done for this VR index at DZ sweep " + "index %d.\n", vrindex[ivr], iswp); + return; + } + /* Make sure vrsweep and dzsweep elevs are approximately the same. */ + if (fabsf(dzsweep->h.elev - vrsweep->h.elev) > 0.1) { + fprintf(stderr,"align_this_split_cut: elevation angles for DZ and " + "VR don't match:\n" + "radar->v[DZ_INDEX=%d]->sweep[%d]->h.elev = %f\n" + "radar->v[VR_INDEX=%d]->sweep[%d]->h.elev = %f\n", + DZ_INDEX, iswp, dzsweep->h.elev, + VR_INDEX, iswp, vrsweep->h.elev); + return; + } + if (radar->v[swindex[ivr]]) + swsweep = radar->v[swindex[ivr]]->sweep[iswp]; + else + swsweep = NULL; + reorder_rays(dzsweep, vrsweep, swsweep); + } +} + + +static Radar *copy_radar_for_split_cut_adjustments(Radar *r) +{ + /* Return a copy of the radar structure for modification at split cuts. */ + + Radar *r_new; + int i,j,k; + + r_new = RSL_new_radar(MAX_RADAR_VOLUMES); + r_new->h = r->h; + /* Copy the contents of VR and SW volumes, which will be modified. For other + * fields, just copy pointers. + */ + for (k=0; k < MAX_RADAR_VOLUMES; k++) { + if (r->v[k] == NULL) continue; + switch(k) { + case VR_INDEX: case SW_INDEX: case V2_INDEX: case S2_INDEX: + case V3_INDEX: case S3_INDEX: + r_new->v[k] = RSL_new_volume(r->v[k]->h.nsweeps); + r_new->v[k]->h = r->v[k]->h; + for (i=0; i < r->v[k]->h.nsweeps; i++) { + /* Add new sweeps to volume and copy ray pointers. */ + r_new->v[k]->sweep[i] = RSL_new_sweep(r->v[k]->sweep[i]->h.nrays); + r_new->v[k]->sweep[i]->h = r->v[k]->sweep[i]->h; + for (j=0; j < r->v[k]->sweep[i]->h.nrays; j++) { + r_new->v[k]->sweep[i]->ray[j] = r->v[k]->sweep[i]->ray[j]; + } + } + break; + default: + /* For fields other than VR or SW types, simply copy volume pointer. + */ + r_new->v[k] = r->v[k]; + break; + } + } + return r_new; +} + + +Radar *wsr88d_align_split_cut_rays(Radar *radar_in) +{ +/* This is the driver function for aligning rays in split cuts by azimuth. */ + + int iswp; + int nsplit, vcp; + Radar *radar; + + /* If VR is not present, nothing more to do. */ + if (radar_in->v[VR_INDEX] == NULL) { + /* If SW without VR, notify user that VR must be included. */ + if (radar_in->v[SW_INDEX]) + fprintf(stderr,"wsr88d_align_split_cut_rays: VR must be selected in" + " order to do ray alignment\n" + " on SW. No ray alignment for SW in split cuts.\n"); + return radar_in; + } + + vcp = radar_in->h.vcp; + + /* How many split cuts are in this VCP? */ + nsplit = 0; + switch (vcp) { + case 11: case 21: case 32: case 211: case 221: + nsplit = 2; + break; + case 12: case 31: case 212: + nsplit = 3; + break; + case 121: + nsplit = 5; + break; + } + + if (nsplit == 0) { + fprintf(stderr, "wsr88d_align_split_cut_rays: Unknown VCP: %d\n", vcp); + return radar_in; + } + + /* Make a copy of radar. It is the copy that we'll modify. */ + radar = copy_radar_for_split_cut_adjustments(radar_in); + + /* Align rays in each split cut. */ + + for (iswp = 0; iswp < nsplit; iswp++) { + if (radar->v[DZ_INDEX]->sweep[iswp] == NULL) { + fprintf(stderr, "wsr88d_align_split_cut_rays: Encountered " + "unexpected NULL DZ sweep\n"); + fprintf(stderr, "iswp (sweep index) = %d\n", iswp); + return radar; + } + align_this_split_cut(radar, iswp); + } + + return radar; +} diff --git a/wsr88d_locations.dat b/wsr88d_locations.dat index cc67354..d2f9304 100644 --- a/wsr88d_locations.dat +++ b/wsr88d_locations.dat @@ -1,24 +1,24 @@ 14929 KABR ABERDEEN SD 45 27 21 -98 24 47 397 -3019 KABX ALBUQUERQUE NM 35 8 59 -106 49 26 1789 +03019 KABX ALBUQUERQUE NM 35 8 59 -106 49 26 1789 93773 KAKQ NORFOLK/RICH VA 36 59 2 -77 0 26 34 23047 KAMA AMARILLO TX 35 14 0 -101 42 33 1093 -12899 KAMX MIAMI FL 25 36 40 -80 24 46 4 -4837 KAPX GAYLORD MI 44 54 26 -84 43 11 446 +12899 KAMX MIAMI FL 25 36 40 -80 24 46 4 +04837 KAPX GAYLORD MI 44 54 26 -84 43 11 446 94987 KARX LA_CROSSE WI 43 49 22 -91 11 28 389 -94287 KATX SEATTLE WA 48 11 40 -122 29 45 151 +94287 KATX SEATTLE WA 48 11 40 -122 29 45 151 93240 KBBX BEALE_AFB CA 39 29 46 -121 37 54 53 -4725 KBGM BINGHAMTON NY 42 11 59 -75 59 5 490 -94289 KBHX EUREKA CA 40 29 54 -124 17 31 732 +04725 KBGM BINGHAMTON NY 42 11 59 -75 59 5 490 +94289 KBHX EUREKA CA 40 29 54 -124 17 31 732 24011 KBIS BISMARCK ND 46 46 15 -100 45 38 505 94046 KBLX BILLINGS MT 45 51 14 -108 36 24 1097 53823 KBMX BIRMINGHAM AL 33 10 20 -86 46 12 197 -54765 KBOX BOSTON MA 41 57 21 -71 8 13 36 +54765 KBOX BOSTON MA 41 57 21 -71 8 13 36 12919 KBRO BROWNSVILLE TX 25 54 58 -97 25 8 7 -14733 KBUF BUFFALO NY 42 56 56 -78 44 12 211 +14733 KBUF BUFFALO NY 42 56 56 -78 44 12 211 92804 KBYX KEY_WEST FL 24 35 51 -81 42 11 3 13883 KCAE COLUMBIA SC 33 56 55 -81 7 6 70 -94625 KCBW HOULTON ME 46 2 22 -67 48 24 227 -4101 KCBX BOISE ID 43 29 27 -116 14 8 933 +94625 KCBW HOULTON ME 46 2 22 -67 48 24 227 +04101 KCBX BOISE ID 43 29 27 -116 14 8 933 54764 KCCX STATE_COLLEGE PA 40 55 23 -78 0 13 733 14820 KCLE CLEVELAND OH 41 24 47 -81 51 35 233 53845 KCLX CHARLESTON SC 32 39 20 -81 2 31 30 @@ -28,129 +28,138 @@ 93235 KDAX SACRAMENTO CA 38 30 4 -121 40 40 9 13985 KDDC DODGE_CITY KS 37 45 39 -99 58 7 789 22015 KDFX LAUGHLIN_AFB TX 29 16 22 -100 16 50 345 +33771 KDGX JACKSON MS 32 16 48 -89 59 04 151 93771 KDIX PHILADELPHIA PA 39 56 49 -74 24 39 45 -14913 KDLH DULUTH MN 46 50 13 -92 12 35 435 +14913 KDLH DULUTH MN 46 50 13 -92 12 35 435 94984 KDMX DES_MOINES IA 41 43 53 -93 43 22 299 93770 KDOX DOVER_AFB DE 38 49 32 -75 26 24 15 -4830 KDTX DETROIT MI 42 41 59 -83 28 18 327 +04830 KDTX DETROIT MI 42 41 59 -83 28 18 327 94982 KDVN DAVENPORT IA 41 36 42 -90 34 51 230 -3987 KDYX DYESS_AFB TX 32 32 18 -99 15 15 462 -3983 KEAX KANSAS_CITY MO 38 48 37 -94 15 52 303 -53112 KEMX TUCSON AZ 31 53 37 -110 37 49 1586 -54766 KENX ALBANY NY 42 35 11 -74 3 50 557 +03987 KDYX DYESS_AFB TX 32 32 18 -99 15 15 462 +03983 KEAX KANSAS_CITY MO 38 48 37 -94 15 52 303 +53112 KEMX TUCSON AZ 31 53 37 -110 37 49 1586 +54766 KENX ALBANY NY 42 35 11 -74 3 50 557 53851 KEOX FORT_RUCKER AL 31 27 38 -85 27 34 132 -3020 KEPZ EL_PASO TX 31 52 23 -106 41 53 1251 +03020 KEPZ EL_PASO TX 31 52 23 -106 41 53 1251 53110 KESX LAS_VEGAS NV 35 42 4 -114 53 29 1483 53825 KEVX EGLIN_AFB FL 30 33 52 -85 55 17 43 12971 KEWX AUSTIN/S_ANT TX 29 42 14 -98 1 42 193 53114 KEYX EDWARDS_AFB CA 35 5 52 -117 33 39 840 -53831 KFCX ROANOKE VA 37 1 28 -80 16 26 874 -3981 KFDR ALTUS_AFB OK 34 21 44 -98 58 35 386 -3022 KFDX CANNON_AFB NM 34 38 7 -103 37 48 1417 -53819 KFFC ATLANTA GA 33 21 49 -84 33 57 262 +53831 KFCX ROANOKE VA 37 1 28 -80 16 26 874 +03981 KFDR ALTUS_AFB OK 34 21 44 -98 58 35 386 +03022 KFDX CANNON_AFB NM 34 38 7 -103 37 48 1417 +53819 KFFC ATLANTA GA 33 21 49 -84 33 57 262 14944 KFSD SIOUX_FALLS SD 43 35 16 -96 43 46 436 53113 KFSX FLAGSTAFF AZ 34 34 28 -111 11 52 2261 -3018 KFTG DENVER CO 39 47 12 -104 32 45 1675 -3985 KFWS DALLAS/FTW TX 32 34 23 -97 18 11 208 -94008 KGGW GLASGOW MT 48 12 23 -106 37 30 694 -3025 KGJX GRAND_JUNCT CO 39 3 44 -108 12 50 3046 +03018 KFTG DENVER CO 39 47 12 -104 32 45 1675 +03985 KFWS DALLAS/FTW TX 32 34 23 -97 18 11 208 +94008 KGGW GLASGOW MT 48 12 23 -106 37 30 694 +03025 KGJX GRAND_JUNCT CO 39 3 44 -108 12 50 3046 23065 KGLD GOODLAND KS 39 21 59 -101 42 2 1113 14898 KGRB GREEN_BAY WI 44 29 54 -88 6 41 208 -3992 KGRK FORT_HOOD TX 30 43 19 -97 22 59 164 +03992 KGRK FORT_HOOD TX 30 43 19 -97 22 59 164 94860 KGRR GRAND_RAPIDS MI 42 53 38 -85 32 41 237 -3870 KGSP GREER SC 34 52 0 -82 13 12 287 +03870 KGSP GREER SC 34 53 0 -82 13 12 287 53837 KGWX COLUMBUS_AFB MS 33 53 48 -88 19 44 145 54762 KGYX PORTLAND ME 43 53 29 -70 15 24 125 -3023 KHDX HOLLOMAN_AFB NM 33 4 35 -106 7 22 1287 -3980 KHGX HOUSTON TX 29 28 19 -95 4 45 5 +03023 KHDX HOLLOMAN_AFB NM 33 4 35 -106 7 22 1287 +03980 KHGX HOUSTON TX 29 28 19 -95 4 45 5 53111 KHNX SAN_JOAQUIN_V CA 36 18 51 -119 37 56 75 53839 KHPX FORT_CAMPBELL KY 36 44 12 -87 17 6 176 53857 KHTX HUNTSVILLE AL 34 55 50 -86 05 00 537 -3928 KICT WICHITA KS 37 39 17 -97 26 34 407 +03928 KICT WICHITA KS 37 39 17 -97 26 34 407 +53118 KICX CEDAR_CITY UT 37 35 27 -112 51 44 3230 13841 KILN CINCINNATI OH 39 25 13 -83 49 18 322 -4833 KILX LINCOLN IL 40 9 2 -89 20 13 177 +04833 KILX LINCOLN IL 40 9 2 -89 20 13 177 93819 KIND INDIANAPOLIS IN 39 42 27 -86 16 49 241 -3984 KINX TULSA OK 36 10 30 -95 33 53 204 -23104 KIWA PHOENIX AZ 33 17 21 -111 40 12 412 -3940 KJAN JACKSON MS 32 19 4 -90 4 48 91 +03984 KINX TULSA OK 36 10 30 -95 33 53 204 +23104 KIWA PHOENIX AZ 33 17 21 -111 40 12 412 +04844 KIWX FORT_WAYNE IN 41 21 32 -85 41 59 290 13889 KJAX JACKSONVILLE FL 30 29 5 -81 42 7 10 53824 KJGX ROBINS_AFB GA 32 40 30 -83 21 4 159 -3889 KJKL JACKSON KY 37 35 27 -83 18 47 415 -23042 KLBB LUBBOCK TX 33 39 14 -101 48 51 993 -3937 KLCH LAKE_CHARLES LA 30 7 31 -93 12 57 4 +03889 KJKL JACKSON KY 37 35 27 -83 18 47 415 +23042 KLBB LUBBOCK TX 33 39 14 -101 48 51 993 +03937 KLCH LAKE_CHARLES LA 30 7 31 -93 12 57 4 +94288 KLGX LANGLEY_HILL WA 47 7 1 -124 6 24 73 53813 KLIX NEW_ORLEANS LA 30 20 12 -89 49 32 7 94049 KLNX NORTH_PLATTE NE 41 57 28 -100 34 35 905 -4831 KLOT CHICAGO IL 41 36 17 -88 5 5 202 -4108 KLRX ELKO NV 40 44 23 -116 48 10 2056 -3982 KLSX ST_LOUIS MO 38 41 56 -90 40 58 185 +04831 KLOT CHICAGO IL 41 36 17 -88 5 5 202 +04108 KLRX ELKO NV 40 44 23 -116 48 10 2056 +03982 KLSX ST_LOUIS MO 38 41 56 -90 40 58 185 93774 KLTX WILMINGTON NC 33 59 22 -78 25 44 20 53827 KLVX LOUISVILLE KY 37 58 31 -85 56 38 219 93767 KLWX STERLING VA 38 58 31 -77 28 40 83 -3952 KLZK LITTLE_ROCK AR 34 50 12 -92 15 44 173 +03952 KLZK LITTLE_ROCK AR 34 50 12 -92 15 44 173 23023 KMAF MIDLAND/ODSSA TX 31 56 36 -102 11 21 874 -94296 KMAX MEDFORD OR 42 4 52 -122 43 2 2290 +94296 KMAX MEDFORD OR 42 4 52 -122 43 2 2290 94048 KMBX MINOT_AFB ND 48 23 33 -100 51 54 455 93768 KMHX MOREHEAD_CITY NC 34 46 34 -76 52 34 9 -4834 KMKX MILWAUKEE WI 42 58 4 -88 33 2 292 +04834 KMKX MILWAUKEE WI 42 58 4 -88 33 2 292 12838 KMLB MELBOURNE FL 28 6 48 -80 39 15 11 -13894 KMOB MOBILE AL 30 40 46 -88 14 23 63 +13894 KMOB MOBILE AL 30 40 46 -88 14 23 63 94983 KMPX MINNEAPOLIS MN 44 50 56 -93 33 56 288 94850 KMQT MARQUETTE MI 46 31 52 -87 32 54 430 53832 KMRX KNOXVILLE TN 36 10 7 -83 24 6 408 -4103 KMSX MISSOULA MT 47 2 28 -113 59 10 2394 -4104 KMTX SALT_LAKE_CTY UT 41 15 46 -112 26 52 1969 +04103 KMSX MISSOULA MT 47 2 28 -113 59 10 2394 +04104 KMTX SALT_LAKE_CTY UT 41 15 46 -112 26 52 1969 93236 KMUX SAN_FRANCISCO CA 37 9 18 -121 53 54 1057 94986 KMVX GRAND_FORKS ND 47 31 40 -97 19 32 300 53826 KMXX MAXWELL_AFB AL 32 32 12 -85 47 23 122 53115 KNKX SAN_DIEGO CA 32 55 8 -117 2 31 291 -93839 KNQA MEMPHIS TN 35 20 41 -89 52 24 86 -94980 KOAX OMAHA NE 41 19 13 -96 22 0 350 +93839 KNQA MEMPHIS TN 35 20 41 -89 52 24 86 +94980 KOAX OMAHA NE 41 19 13 -96 22 0 350 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 +04106 KOTX SPOKANE WA 47 40 49 -117 37 36 727 +03816 KPAH PADUCAH KY 37 4 6 -88 46 19 119 +04832 KPBZ PITTSBURGH PA 40 31 54 -80 13 6 361 24155 KPDT PENDLETON OR 45 41 26 -118 51 10 462 -3993 KPOE FORT_POLK LA 31 9 20 -92 58 33 124 -3021 KPUX PUEBLO CO 38 27 34 -104 10 53 1600 +03993 KPOE FORT_POLK LA 31 9 20 -92 58 33 124 +03021 KPUX PUEBLO CO 38 27 34 -104 10 53 1600 93772 KRAX RALEIGH/DUR NC 35 39 56 -78 29 23 106 -53104 KRGX RENO NV 39 45 19 -119 27 44 2530 +53104 KRGX RENO NV 39 45 19 -119 27 44 2530 24061 KRIW RIVERTON WY 43 3 58 -108 28 38 1697 53834 KRLX CHARLESTON WV 38 18 40 -81 43 23 329 -54763 KRMX GRIFFISS_AFB NY 43 28 4 -75 27 29 462 94288 KRTX PORTLAND OR 45 42 53 -122 57 56 479 -4107 KSFX POCATELLO ID 43 6 21 -112 41 10 1364 +04107 KSFX POCATELLO ID 43 6 21 -112 41 10 1364 13995 KSGF SPRINGFIELD MO 37 14 7 -93 24 2 390 13957 KSHV SHREVEPORT LA 32 27 3 -93 50 29 83 23034 KSJT SAN_ANGELO TX 31 22 17 -100 29 33 576 53117 KSOX SANTA_ANA_MT CA 33 49 4 -117 38 9 923 -92801 KTBW TAMPA FL 27 42 20 -82 24 6 12 -4102 KTFX GREAT_FALLS MT 47 27 35 -111 23 7 1132 +53906 KSRX FORT_SMITH AR 35 17 26 -94 21 42 195 +92801 KTBW TAMPA FL 27 42 20 -82 24 6 12 +04102 KTFX GREAT_FALLS MT 47 27 35 -111 23 7 1132 93805 KTLH TALLAHASSEE FL 30 23 51 -84 19 44 19 -3979 KTLX OKLAHOMA_CITY OK 35 19 59 -97 16 40 370 -3986 KTWX TOPEKA KS 38 59 49 -96 13 57 417 +03979 KTLX OKLAHOMA_CITY OK 35 19 59 -97 16 40 370 +03986 KTWX TOPEKA KS 38 59 49 -96 13 57 417 +54776 KTYX MONTAGUE NY 43 45 21 -75 40 48 563 94047 KUDX RAPID_CITY SD 44 7 30 -102 49 47 919 94981 KUEX HASTINGS NE 40 19 15 -98 26 31 602 53856 KVAX MOODY_AFB GA 30 53 25 -83 0 6 54 94234 KVBX VANDENBRG_AFB CA 34 50 17 -120 23 45 376 -3995 KVNX VANCE_AFB OK 36 44 27 -98 7 40 369 +03995 KVNX VANCE_AFB OK 36 44 27 -98 7 40 369 53101 KVTX LOS_ANGELES CA 34 24 42 -119 10 47 831 -53116 KYUX YUMA AZ 32 29 43 -114 39 24 53 +63845 KVWX EVANSVILLE IN 38 15 37 -87 43 28 58 +53116 KYUX YUMA AZ 32 29 43 -114 39 24 53 26548 PAHG ANCHORAGE AK 60 43 33 -151 21 5 74 -25404 PAIH MIDDLETON_IS AK 59 27 41 -146 18 11 20 +25404 PAIH MIDDLETON_IS AK 59 27 43 -146 18 04 20 24690 PAPD FAIRBANKS AK 65 2 6 -147 30 6 790 +00000 PABC BETHEL AK 60 47 34 -161 52 27 49 +25361 PACG SITKA AK 56 51 10 -135 31 45 82 +26641 PAEC NOME AK 64 30 41 -165 17 42 16 +25503 PAKC KING_SALMON AK 58 40 46 -156 37 46 19 41417 PGUA ANDERSEN_AFB GU 13 27 16 144 48 30 80 22548 PHKI SOUTH_KAUAI HI 21 53 39 -159 33 8 55 -43216 RKSG CMP_HUMPHRYS KOR 36 57 21 127 1 16 16 -42219 RODN KADENA OKI 26 18 7 127 54 35 66 -33771 KDGX JACKSON MS 32 16 33 -89 58 48 151 -53118 KICX CEDAR_CITY UT 37 35 9 -112 51 21 3231 -4844 KIWX N_INDIANA IN 41 21 0 -85 42 0 293 -53906 KSRX W_ARKANSAS AR 35 16 48 -94 21 0 195 -54776 KTYX MONTAGUE NY 43 45 3 -75 40 30 563 -0 KVWX EVANSVILLE IN 38 16 12 -87 43 12 190 -0 PHKM KOHALA HI 20 7 12 -155 46 12 1162 -0 PHMO MOLOKAI HI 21 7 12 -157 10 12 415 -0 PHWA HAWAII HI 19 4 48 -155 34 12 421 -0 TJUA SAN_JUAN PR 18 7 12 -66 4 48 931 +22550 PHKM KAMUELA HI 20 7 32 -155 46 40 1161 +22547 PHMO MOLOKAI HI 21 7 58 -157 10 48 415 +22513 PHWA HAWAII HI 19 5 42 -155 34 08 417 +43216 RKSG CMP_HUMPHRYS KR 36 57 34 127 1 32 16 +43219 RKJK KUNSAN KR 35 55 27 126 37 20 24 +13201 LPLA LAJES PT 38 43 49 -27 19 18 1016 +42219 RODN KADENA JP 26 20 50 127 45 58 43 +11655 TJUA SAN_JUAN PR 18 06 56 -66 04 41 852 +00000 NOP4 ROC_4 OK 35 14 10 -97 27 44 390 +00000 NOP3 ROC_3 OK 35 14 10 -97 27 44 390 +03948 KOUN NORMAN OK 35 14 10 -97 27 44 390 +03940 KJAN JACKSON MS 32 19 4 -90 4 48 91 +54763 KRMX GRIFFISS_AFB NY 43 28 4 -75 27 29 462 diff --git a/wsr88d_m31.c b/wsr88d_m31.c index 3cd4f4c..2ae2143 100644 --- a/wsr88d_m31.c +++ b/wsr88d_m31.c @@ -392,11 +392,14 @@ void wsr88d_load_ray_into_radar(Wsr88d_ray_m31 *wsr88d_ray, int isweep, int vol_index, waveform; char *type_str; - int keep_hi_prf_dz = 0; /* TODO: implement an interface for this. */ + extern int rsl_qfield[]; /* See RSL_select_fields in volume.c */ enum waveforms {surveillance=1, doppler_w_amb_res, doppler_no_amb_res, batch}; + int merging_split_cuts; + + merging_split_cuts = wsr88d_merge_split_cuts_is_set(); nfields = wsr88d_ray->ray_hdr.data_block_count - nconstblocks; field_offset = (int *) &wsr88d_ray->ray_hdr.radial_const; do_swap = little_endian(); @@ -418,6 +421,10 @@ void wsr88d_load_ray_into_radar(Wsr88d_ray_m31 *wsr88d_ray, int isweep, iray); return; } + + /* Is this field in the selected fields list? */ + if (!rsl_qfield[vol_index]) continue; + switch (vol_index) { case DZ_INDEX: f = DZ_F; invf = DZ_INVF; type_str = "Reflectivity"; break; @@ -435,13 +442,17 @@ void wsr88d_load_ray_into_radar(Wsr88d_ray_m31 *wsr88d_ray, int isweep, waveform = vcp_data.waveform[isweep]; - /* Ignore short-range reflectivity from velocity split cuts unless - * keep_hi_prf_dz is set. The indicators for this type of - * reflectivity are surveillance mode of 0 and elevation angle - * below 6 degrees. + /* If this field is reflectivity, check to see if it's from the velocity + * sweep in a split cut. If so, we normally skip it since we already + * have reflectivity from the surveillance sweep. It is kept only when + * the user has turned off merging of split cuts. We skip over this + * field if all of the following are true: surveillance PRF number is 0, + * waveform is Contiguous Doppler with Ambiguity Resolution (range + * unfolding), and we're merging split cuts. */ if (vol_index == DZ_INDEX && (vcp_data.surveil_prf_num[isweep] == 0 && - vcp_data.fixed_angle[isweep] < 6.0 && !keep_hi_prf_dz)) + vcp_data.waveform[isweep] == doppler_w_amb_res && + merging_split_cuts)) continue; /* Load the data for this field. */ @@ -543,6 +554,7 @@ Radar *wsr88d_load_m31_into_radar(Wsr88d_file *wf) msg_hdr_size = sizeof(Wsr88d_msg_hdr) - sizeof(msghdr.rpg); radar = RSL_new_radar(MAX_RADAR_VOLUMES); + memset(&wsr88d_ray, 0, sizeof(Wsr88d_ray_m31)); /* Initialize to be safe. */ while (! end_of_vos) { if (msghdr.msg_type == 31) { diff --git a/wsr88d_merge_split_cuts.c b/wsr88d_merge_split_cuts.c new file mode 100644 index 0000000..fe8f1d8 --- /dev/null +++ b/wsr88d_merge_split_cuts.c @@ -0,0 +1,129 @@ +#include +#ifndef _rsl_h +#include "rsl.h" +#endif + +static int merge_split_cuts = 1; + +void RSL_wsr88d_merge_split_cuts_on() +{ + merge_split_cuts = 1; +} + +void RSL_wsr88d_merge_split_cuts_off() +{ + merge_split_cuts = 0; +} + +void RSL_wsr88d_keep_short_refl() +{ + RSL_wsr88d_merge_split_cuts_off(); +} + +int wsr88d_merge_split_cuts_is_set() +{ + return merge_split_cuts; +} + +void wsr88d_remove_extra_refl(Radar *radar) +{ + /* This function removes any extra reflectivity for an elevation angle. + * I.e., only keep reflectivity from the surveillance sweep of a split cut. + */ + + int i; + float prev_elev; + + prev_elev = radar->v[DZ_INDEX]->sweep[0]->h.elev; + + for (i=1; i < radar->v[DZ_INDEX]->h.nsweeps; i++) { + if (radar->v[DZ_INDEX]->sweep[i]) { + if (fabsf(radar->v[DZ_INDEX]->sweep[i]->h.elev - prev_elev) < .2) { + RSL_free_sweep(radar->v[DZ_INDEX]->sweep[i]); + radar->v[DZ_INDEX]->sweep[i] = NULL; + } + else prev_elev = radar->v[DZ_INDEX]->sweep[i]->h.elev; + } + } +} + +void wsr88d_move_vcp121_extra_velsweeps(Radar *radar) +{ + /* Move the extra velocity and spectrum-width sweeps in VCP 121 split cuts + * to other volume indexes. We do this in order to allow all data moments + * at a particular split cut elevation to later be moved to the same sweep + * level in the radar structure. + */ + + int iswp; + + if (radar->v[VR_INDEX] != NULL) { + if (radar->v[V2_INDEX] == NULL) { + radar->v[V2_INDEX] = RSL_new_volume(radar->v[VR_INDEX]->h.nsweeps); + radar->v[V2_INDEX]->h.type_str = "Velocity 2"; + radar->v[V2_INDEX]->h.f = VR_F; + radar->v[V2_INDEX]->h.invf = VR_INVF; + } + if (radar->v[V3_INDEX] == NULL) { + radar->v[V3_INDEX] = RSL_new_volume(radar->v[VR_INDEX]->h.nsweeps); + radar->v[V3_INDEX]->h.type_str = "Velocity 3"; + radar->v[V3_INDEX]->h.f = VR_F; + radar->v[V3_INDEX]->h.invf = VR_INVF; + } + } + if (radar->v[SW_INDEX] != NULL) { + if (radar->v[S2_INDEX] == NULL) { + radar->v[S2_INDEX] = RSL_new_volume(radar->v[SW_INDEX]->h.nsweeps); + radar->v[S2_INDEX]->h.type_str = "Spectrum width 2"; + radar->v[S2_INDEX]->h.f = SW_F; + radar->v[S2_INDEX]->h.invf = SW_INVF; + } + if (radar->v[S3_INDEX] == NULL) { + radar->v[S3_INDEX] = RSL_new_volume(radar->v[SW_INDEX]->h.nsweeps); + radar->v[S3_INDEX]->h.type_str = "Spectrum width 3"; + radar->v[S3_INDEX]->h.f = SW_F; + radar->v[S3_INDEX]->h.invf = SW_INVF; + } + } + + for (iswp=2; iswp < 16; iswp++) { + switch (iswp) { + case 2: case 6: case 9: case 12: case 15: + if (radar->v[VR_INDEX] == NULL) break; + if (radar->v[VR_INDEX]->sweep[iswp] == NULL) break; + radar->v[V2_INDEX]->sweep[iswp] = radar->v[VR_INDEX]->sweep[iswp]; + radar->v[VR_INDEX]->sweep[iswp] = NULL; + if (radar->v[SW_INDEX] == NULL) break; + if (radar->v[SW_INDEX]->sweep[iswp] == NULL) break; + radar->v[S2_INDEX]->sweep[iswp] = radar->v[SW_INDEX]->sweep[iswp]; + radar->v[SW_INDEX]->sweep[iswp] = NULL; + break; + case 3: case 7: case 10: case 13: + if (radar->v[VR_INDEX] == NULL) break; + if (radar->v[VR_INDEX]->sweep[iswp] == NULL) break; + radar->v[V3_INDEX]->sweep[iswp] = radar->v[VR_INDEX]->sweep[iswp]; + radar->v[VR_INDEX]->sweep[iswp] = NULL; + if (radar->v[SW_INDEX] == NULL) break; + if (radar->v[SW_INDEX]->sweep[iswp] == NULL) break; + radar->v[S3_INDEX]->sweep[iswp] = radar->v[SW_INDEX]->sweep[iswp]; + radar->v[SW_INDEX]->sweep[iswp] = NULL; + break; + } + } +} + + +Radar *wsr88d_merge_split_cuts(Radar *radar) +{ + /* Move the split-cut sweeps so that all sweeps of an elevation angle + * are at the same sweep index in the Radar structure. + */ + + if (!wsr88d_merge_split_cuts_is_set()) RSL_wsr88d_merge_split_cuts_on(); + + wsr88d_remove_extra_refl(radar); + if (radar->h.vcp == 121) wsr88d_move_vcp121_extra_velsweeps(radar); + radar = RSL_prune_radar(radar); + + return radar; +} diff --git a/wsr88d_to_radar.c b/wsr88d_to_radar.c index 6966143..a290545 100644 --- a/wsr88d_to_radar.c +++ b/wsr88d_to_radar.c @@ -234,6 +234,7 @@ Radar *RSL_wsr88d_to_radar(char *infile, char *call_or_first_tape_file) char *the_file; int expected_msgtype = 0; char version[8]; + int vnum; extern int rsl_qfield[]; /* See RSL_select_fields in volume.c */ extern int *rsl_qsweep; /* See RSL_read_these_sweeps in volume.c */ @@ -295,24 +296,16 @@ Radar *RSL_wsr88d_to_radar(char *infile, char *call_or_first_tape_file) /* * Get the expected digital radar message type based on version string * from the Archive II header. The message type is 31 for Build 10, and 1 - * for prior builds. Note that we consider AR2V0001 to be message type 1, - * because it has been in the past, but with Build 10 this officially - * becomes the version number for Evansville (KVWX), which will use message - * type 31. This could be a problem if RSL is used to process KVWX. + * for prior builds. */ if (n > 0) { strncpy(version, wsr88d_file_header.title.filename, 8); - if (strncmp(version,"AR2V0007",8) == 0 || - 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; + if (strncmp(version,"AR2V",4) == 0) { + sscanf(version, "AR2V%4d", &vnum); + if (vnum > 1) expected_msgtype = 31; + else expected_msgtype = 1; } + else if (strncmp(version,"ARCHIVE2",8) == 0) expected_msgtype = 1; } if (n <= 0 || expected_msgtype == 0) { @@ -331,9 +324,7 @@ Radar *RSL_wsr88d_to_radar(char *infile, char *call_or_first_tape_file) if (expected_msgtype == 31) { - /* Get radar for message type 31. */ - nvolumes = 6; radar = wsr88d_load_m31_into_radar(wf); if (radar == NULL) return NULL; } @@ -441,6 +432,6 @@ Radar *RSL_wsr88d_to_radar(char *infile, char *call_or_first_tape_file) free(sitep); - radar = RSL_prune_radar(radar); + if (wsr88d_merge_split_cuts_is_set()) radar = wsr88d_merge_split_cuts(radar); return radar; }