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 i;
 227   int iv;
 228   int nvolumes;
 229   int volume_mask[] = {WSR88D_DZ, WSR88D_VR, WSR88D_SW};
 230   char *field_str[] = {"Reflectivity", "Velocity", "Spectrum width"};
 231   Wsr88d_site_info *sitep;
 232   char site_id_str[5];
 233   char *the_file;
 234   int expected_msgtype = 0;
 235   char version[8];
 236
 237   extern int rsl_qfield[]; /* See RSL_select_fields in volume.c */
 238   extern int *rsl_qsweep; /* See RSL_read_these_sweeps in volume.c */
 239   extern int rsl_qsweep_max;
 240
 241   Radar *wsr88d_load_m31_into_radar(Wsr88d_file *wf);
 242
 243   sitep = NULL;
 244 /* Determine the site quasi automatically.  Here is the procedure:
 245  *    1. Determine if we have a call sign.
 246  *    2. Try reading 'call_or_first_tape_file' from disk.  This is done via
 247  *       wsr88d_read_tape_header.
 248  *    3. If no valid site info, abort.
 249  */
 250   if (call_or_first_tape_file == NULL) {
 251     fprintf(stderr, "wsr88d_to_radar: No valid site ID info provided.\n");
 252     return(NULL);
 253   } else if (strlen(call_or_first_tape_file) == 4)
 254     sitep =  wsr88d_get_site(call_or_first_tape_file);
 255   else if (strlen(call_or_first_tape_file) == 0) {
 256     fprintf(stderr, "wsr88d_to_radar: No valid site ID info provided.\n");
 257     return(NULL);
 258   }
 259
 260   if (sitep == NULL)
 261     if (wsr88d_read_tape_header(call_or_first_tape_file, &wsr88d_tape_header) > 0) {
 262       memcpy(site_id_str, wsr88d_tape_header.site_id, 4);
 263       sitep  = wsr88d_get_site(site_id_str);
 264     }
 265   if (sitep == NULL) {
 266       fprintf(stderr,"wsr88d_to_radar: No valid site ID info found.\n");
 267         return(NULL);
 268   }
 269     if (radar_verbose_flag)
 270       fprintf(stderr,"SITE: %c%c%c%c\n", sitep->name[0], sitep->name[1],
 271              sitep->name[2], sitep->name[3]);
 272
 273
 274   memset(&wsr88d_sweep, 0, sizeof(Wsr88d_sweep)); /* Initialize to 0 a
 275                                                    * heavily used variable.
 276                                                    */
 277
 278 /* 1. Open the input wsr88d file. */
 279   if (infile == NULL) the_file = "stdin";  /* wsr88d.c understands this to
 280                                             * mean read from stdin.
 281                                             */
 282   else the_file = infile;
 283
 284   if ((wf = wsr88d_open(the_file)) == NULL) {
 285     wsr88d_perror(the_file);
 286     return NULL;
 287   }
 288
 289
 290
 291 /* 2. Read wsr88d headers. */
 292   /* Return # bytes, 0 or neg. on fail. */
 293   n = wsr88d_read_file_header(wf, &wsr88d_file_header);
 294   /*
 295    * Get the expected digital radar message type based on version string
 296    * from the Archive II header.  The message type is 31 for Build 10, and 1
 297    * for prior builds.  Note that we consider AR2V0001 to be message type 1,
 298    * because it has been in the past, but with Build 10 this officially
 299    * becomes the version number for Evansville (KVWX), which will use message
 300    * type 31.  This could be a problem if RSL is used to process KVWX.
 301    */
 302   if (n > 0) {
 303       strncpy(version, wsr88d_file_header.title.filename, 8);
 304       if (strncmp(version,"AR2V0006",8) == 0 ||
 305           strncmp(version,"AR2V0004",8) == 0 ||
 306           strncmp(version,"AR2V0003",8) == 0 ||
 307           strncmp(version,"AR2V0002",8) == 0) {
 308           expected_msgtype = 31;
 309       }
 310       else if (strncmp(version,"ARCHIVE2",8) == 0 ||
 311           strncmp(version,"AR2V0001",8) == 0) {
 312           expected_msgtype = 1;
 313       }
 314   }
 315
 316   if (n <= 0 || expected_msgtype == 0) {
 317       fprintf(stderr,"RSL_wsr88d_to_radar: ");
 318       if (n <= 0)
 319       fprintf(stderr,"wsr88d_read_file_header failed\n");
 320       else
 321       fprintf(stderr,"Archive II header contains unknown version "
 322           ": '%s'\n", version);
 323       wsr88d_close(wf);
 324       return NULL;
 325   }
 326
 327   if (radar_verbose_flag)
 328     print_head(wsr88d_file_header);
 329
 330
 331   if (expected_msgtype == 31) {
 332
 333       /* Get radar for message type 31. */
 334       nvolumes = 6;
 335       radar = wsr88d_load_m31_into_radar(wf);
 336       if (radar == NULL) return NULL;
 337   }
 338   else {
 339       /* Get radar for message type 1. */
 340       nvolumes = 3;
 341       /* Allocate all Volume pointers. */
 342       radar = RSL_new_radar(MAX_RADAR_VOLUMES);
 343       if (radar == NULL) return NULL;
 344
 345     /* Clear the sweep pointers. */
 346       clear_sweep(&wsr88d_sweep, 0, MAX_RAYS_IN_SWEEP);
 347
 348     /* Allocate a maximum of 30 sweeps for the volume. */
 349     /* Order is important.  WSR88D_DZ, WSR88D_VR, WSR88D_SW, is
 350      * assigned to the indexes DZ_INDEX, VR_INDEX and SW_INDEX respectively.
 351      */
 352
 353       for (iv=0; iv<nvolumes; iv++)
 354         if (rsl_qfield[iv]) radar->v[iv] = RSL_new_volume(20);
 355
 356
 357     /* LOOP until EOF */
 358       nsweep = 0;
 359       for (;(n = wsr88d_read_sweep(wf, &wsr88d_sweep)) > 0; nsweep++) {
 360         if (rsl_qsweep != NULL) {
 361           if (nsweep > rsl_qsweep_max) break;
 362           if (rsl_qsweep[nsweep] == 0) continue;
 363         }
 364         if (radar_verbose_flag)
 365         fprintf(stderr,"Processing for SWEEP # %d\n", nsweep);
 366
 367           /*  wsr88d_print_sweep_info(&wsr88d_sweep); */
 368
 369         for (iv=0; iv<nvolumes; iv++) {
 370           if (rsl_qfield[iv]) {
 371             /* Exceeded sweep limit.
 372              * Allocate more sweeps.
 373              * Copy all previous sweeps.
 374              */
 375             if (nsweep >= radar->v[iv]->h.nsweeps) {
 376               if (radar_verbose_flag)
 377                 fprintf(stderr,"Exceeded sweep allocation of %d. "
 378                     "Adding 20 more.\n", nsweep);
 379               new_volume = RSL_new_volume(radar->v[iv]->h.nsweeps+20);
 380               new_volume = copy_sweeps_into_volume(new_volume, radar->v[iv]);
 381               radar->v[iv] = new_volume;
 382             }
 383             if (wsr88d_load_sweep_into_volume(wsr88d_sweep,
 384                radar->v[iv], nsweep, volume_mask[iv]) != 0) {
 385               RSL_free_radar(radar);
 386               return NULL;
 387             }
 388           }
 389         }
 390         if (nsweep == 0) {
 391           /* Get Volume Coverage Pattern number for radar header. */
 392           i=0;
 393           while (i < MAX_RAYS_IN_SWEEP && wsr88d_sweep.ray[i] == NULL) i++;
 394           if (i < MAX_RAYS_IN_SWEEP) radar->h.vcp = wsr88d_get_volume_coverage(
 395             wsr88d_sweep.ray[i]);
 396         }
 397
 398         free_and_clear_sweep(&wsr88d_sweep, 0, MAX_RAYS_IN_SWEEP);
 399       }
 400
 401       for (iv=0; iv<nvolumes; iv++) {
 402         if (rsl_qfield[iv]) {
 403           radar->v[iv]->h.type_str = strdup(field_str[iv]);
 404           radar->v[iv]->h.nsweeps = nsweep;
 405         }
 406       }
 407   }
 408   wsr88d_close(wf);
 409
 410 /*
 411  * Here we will assign the Radar_header information.  Take most of it
 412  * from an existing volume's header.
 413  */
 414   radar_load_date_time(radar);  /* Magic :-) */
 415
 416     radar->h.number = sitep->number;
 417     memcpy(&radar->h.name, sitep->name, sizeof(sitep->name));
 418     memcpy(&radar->h.radar_name, sitep->name, sizeof(sitep->name)); /* Redundant */
 419     memcpy(&radar->h.city, sitep->city, sizeof(sitep->city));
 420     memcpy(&radar->h.state, sitep->state, sizeof(sitep->state));
 421     strcpy(radar->h.radar_type, "wsr88d");
 422     radar->h.latd = sitep->latd;
 423     radar->h.latm = sitep->latm;
 424     radar->h.lats = sitep->lats;
 425     if (radar->h.latd < 0) { /* Degree/min/sec  all the same sign */
 426       radar->h.latm *= -1;
 427       radar->h.lats *= -1;
 428     }
 429     radar->h.lond = sitep->lond;
 430     radar->h.lonm = sitep->lonm;
 431     radar->h.lons = sitep->lons;
 432     if (radar->h.lond < 0) { /* Degree/min/sec  all the same sign */
 433       radar->h.lonm *= -1;
 434       radar->h.lons *= -1;
 435     }
 436     radar->h.height = sitep->height;
 437     radar->h.spulse = sitep->spulse;
 438     radar->h.lpulse = sitep->lpulse;
 439
 440     free(sitep);
 441
 442   radar = RSL_prune_radar(radar);
 443   return radar;
 444 }