Pileus Git - ~andy/rsl/blob - wsr88d_to_radar.c

   1 /*
   2     NASA/TRMM, Code 910.1.
   3     This is the TRMM Office Radar Software Library.
   4     Copyright (C) 1996, 1997
   5             John H. Merritt
   6             Space Applications Corporation
   7             Vienna, Virginia
   8
   9     This library is free software; you can redistribute it and/or
  10     modify it under the terms of the GNU Library General Public
  11     License as published by the Free Software Foundation; either
  12     version 2 of the License, or (at your option) any later version.
  13
  14     This library is distributed in the hope that it will be useful,
  15     but WITHOUT ANY WARRANTY; without even the implied warranty of
  16     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17     Library General Public License for more details.
  18
  19     You should have received a copy of the GNU Library General Public
  20     License along with this library; if not, write to the Free
  21     Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  22 */
  23 #include <stdio.h>
  24 #include <string.h>
  25
  26 #include <stdlib.h>
  27
  28 #include "rsl.h"
  29 #include "wsr88d.h"
  30
  31 extern int radar_verbose_flag;
  32 /*
  33  * These externals can be found in the wsr88d library; secret code.
  34  */
  35 void print_head(Wsr88d_file_header wsr88d_file_header);
  36 void clear_sweep(Wsr88d_sweep *wsr88d_sweep, int x, int n);
  37 void free_and_clear_sweep(Wsr88d_sweep *wsr88d_sweep, int x, int n);
  38 /*
  39  * Secretly in uf_to_radar.c
  40  */
  41 Volume *copy_sweeps_into_volume(Volume *new_volume, Volume *old_volume);
  42
  43
  44 void float_to_range(float *x, Range *c, int n, Range (*function)(float x) )
  45 {
  46   while (n--) {
  47         if (*x == WSR88D_BADVAL) *c = function(BADVAL);
  48         else if (*x == WSR88D_RFVAL) *c = function(RFVAL);
  49         else *c = function(*x);
  50         c++; x++;
  51   }
  52 }
  53
  54 /**********************************************************************/
  55 /*                                                                    */
  56 /* done 3/23         wsr88d_load_sweep_into_volume                    */
  57 /*                                                                    */
  58 /*  By: John Merritt                                                  */
  59 /*      Space Applications Corporation                                */
  60 /*      March 3, 1994                                                 */
  61 /**********************************************************************/
  62 int wsr88d_load_sweep_into_volume(Wsr88d_sweep ws,
  63                                            Volume *v, int nsweep, unsigned int vmask)
  64 {
  65   int i;
  66   int iray;
  67   float v_data[1000];
  68   Range  c_data[1000];
  69   int n;
  70
  71   int mon, day, year;
  72   int hh, mm, ss;
  73   float fsec;
  74   int vol_cpat;
  75
  76   Ray *ray_ptr;
  77   Range (*invf)(float x);
  78   float (*f)(Range x);
  79
  80   /* Allocate memory for MAX_RAYS_IN_SWEEP rays. */
  81   v->sweep[nsweep] = RSL_new_sweep(MAX_RAYS_IN_SWEEP);
  82   if (v->sweep[nsweep] == NULL) {
  83         perror("wsr88d_load_sweep_into_volume: RSL_new_sweep");
  84         return -1;
  85   }
  86
  87   v->sweep[nsweep]->h.nrays = 0;
  88   f = (float (*)(Range x))NULL;
  89   invf = (Range (*)(float x))NULL;
  90   if (vmask & WSR88D_DZ) { invf = DZ_INVF; f = DZ_F; }
  91   if (vmask & WSR88D_VR) { invf = VR_INVF; f = VR_F; }
  92   if (vmask & WSR88D_SW) { invf = SW_INVF; f = SW_F; }
  93
  94   v->h.invf = invf;
  95   v->h.f    = f;
  96   v->sweep[nsweep]->h.invf = invf;
  97   v->sweep[nsweep]->h.f    = f;
  98
  99   for (i=0,iray=0; i<MAX_RAYS_IN_SWEEP; i++) {
 100         if (ws.ray[i] != NULL) {
 101           wsr88d_ray_to_float(ws.ray[i], vmask, v_data, &n);
 102           float_to_range(v_data, c_data, n, invf);
 103           if (n > 0) {
 104                 wsr88d_get_date(ws.ray[i], &mon, &day, &year);
 105                 wsr88d_get_time(ws.ray[i], &hh, &mm, &ss, &fsec);
 106                 /*
 107           fprintf(stderr,"n %d, mon %d, day %d, year %d,  hour %d, min %d, sec %d, fsec %f\n",
 108                         n, mon, day, year, hh, mm, ss, fsec);
 109                         */
 110                 /*
 111                  * Load the sweep/ray headar information.
 112                  */
 113
 114                 v->sweep[nsweep]->ray[iray] = RSL_new_ray(n);
 115                 /*(Range *)calloc(n, sizeof(Range)); */
 116
 117                 ray_ptr = v->sweep[nsweep]->ray[iray]; /* Make code below readable. */
 118                 ray_ptr->h.f        = f;
 119                 ray_ptr->h.invf     = invf;
 120                 ray_ptr->h.month    = mon;
 121                 ray_ptr->h.day      = day;
 122                 ray_ptr->h.year     = year + 1900; /* Yes 1900 makes this year 2000 compliant, due to wsr88d using unix time(). */
 123                 ray_ptr->h.hour     = hh;
 124                 ray_ptr->h.minute   = mm;
 125                 ray_ptr->h.sec      = ss + fsec;
 126                 ray_ptr->h.unam_rng = wsr88d_get_range   (ws.ray[i]);
 127                 ray_ptr->h.azimuth  = wsr88d_get_azimuth (ws.ray[i]);
 128 /* -180 to +180 is converted to 0 to 360 */
 129                 if (ray_ptr->h.azimuth < 0) ray_ptr->h.azimuth += 360;
 130                 ray_ptr->h.ray_num  = ws.ray[i]->ray_num;
 131                 ray_ptr->h.elev       = wsr88d_get_elevation_angle(ws.ray[i]);
 132                 ray_ptr->h.elev_num   = ws.ray[i]->elev_num;
 133                 if (vmask & WSR88D_DZ) {
 134                   ray_ptr->h.range_bin1 = ws.ray[i]->refl_rng;
 135                   ray_ptr->h.gate_size  = ws.ray[i]->refl_size;
 136                 } else {
 137                   ray_ptr->h.range_bin1 = ws.ray[i]->dop_rng;
 138                   ray_ptr->h.gate_size  = ws.ray[i]->dop_size;
 139                 }
 140                 ray_ptr->h.vel_res  = wsr88d_get_velocity_resolution(ws.ray[i]);
 141                 vol_cpat = wsr88d_get_volume_coverage(ws.ray[i]);
 142                 switch (vol_cpat) {
 143                 case 11: ray_ptr->h.sweep_rate = 16.0/5.0;  break;
 144                 case 12: ray_ptr->h.sweep_rate = 17.0/4.2;  break;
 145                 case 21: ray_ptr->h.sweep_rate = 11.0/6.0;  break;
 146                 case 31: ray_ptr->h.sweep_rate =  8.0/10.0; break;
 147                 case 32: ray_ptr->h.sweep_rate =  7.0/10.0; break;
 148                 case 121:ray_ptr->h.sweep_rate = 20.0/5.5;  break;
 149                 default: ray_ptr->h.sweep_rate =  0.0; break;
 150                 }
 151
 152                 ray_ptr->h.nyq_vel  = wsr88d_get_nyquist(ws.ray[i]);
 153                 ray_ptr->h.azim_rate   = wsr88d_get_azimuth_rate(ws.ray[i]);
 154                 ray_ptr->h.fix_angle   = wsr88d_get_fix_angle(ws.ray[i]);
 155                 ray_ptr->h.pulse_count = wsr88d_get_pulse_count(ws.ray[i]);
 156                 ray_ptr->h.pulse_width = wsr88d_get_pulse_width(ws.ray[i]);
 157                 ray_ptr->h.beam_width  = .95;
 158                 ray_ptr->h.prf         = wsr88d_get_prf(ws.ray[i]);
 159                 ray_ptr->h.frequency   = wsr88d_get_frequency(ws.ray[i]);
 160                 ray_ptr->h.wavelength = 0.1071; /* Previously called
 161                                              * wsr88d_get_wavelength(ws.ray[i]).
 162                                              * See wsr88d.c for explanation.
 163                                              */
 164
 165                 /* It is no coincidence that the 'vmask' and wsr88d datatype
 166                  * values are the same.  We expect 'vmask' to be one of
 167                  * REFL_MASK, VEL_MASK, or SW_MASK.  These match WSR88D_DZ,
 168                  * WSR88D_VR, and WSR88D_SW in the wsr88d library.
 169                  */
 170                 ray_ptr->h.nbins = n;
 171                 memcpy(ray_ptr->range, c_data, n*sizeof(Range));
 172                 v->sweep[nsweep]->h.nrays = iray+1;
 173                 v->sweep[nsweep]->h.elev += ray_ptr->h.elev;
 174                 v->sweep[nsweep]->h.sweep_num = ray_ptr->h.elev_num;
 175                 iray++;
 176           }
 177         }
 178   }
 179   v->sweep[nsweep]->h.beam_width = .95;
 180   v->sweep[nsweep]->h.vert_half_bw = .475;
 181   v->sweep[nsweep]->h.horz_half_bw = .475;
 182   /* Now calculate the mean elevation angle for this sweep. */
 183   if (v->sweep[nsweep]->h.nrays > 0)
 184         v->sweep[nsweep]->h.elev /= v->sweep[nsweep]->h.nrays;
 185   else {
 186         free(v->sweep[nsweep]);  /* No rays loaded, free this sweep. */
 187         v->sweep[nsweep] = NULL;
 188   }
 189
 190   return 0;
 191 }
 192
 193 /**********************************************************************/
 194 /*                                                                    */
 195 /*                     RSL_wsr88d_to_radar                            */
 196 /*                                                                    */
 197 /*  By: John Merritt                                                  */
 198 /*      Space Applications Corporation                                */
 199 /*      March 3, 1994                                                 */
 200 /**********************************************************************/
 201
 202 Radar *RSL_wsr88d_to_radar(char *infile, char *call_or_first_tape_file)
 203 /*
 204  * Gets all volumes from the nexrad file.  Input file is 'infile'.
 205  * Site information is extracted from 'call_or_first_tape_file'; this
 206  * is typically a disk file called 'nex.file.1'.
 207  *
 208  *  -or-
 209  *
 210  * Uses the string in 'call_or_first_tape_file' as the 4 character call sign
 211  * for the sight.  All UPPERCASE characters.  Normally, this call sign
 212  * is extracted from the file 'nex.file.1'.
 213  *
 214  * Returns a pointer to a Radar structure; that contains the different
 215  * Volumes of data.
 216  */
 217 {
 218   Radar *radar;
 219   Volume *new_volume;
 220   Wsr88d_file *wf;
 221   Wsr88d_sweep wsr88d_sweep;
 222   Wsr88d_file_header wsr88d_file_header;
 223   Wsr88d_tape_header wsr88d_tape_header;
 224   int n;
 225   int nsweep;
 226   int iv;
 227   int nvolumes;
 228   int volume_mask[] = {WSR88D_DZ, WSR88D_VR, WSR88D_SW};
 229   char *field_str[] = {"Reflectivity", "Velocity", "Spectrum width"};
 230   Wsr88d_site_info *sitep;
 231   char site_id_str[5];
 232   char *the_file;
 233   int expected_msgtype = 0;
 234   char version[8];
 235
 236   extern int rsl_qfield[]; /* See RSL_select_fields in volume.c */
 237   extern int *rsl_qsweep; /* See RSL_read_these_sweeps in volume.c */
 238   extern int rsl_qsweep_max;
 239
 240   Radar *load_wsr88d_m31_into_radar(Wsr88d_file *wf);
 241
 242   sitep = NULL;
 243 /* Determine the site quasi automatically.  Here is the procedure:
 244  *    1. Determine if we have a call sign.
 245  *    2. Try reading 'call_or_first_tape_file' from disk.  This is done via
 246  *       wsr88d_read_tape_header.
 247  *    3. If no valid site info, abort.
 248  */
 249   if (call_or_first_tape_file == NULL) {
 250         fprintf(stderr, "wsr88d_to_radar: No valid site ID info provided.\n");
 251         return(NULL);
 252   }     else if (strlen(call_or_first_tape_file) == 4)
 253         sitep =  wsr88d_get_site(call_or_first_tape_file);
 254   else if (strlen(call_or_first_tape_file) == 0) {
 255         fprintf(stderr, "wsr88d_to_radar: No valid site ID info provided.\n");
 256         return(NULL);
 257   }
 258
 259   if (sitep == NULL)
 260         if (wsr88d_read_tape_header(call_or_first_tape_file, &wsr88d_tape_header) > 0) {
 261           memcpy(site_id_str, wsr88d_tape_header.site_id, 4);
 262           sitep  = wsr88d_get_site(site_id_str);
 263         }
 264   if (sitep == NULL) {
 265           fprintf(stderr,"wsr88d_to_radar: No valid site ID info found.\n");
 266                 return(NULL);
 267   }
 268         if (radar_verbose_flag)
 269           fprintf(stderr,"SITE: %c%c%c%c\n", sitep->name[0], sitep->name[1],
 270                          sitep->name[2], sitep->name[3]);
 271
 272
 273   memset(&wsr88d_sweep, 0, sizeof(Wsr88d_sweep)); /* Initialize to 0 a
 274                                                                                                    * heavily used variable.
 275                                                                                                    */
 276
 277 /* 1. Open the input wsr88d file. */
 278   if (infile == NULL) the_file = "stdin";  /* wsr88d.c understands this to
 279                                                                                         * mean read from stdin.
 280                                                                                         */
 281   else the_file = infile;
 282
 283   if ((wf = wsr88d_open(the_file)) == NULL) {
 284         wsr88d_perror(the_file);
 285         return NULL;
 286   }
 287
 288
 289
 290 /* 2. Read wsr88d headers. */
 291   /* Return # bytes, 0 or neg. on fail. */
 292   n = wsr88d_read_file_header(wf, &wsr88d_file_header);
 293   /*
 294    * Get the expected digital radar message type based on version string
 295    * from the Archive II header.  The message type is 31 for Build 10, and 1
 296    * for prior builds.  Note that we consider AR2V0001 to be message type 1,
 297    * because it has been in the past, but with Build 10 this officially
 298    * becomes the version number for Evansville (KVWX), which will use message
 299    * type 31.  This could be a problem if RSL is used to process KVWX.
 300    */
 301   if (n > 0) {
 302       strncpy(version, wsr88d_file_header.title.filename, 8);
 303       if (strncmp(version,"AR2V0004",8) == 0 ||
 304               strncmp(version,"AR2V0003",8) ==0 ||
 305               strncmp(version,"AR2V0002",8) == 0) {
 306           expected_msgtype = 31;
 307       }
 308       else if (strncmp(version,"ARCHIVE2",8) == 0 ||
 309               strncmp(version,"AR2V0001",8) == 0) {
 310           expected_msgtype = 1;
 311       }
 312   }
 313
 314   if (n <= 0 || expected_msgtype == 0) {
 315       fprintf(stderr,"RSL_wsr88d_to_radar: ");
 316       if (n <= 0)
 317           fprintf(stderr,"wsr88d_read_file_header failed\n");
 318       else
 319           fprintf(stderr,"Archive II header contains unknown version "
 320                   ": '%s'\n", version);
 321       wsr88d_close(wf);
 322       free(radar);
 323       return NULL;
 324   }
 325
 326   if (radar_verbose_flag)
 327         print_head(wsr88d_file_header);
 328
 329
 330   if (expected_msgtype == 31) {
 331
 332       /* Get radar for message type 31. */
 333       nvolumes = 6;
 334       radar = load_wsr88d_m31_into_radar(wf);
 335       if (radar == NULL) return NULL;
 336   }
 337   else {
 338       /* Get radar for message type 1. */
 339       nvolumes = 3;
 340       /* Allocate all Volume pointers. */
 341       radar = RSL_new_radar(MAX_RADAR_VOLUMES);
 342       if (radar == NULL) return NULL;
 343
 344     /* Clear the sweep pointers. */
 345       clear_sweep(&wsr88d_sweep, 0, MAX_RAYS_IN_SWEEP);
 346
 347     /* Allocate a maximum of 30 sweeps for the volume. */
 348     /* Order is important.  WSR88D_DZ, WSR88D_VR, WSR88D_SW, is
 349      * assigned to the indexes DZ_INDEX, VR_INDEX and SW_INDEX respectively.
 350      */
 351
 352       for (iv=0; iv<nvolumes; iv++)
 353             if (rsl_qfield[iv]) radar->v[iv] = RSL_new_volume(20);
 354
 355
 356     /* LOOP until EOF */
 357       nsweep = 0;
 358       for (;(n = wsr88d_read_sweep(wf, &wsr88d_sweep)) > 0; nsweep++) {
 359             if (rsl_qsweep != NULL) {
 360               if (nsweep > rsl_qsweep_max) break;
 361               if (rsl_qsweep[nsweep] == 0) continue;
 362             }
 363             if (radar_verbose_flag)
 364             fprintf(stderr,"Processing for SWEEP # %d\n", nsweep);
 365
 366               /*  wsr88d_print_sweep_info(&wsr88d_sweep); */
 367
 368             for (iv=0; iv<nvolumes; iv++) {
 369               if (rsl_qfield[iv]) {
 370                     /* Exceeded sweep limit.
 371                      * Allocate more sweeps.
 372                      * Copy all previous sweeps.
 373                      */
 374                     if (nsweep >= radar->v[iv]->h.nsweeps) {
 375                       if (radar_verbose_flag)
 376                             fprintf(stderr,"Exceeded sweep allocation of %d. "
 377                                     "Adding 20 more.\n", nsweep);
 378                       new_volume = RSL_new_volume(radar->v[iv]->h.nsweeps+20);
 379                       new_volume = copy_sweeps_into_volume(new_volume, radar->v[iv]);
 380                       radar->v[iv] = new_volume;
 381                     }
 382                     if (wsr88d_load_sweep_into_volume(wsr88d_sweep,
 383                            radar->v[iv], nsweep, volume_mask[iv]) != 0) {
 384                       RSL_free_radar(radar);
 385                       return NULL;
 386                     }
 387               }
 388             }
 389             free_and_clear_sweep(&wsr88d_sweep, 0, MAX_RAYS_IN_SWEEP);
 390       }
 391
 392       for (iv=0; iv<nvolumes; iv++) {
 393             if (rsl_qfield[iv]) {
 394               radar->v[iv]->h.type_str = strdup(field_str[iv]);
 395               radar->v[iv]->h.nsweeps = nsweep;
 396             }
 397       }
 398   }
 399   wsr88d_close(wf);
 400
 401 /*
 402  * Here we will assign the Radar_header information.  Take most of it
 403  * from an existing volume's header.
 404  */
 405   radar_load_date_time(radar);  /* Magic :-) */
 406
 407         radar->h.number = sitep->number;
 408         memcpy(&radar->h.name, sitep->name, sizeof(sitep->name));
 409         memcpy(&radar->h.radar_name, sitep->name, sizeof(sitep->name)); /* Redundant */
 410         memcpy(&radar->h.city, sitep->city, sizeof(sitep->city));
 411         memcpy(&radar->h.state, sitep->state, sizeof(sitep->state));
 412         strcpy(radar->h.radar_type, "wsr88d");
 413         radar->h.latd = sitep->latd;
 414         radar->h.latm = sitep->latm;
 415         radar->h.lats = sitep->lats;
 416         if (radar->h.latd < 0) { /* Degree/min/sec  all the same sign */
 417           radar->h.latm *= -1;
 418           radar->h.lats *= -1;
 419         }
 420         radar->h.lond = sitep->lond;
 421         radar->h.lonm = sitep->lonm;
 422         radar->h.lons = sitep->lons;
 423         if (radar->h.lond < 0) { /* Degree/min/sec  all the same sign */
 424           radar->h.lonm *= -1;
 425           radar->h.lons *= -1;
 426         }
 427         radar->h.height = sitep->height;
 428         radar->h.spulse = sitep->spulse;
 429         radar->h.lpulse = sitep->lpulse;
 430
 431   radar = RSL_prune_radar(radar);
 432   return radar;
 433 }