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.elev = 0;
  88   v->sweep[nsweep]->h.nrays = 0;
  89   f = (float (*)(Range x))NULL;
  90   invf = (Range (*)(float x))NULL;
  91   if (vmask & WSR88D_DZ) { invf = DZ_INVF; f = DZ_F; }
  92   if (vmask & WSR88D_VR) { invf = VR_INVF; f = VR_F; }
  93   if (vmask & WSR88D_SW) { invf = SW_INVF; f = SW_F; }
  94
  95   v->h.invf = invf;
  96   v->h.f    = f;
  97   v->sweep[nsweep]->h.invf = invf;
  98   v->sweep[nsweep]->h.f    = f;
  99
 100   for (i=0,iray=0; i<MAX_RAYS_IN_SWEEP; i++) {
 101     if (ws.ray[i] != NULL) {
 102       wsr88d_ray_to_float(ws.ray[i], vmask, v_data, &n);
 103       float_to_range(v_data, c_data, n, invf);
 104       if (n > 0) {
 105         wsr88d_get_date(ws.ray[i], &mon, &day, &year);
 106         wsr88d_get_time(ws.ray[i], &hh, &mm, &ss, &fsec);
 107         /*
 108       fprintf(stderr,"n %d, mon %d, day %d, year %d,  hour %d, min %d, sec %d, fsec %f\n",
 109             n, mon, day, year, hh, mm, ss, fsec);
 110             */
 111         /*
 112          * Load the sweep/ray headar information.
 113          */
 114
 115         v->sweep[nsweep]->ray[iray] = RSL_new_ray(n);
 116         /*(Range *)calloc(n, sizeof(Range)); */
 117
 118         ray_ptr = v->sweep[nsweep]->ray[iray]; /* Make code below readable. */
 119         ray_ptr->h.f        = f;
 120         ray_ptr->h.invf     = invf;
 121         ray_ptr->h.month    = mon;
 122         ray_ptr->h.day      = day;
 123         ray_ptr->h.year     = year + 1900; /* Yes 1900 makes this year 2000 compliant, due to wsr88d using unix time(). */
 124         ray_ptr->h.hour     = hh;
 125         ray_ptr->h.minute   = mm;
 126         ray_ptr->h.sec      = ss + fsec;
 127         ray_ptr->h.unam_rng = wsr88d_get_range   (ws.ray[i]);
 128         ray_ptr->h.azimuth  = wsr88d_get_azimuth (ws.ray[i]);
 129 /* -180 to +180 is converted to 0 to 360 */
 130         if (ray_ptr->h.azimuth < 0) ray_ptr->h.azimuth += 360;
 131         ray_ptr->h.ray_num  = ws.ray[i]->ray_num;
 132         ray_ptr->h.elev       = wsr88d_get_elevation_angle(ws.ray[i]);
 133         ray_ptr->h.elev_num   = ws.ray[i]->elev_num;
 134         if (vmask & WSR88D_DZ) {
 135           ray_ptr->h.range_bin1 = ws.ray[i]->refl_rng;
 136           ray_ptr->h.gate_size  = ws.ray[i]->refl_size;
 137         } else {
 138           ray_ptr->h.range_bin1 = ws.ray[i]->dop_rng;
 139           ray_ptr->h.gate_size  = ws.ray[i]->dop_size;
 140         }
 141         ray_ptr->h.vel_res  = wsr88d_get_velocity_resolution(ws.ray[i]);
 142         vol_cpat = wsr88d_get_volume_coverage(ws.ray[i]);
 143         switch (vol_cpat) {
 144         case 11: ray_ptr->h.sweep_rate = 16.0/5.0;  break;
 145         case 12: ray_ptr->h.sweep_rate = 17.0/4.2;  break;
 146         case 21: ray_ptr->h.sweep_rate = 11.0/6.0;  break;
 147         case 31: ray_ptr->h.sweep_rate =  8.0/10.0; break;
 148         case 32: ray_ptr->h.sweep_rate =  7.0/10.0; break;
 149         case 121:ray_ptr->h.sweep_rate = 20.0/5.5;  break;
 150         default: ray_ptr->h.sweep_rate =  0.0; break;
 151         }
 152
 153         ray_ptr->h.nyq_vel  = wsr88d_get_nyquist(ws.ray[i]);
 154         ray_ptr->h.azim_rate   = wsr88d_get_azimuth_rate(ws.ray[i]);
 155         ray_ptr->h.fix_angle   = wsr88d_get_fix_angle(ws.ray[i]);
 156         ray_ptr->h.pulse_count = wsr88d_get_pulse_count(ws.ray[i]);
 157         ray_ptr->h.pulse_width = wsr88d_get_pulse_width(ws.ray[i]);
 158         ray_ptr->h.beam_width  = .95;
 159         ray_ptr->h.prf         = wsr88d_get_prf(ws.ray[i]);
 160         ray_ptr->h.frequency   = wsr88d_get_frequency(ws.ray[i]);
 161         ray_ptr->h.wavelength = 0.1071; /* Previously called
 162                                      * wsr88d_get_wavelength(ws.ray[i]).
 163                                      * See wsr88d.c for explanation.
 164                          */
 165
 166         /* It is no coincidence that the 'vmask' and wsr88d datatype
 167          * values are the same.  We expect 'vmask' to be one of
 168          * REFL_MASK, VEL_MASK, or SW_MASK.  These match WSR88D_DZ,
 169          * WSR88D_VR, and WSR88D_SW in the wsr88d library.
 170          */
 171         ray_ptr->h.nbins = n;
 172         memcpy(ray_ptr->range, c_data, n*sizeof(Range));
 173         v->sweep[nsweep]->h.nrays = iray+1;
 174         v->sweep[nsweep]->h.elev += ray_ptr->h.elev;
 175         v->sweep[nsweep]->h.sweep_num = ray_ptr->h.elev_num;
 176         iray++;
 177       }
 178     }
 179   }
 180   v->sweep[nsweep]->h.beam_width = .95;
 181   v->sweep[nsweep]->h.vert_half_bw = .475;
 182   v->sweep[nsweep]->h.horz_half_bw = .475;
 183   /* Now calculate the mean elevation angle for this sweep. */
 184   if (v->sweep[nsweep]->h.nrays > 0)
 185     v->sweep[nsweep]->h.elev /= v->sweep[nsweep]->h.nrays;
 186   else {
 187     RSL_free_sweep(v->sweep[nsweep]); /* No rays loaded, free this sweep. */
 188     v->sweep[nsweep] = NULL;
 189   }
 190
 191   return 0;
 192 }
 193
 194 /**********************************************************************/
 195 /*                                                                    */
 196 /*                     RSL_wsr88d_to_radar                            */
 197 /*                                                                    */
 198 /*  By: John Merritt                                                  */
 199 /*      Space Applications Corporation                                */
 200 /*      March 3, 1994                                                 */
 201 /**********************************************************************/
 202
 203 Radar *RSL_wsr88d_to_radar(char *infile, char *call_or_first_tape_file)
 204 /*
 205  * Gets all volumes from the nexrad file.  Input file is 'infile'.
 206  * Site information is extracted from 'call_or_first_tape_file'; this
 207  * is typically a disk file called 'nex.file.1'.
 208  *
 209  *  -or-
 210  *
 211  * Uses the string in 'call_or_first_tape_file' as the 4 character call sign
 212  * for the sight.  All UPPERCASE characters.  Normally, this call sign
 213  * is extracted from the file 'nex.file.1'.
 214  *
 215  * Returns a pointer to a Radar structure; that contains the different
 216  * Volumes of data.
 217  */
 218 {
 219   Radar *radar;
 220   Volume *new_volume;
 221   Wsr88d_file *wf;
 222   Wsr88d_sweep wsr88d_sweep;
 223   Wsr88d_file_header wsr88d_file_header;
 224   Wsr88d_tape_header wsr88d_tape_header;
 225   int n;
 226   int nsweep;
 227   int i;
 228   int iv;
 229   int nvolumes;
 230   int volume_mask[] = {WSR88D_DZ, WSR88D_VR, WSR88D_SW};
 231   char *field_str[] = {"Reflectivity", "Velocity", "Spectrum width"};
 232   Wsr88d_site_info *sitep;
 233   char site_id_str[5];
 234   char *the_file;
 235   int expected_msgtype = 0;
 236   char version[8];
 237   int vnum;
 238
 239   extern int rsl_qfield[]; /* See RSL_select_fields in volume.c */
 240   extern int *rsl_qsweep; /* See RSL_read_these_sweeps in volume.c */
 241   extern int rsl_qsweep_max;
 242
 243   Radar *wsr88d_load_m31_into_radar(Wsr88d_file *wf);
 244
 245   sitep = NULL;
 246 /* Determine the site quasi automatically.  Here is the procedure:
 247  *    1. Determine if we have a call sign.
 248  *    2. Try reading 'call_or_first_tape_file' from disk.  This is done via
 249  *       wsr88d_read_tape_header.
 250  *    3. If no valid site info, abort.
 251  */
 252   if (call_or_first_tape_file == NULL) {
 253     fprintf(stderr, "wsr88d_to_radar: No valid site ID info provided.\n");
 254     return(NULL);
 255   } else if (strlen(call_or_first_tape_file) == 4)
 256     sitep =  wsr88d_get_site(call_or_first_tape_file);
 257   else if (strlen(call_or_first_tape_file) == 0) {
 258     fprintf(stderr, "wsr88d_to_radar: No valid site ID info provided.\n");
 259     return(NULL);
 260   }
 261
 262   if (sitep == NULL)
 263     if (wsr88d_read_tape_header(call_or_first_tape_file, &wsr88d_tape_header) > 0) {
 264       memcpy(site_id_str, wsr88d_tape_header.site_id, 4);
 265       sitep  = wsr88d_get_site(site_id_str);
 266     }
 267   if (sitep == NULL) {
 268       fprintf(stderr,"wsr88d_to_radar: No valid site ID info found.\n");
 269         return(NULL);
 270   }
 271     if (radar_verbose_flag)
 272       fprintf(stderr,"SITE: %c%c%c%c\n", sitep->name[0], sitep->name[1],
 273              sitep->name[2], sitep->name[3]);
 274
 275
 276   memset(&wsr88d_sweep, 0, sizeof(Wsr88d_sweep)); /* Initialize to 0 a
 277                                                    * heavily used variable.
 278                                                    */
 279
 280 /* 1. Open the input wsr88d file. */
 281   if (infile == NULL) the_file = "stdin";  /* wsr88d.c understands this to
 282                                             * mean read from stdin.
 283                                             */
 284   else the_file = infile;
 285
 286   if ((wf = wsr88d_open(the_file)) == NULL) {
 287     wsr88d_perror(the_file);
 288     return NULL;
 289   }
 290
 291
 292
 293 /* 2. Read wsr88d headers. */
 294   /* Return # bytes, 0 or neg. on fail. */
 295   n = wsr88d_read_file_header(wf, &wsr88d_file_header);
 296   /*
 297    * Get the expected digital radar message type based on version string
 298    * from the Archive II header.  The message type is 31 for Build 10, and 1
 299    * for prior builds.
 300    */
 301   if (n > 0) {
 302       strncpy(version, wsr88d_file_header.title.filename, 8);
 303       if (strncmp(version,"AR2V",4) == 0) {
 304           sscanf(version, "AR2V%4d", &vnum);
 305           if (vnum > 1) expected_msgtype = 31;
 306           else expected_msgtype = 1;
 307       }
 308       else if (strncmp(version,"ARCHIVE2",8) == 0) expected_msgtype = 1;
 309   }
 310
 311   if (n <= 0 || expected_msgtype == 0) {
 312       fprintf(stderr,"RSL_wsr88d_to_radar: ");
 313       if (n <= 0)
 314       fprintf(stderr,"wsr88d_read_file_header failed\n");
 315       else
 316       fprintf(stderr,"Archive II header contains unknown version "
 317           ": '%s'\n", version);
 318       wsr88d_close(wf);
 319       return NULL;
 320   }
 321
 322   if (radar_verbose_flag)
 323     print_head(wsr88d_file_header);
 324
 325
 326   if (expected_msgtype == 31) {
 327       /* Get radar for message type 31. */
 328       radar = wsr88d_load_m31_into_radar(wf);
 329       if (radar == NULL) return NULL;
 330   }
 331   else {
 332       /* Get radar for message type 1. */
 333       nvolumes = 3;
 334       /* Allocate all Volume pointers. */
 335       radar = RSL_new_radar(MAX_RADAR_VOLUMES);
 336       if (radar == NULL) return NULL;
 337
 338     /* Clear the sweep pointers. */
 339       clear_sweep(&wsr88d_sweep, 0, MAX_RAYS_IN_SWEEP);
 340
 341     /* Allocate a maximum of 30 sweeps for the volume. */
 342     /* Order is important.  WSR88D_DZ, WSR88D_VR, WSR88D_SW, is
 343      * assigned to the indexes DZ_INDEX, VR_INDEX and SW_INDEX respectively.
 344      */
 345
 346       for (iv=0; iv<nvolumes; iv++)
 347         if (rsl_qfield[iv]) radar->v[iv] = RSL_new_volume(20);
 348
 349
 350     /* LOOP until EOF */
 351       nsweep = 0;
 352       for (;(n = wsr88d_read_sweep(wf, &wsr88d_sweep)) > 0; nsweep++) {
 353         if (rsl_qsweep != NULL) {
 354           if (nsweep > rsl_qsweep_max) break;
 355           if (rsl_qsweep[nsweep] == 0) continue;
 356         }
 357         if (radar_verbose_flag)
 358         fprintf(stderr,"Processing for SWEEP # %d\n", nsweep);
 359
 360           /*  wsr88d_print_sweep_info(&wsr88d_sweep); */
 361
 362         for (iv=0; iv<nvolumes; iv++) {
 363           if (rsl_qfield[iv]) {
 364             /* Exceeded sweep limit.
 365              * Allocate more sweeps.
 366              * Copy all previous sweeps.
 367              */
 368             if (nsweep >= radar->v[iv]->h.nsweeps) {
 369               if (radar_verbose_flag)
 370                 fprintf(stderr,"Exceeded sweep allocation of %d. "
 371                     "Adding 20 more.\n", nsweep);
 372               new_volume = RSL_new_volume(radar->v[iv]->h.nsweeps+20);
 373               new_volume = copy_sweeps_into_volume(new_volume, radar->v[iv]);
 374               radar->v[iv] = new_volume;
 375             }
 376             if (wsr88d_load_sweep_into_volume(wsr88d_sweep,
 377                radar->v[iv], nsweep, volume_mask[iv]) != 0) {
 378               RSL_free_radar(radar);
 379               return NULL;
 380             }
 381           }
 382         }
 383         if (nsweep == 0) {
 384           /* Get Volume Coverage Pattern number for radar header. */
 385           i=0;
 386           while (i < MAX_RAYS_IN_SWEEP && wsr88d_sweep.ray[i] == NULL) i++;
 387           if (i < MAX_RAYS_IN_SWEEP) radar->h.vcp = wsr88d_get_volume_coverage(
 388             wsr88d_sweep.ray[i]);
 389         }
 390
 391         free_and_clear_sweep(&wsr88d_sweep, 0, MAX_RAYS_IN_SWEEP);
 392       }
 393
 394       for (iv=0; iv<nvolumes; iv++) {
 395         if (rsl_qfield[iv]) {
 396           radar->v[iv]->h.type_str = strdup(field_str[iv]);
 397           radar->v[iv]->h.nsweeps = nsweep;
 398         }
 399       }
 400   }
 401   wsr88d_close(wf);
 402
 403 /*
 404  * Here we will assign the Radar_header information.  Take most of it
 405  * from an existing volume's header.
 406  */
 407   radar_load_date_time(radar);  /* Magic :-) */
 408
 409     radar->h.number = sitep->number;
 410     memcpy(&radar->h.name, sitep->name, sizeof(sitep->name));
 411     memcpy(&radar->h.radar_name, sitep->name, sizeof(sitep->name)); /* Redundant */
 412     memcpy(&radar->h.city, sitep->city, sizeof(sitep->city));
 413     memcpy(&radar->h.state, sitep->state, sizeof(sitep->state));
 414     strcpy(radar->h.radar_type, "wsr88d");
 415     radar->h.latd = sitep->latd;
 416     radar->h.latm = sitep->latm;
 417     radar->h.lats = sitep->lats;
 418     if (radar->h.latd < 0) { /* Degree/min/sec  all the same sign */
 419       radar->h.latm *= -1;
 420       radar->h.lats *= -1;
 421     }
 422     radar->h.lond = sitep->lond;
 423     radar->h.lonm = sitep->lonm;
 424     radar->h.lons = sitep->lons;
 425     if (radar->h.lond < 0) { /* Degree/min/sec  all the same sign */
 426       radar->h.lonm *= -1;
 427       radar->h.lons *= -1;
 428     }
 429     radar->h.height = sitep->height;
 430     radar->h.spulse = sitep->spulse;
 431     radar->h.lpulse = sitep->lpulse;
 432
 433     free(sitep);
 434
 435   if (wsr88d_merge_split_cuts_is_set()) radar = wsr88d_merge_split_cuts(radar);
 436   return radar;
 437 }