RSL v1.44

[~andy/rsl] / uf_to_radar.c

diff --git a/uf_to_radar.c b/uf_to_radar.c
index 341ec339e587de8d943ad8bf4b4edb7d50dd86f3..9835a1945a17fcf23a4ac3c9f5109794ffd042f0 100644 (file)
--- a/uf_to_radar.c
+++ b/uf_to_radar.c
@@ -24,13 +24,15 @@
 #include <string.h>
 #include <stdlib.h>
 #include <unistd.h>
+#include <math.h>
 
 /* 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
+/* Changed old buffer size (16384) for larger dualpol files.  BLK 5/18/2011 */
+typedef short UF_buffer[20000]; /* Some UF files are bigger than 4096
                                  * that the UF doc's specify.
                                  */
 
@@ -86,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; i<sweep->h.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.);
+}
 
 /********************************************************************/
 /*********************************************************************/
@@ -149,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;
@@ -335,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;
@@ -391,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) {
@@ -485,7 +556,7 @@ Radar *RSL_uf_to_radar_fp(FILE *fp)
 
 
   radar = NULL;
-  setvbuf(fp,NULL,_IOFBF,(size_t)NEW_BUFSIZ); /* Faster i/o? */
+  /* setvbuf(fp,NULL,_IOFBF,(size_t)NEW_BUFSIZ); * Faster i/o? */
   if (fread(magic.buf, sizeof(char), 6, fp) <= 0) return NULL;
 /*
  * Check for fortran record length delimeters, NCAR kludge.
@@ -567,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;