Pileus Git - ~andy/rsl/blob - radar_to_uf.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 #include <time.h>
  26 #include <stdlib.h>
  27
  28 #define USE_RSL_VARS
  29 #include "rsl.h"
  30 extern int radar_verbose_flag;
  31 /* Missing data flag : -32768 when a signed short. */
  32 #define UF_NO_DATA 0X8000
  33
  34
  35 /* Field names. Any convensions may be observed. */
  36 /* Typically:
  37  *    DZ = Reflectivity (dBZ).
  38  *    VR = Radial Velocity.
  39  *    SW = Spectrum Width.
  40  *    CZ = Corrected Reflectivity. (Quality controlled: AP removed, etc.)
  41  *    ZT = Total Reflectivity (dB(mW)).  Becomes UZ in UF files.
  42  *    DR = Differential Reflectivity.
  43  *    LR = Another DR (LDR).
  44  *    ZD = Tina Johnson use this one.
  45  *    DM = Received power.
  46  *    RH = Rho coefficient.
  47  *    PH = Phi (MCTEX parameter).
  48  *    XZ = X-band reflectivity.
  49  *    CD = Corrected ZD.
  50  *    MZ = DZ mask for 1C-51 HDF.
  51  *    MD = ZD mask for 1C-51 HDF.
  52  *    ZE = Edited reflectivity.
  53  *    VE = Edited velocity.
  54  *    KD = KDP wavelength*deg/km
  55  *    TI = TIME (units unknown).
  56  * These fields may appear in any order in the UF file.
  57  * There are more fields than appear here.  See rsl.h.
  58  */
  59
  60
  61 /* Changed old buffer size (16384) for larger dualpol files.  BLK 5/20/2011 */
  62 typedef short UF_buffer[20000]; /* Bigger than documented 4096. */
  63
  64 void swap_uf_buffer(UF_buffer uf);
  65 void swap2(short *buf, int n);
  66 Radar *wsr88d_align_split_cut_rays(Radar *radar);
  67
  68
  69 /**********************************************************************/
  70 /*                                                                    */
  71 /*                     RSL_radar_to_uf_fp                             */
  72 /*                                                                    */
  73 /*  By: John Merritt                                                  */
  74 /*      Space Applications Corporation                                */
  75 /*      May  20, 1994                                                 */
  76 /**********************************************************************/
  77 void RSL_radar_to_uf_fp(Radar *r, FILE *fp)
  78
  79 {
  80   /*
  81    * 1. Fill the UF buffers with data from the Radar structure.
  82    * 2. Write to a stream.  Assume open and leave it so.
  83    */
  84
  85
  86   UF_buffer uf;
  87
  88 /* These are pointers to various locations within the UF buffer 'uf'.
  89  * They are used to index the different components of the UF structure in
  90  * a manor consistant with the UF documentation.  For instance, uf_ma[1]
  91  * will be equivalenced to the second word (2 bytes/each) of the UF
  92  * buffer.
  93  */
  94   short *uf_ma;  /* Mandatory header block. */
  95   short *uf_op;  /* Optional header block.  */
  96   short *uf_lu;  /* Local Use header block.  */
  97   short *uf_dh;  /* Data header.  */
  98   short *uf_fh;  /* Field header. */
  99   short *uf_data; /* Data. */
 100
 101 /* The length of each header. */
 102   int len_ma, len_op, len_lu, len_dh, len_fh, len_data;
 103
 104 /* Booleans to flag inclusion of headers. */
 105   int q_op, q_lu, q_dh, q_fh;
 106
 107   int current_fh_index;
 108   int scale_factor;
 109   int rec_len, save_rec_len;
 110   int nfield;
 111   float vr_az;
 112   int max_field_names;
 113
 114   struct tm *tm;
 115   time_t the_time;
 116
 117   int i,j,k,m;
 118   int degree, minute;
 119   float second;
 120
 121   int uf_sweep_mode = 1;  /* default PPI */
 122
 123 /* Here are the arrays for each field type.  Each dimension is the number
 124  * of fields in the radar structure.  I do this because the radar organization
 125  * is by volumes (field types) and the UF demands that each ray contain
 126  * all the field types.
 127  */
 128   Volume **volume;
 129   Sweep  **sweep;
 130   Ray     *ray;
 131   int     *nsweeps;
 132   int nvolumes, maxsweeps, nrays;
 133   int true_nvolumes;
 134   int sweep_num, ray_num, rec_num;
 135   float x;
 136
 137   Radar *r_save = NULL;
 138
 139   if (r == NULL) {
 140     fprintf(stderr, "radar_to_uf_fp: radar pointer NULL\n");
 141     return;
 142   }
 143
 144   /* If this is WSR-88D data, reorder VR and SW rays in split cuts so that
 145    * their azimuths agree with DZ rays.
 146    */
 147   if (r->h.vcp > 0) {
 148     if (wsr88d_merge_split_cuts_is_set()) {
 149       /* Save a copy of the input radar; we'll restore it later. */
 150       r_save = r;
 151       r = wsr88d_align_split_cut_rays(r);
 152     }
 153   }
 154
 155 /* Do all the headers first time around.  Then, prune OP and LU. */
 156   q_op = q_lu = q_dh = q_fh = 1;
 157
 158   memset(&uf, 0, sizeof(uf)); /* Init to 0 or NULL for pointers. */
 159
 160   sweep_num = ray_num = rec_num = 0;
 161   true_nvolumes = nvolumes = maxsweeps = nrays = 0;
 162
 163   /* PPI and RHI are enum constants defined in rsl.h */
 164   if (r->h.scan_mode == PPI) uf_sweep_mode = 1;
 165   else if (r->h.scan_mode == RHI) uf_sweep_mode = 3;
 166
 167 /*
 168  * The organization of the Radar structure is by volumes, then sweeps, then
 169  * rays, then gates.  This is different from the UF file organization.
 170  * The UF format wants sweeps, rays, then gates for all field types (volumes).
 171  * So, we have to do a back flip, here.  This is achieved by maintaining
 172  * an array of volume pointers and sweep pointers, each dimensioned by
 173  * 'nvolumes', which contains the data for the different field types; this
 174  * is our innermost loop.  The variables are 'volume[i]' and 'sweep[i]' where
 175  * 'i' is the volume index.
 176  *
 177  * In other words, we are getting all the field types together, when we
 178  * are looping on the number of rays in a sweep, so we can load the UF_buffer
 179  * appropriately.
 180  */
 181
 182   nvolumes = r->h.nvolumes;
 183   volume   = (Volume **) calloc(nvolumes, sizeof(Volume *));
 184   sweep    = (Sweep  **) calloc(nvolumes, sizeof(Sweep  *));
 185   nsweeps  = (int *)     calloc(nvolumes, sizeof(int));
 186
 187 /* Get the the number of sweeps in the radar structure.  This will be
 188  * the main controlling loop variable.
 189  */
 190   for (i=0; i<nvolumes; i++) {
 191     volume[i] = r->v[i];
 192     if(volume[i]) {
 193       nsweeps[i] = volume[i]->h.nsweeps;
 194       if (nsweeps[i] > maxsweeps) maxsweeps = nsweeps[i];
 195       true_nvolumes++;
 196     }
 197   }
 198
 199   if (radar_verbose_flag) {
 200     fprintf(stderr,"True number of volumes for UF is %d\n", true_nvolumes);
 201     fprintf(stderr,"Maximum #   of volumes for UF is %d\n", nvolumes);
 202   }
 203
 204   max_field_names = sizeof(RSL_ftype) / 4;
 205
 206 /*--------
 207  *   LOOP for all sweeps (typically 11 or 16 for wsr88d data.
 208  *
 209  */
 210   for (i=0; i<maxsweeps; i++) {
 211     /* Get the array of volume and sweep pointers; one for each field type. */
 212     nrays = 0;
 213     for (k=0; k<nvolumes; k++) {
 214       if (volume[k]) sweep[k] = volume[k]->sweep[i];
 215
 216       /* Check if we really can access this sweep.  Paul discovered that
 217        * if the actual number of sweeps is less than the maximum that we
 218        * could be chasing a bad pointer (a NON-NULL garbage pointer).
 219        */
 220       if (i >= nsweeps[k]) sweep[k] = NULL;
 221
 222       if (sweep[k]) if (sweep[k]->h.nrays > nrays) nrays = sweep[k]->h.nrays;
 223     }
 224
 225     sweep_num++;  /* I guess it will be ok to count NULL sweeps. */
 226     ray_num = 0;
 227   if (radar_verbose_flag)
 228     fprintf(stderr,"Processing sweep %d for %d rays.", i, nrays);
 229   if (radar_verbose_flag)
 230     if (little_endian()) fprintf(stderr," ... On Little endian.\n");
 231     else fprintf(stderr,"\n");
 232
 233
 234 /* Now LOOP for all rays within this particular sweep (i).
 235  *    Get all the field types together for the ray, see ray[k], and
 236  *    fill the UF data buffer appropriately.
 237  */
 238     for (j=0; j<nrays; j++) {
 239       memset(uf, 0, sizeof(uf));
 240       nfield = 0;
 241       ray_num++;  /* And counting, possibly, NULL rays. */
 242       current_fh_index = 0;
 243
 244
 245       /* Find any ray for header information. It does not matter which
 246        * ray, since the information for the MANDITORY, OPTIONAL, and LOCAL
 247        * USE headers is common to any field type ray.
 248        */
 249       ray = NULL;
 250       for (k=0; k<nvolumes; k++) {
 251         if (sweep[k])
 252           if (j < sweep[k]->h.nrays)
 253             if (sweep[k]->ray)
 254               if ((ray = sweep[k]->ray[j])) break;
 255       }
 256
 257       /* If there is no such ray, then continue on to the next ray. */
 258       if (ray) {
 259 /*
 260                   fprintf(stderr,"Ray: %.4d, Time: %2.2d:%2.2d:%f  %.2d/%.2d/%.4d\n", ray_num, ray->h.hour, ray->h.minute, ray->h.sec, ray->h.month, ray->h.day, ray->h.year);
 261 */
 262
 263         /*
 264          * ---- Begining of MANDITORY HEADER BLOCK.
 265          */
 266         uf_ma = uf;
 267         memcpy(&uf_ma[0], "UF", 2);
 268         if (little_endian()) memcpy(&uf_ma[0], "FU", 2);
 269         uf_ma[1]  = 0;  /* Not known yet. */
 270         uf_ma[2]  = 0;  /* Not known yet. Really, I do. */
 271         uf_ma[3]  = 0;  /* Not known yet. */
 272         uf_ma[4]  = 0;  /* Not known yet. */
 273
 274         uf_ma[6]  = 1;
 275         uf_ma[7]  = ray_num;
 276         uf_ma[8 ] = 1;
 277         uf_ma[9 ] = sweep_num;
 278         memcpy(&uf_ma[10], r->h.radar_name, 8);
 279         if (little_endian()) swap2(&uf_ma[10], 8/2);
 280         memcpy(&uf_ma[14], r->h.name, 8);
 281         if (little_endian()) swap2(&uf_ma[14], 8/2);
 282         /* Convert decimal lat/lon to d:m:s */
 283
 284         if (ray->h.lat != 0.0) {
 285           degree = (int)ray->h.lat;
 286           minute = (int)((ray->h.lat - degree) * 60);
 287           second = (ray->h.lat - degree - minute/60.0) * 3600.0;
 288         } else {
 289           degree = r->h.latd;
 290           minute = r->h.latm;
 291           second = r->h.lats;
 292         }
 293         uf_ma[18] = degree;
 294         uf_ma[19] = minute;
 295         if (second > 0.0) uf_ma[20] = second*64 + 0.5;
 296         else uf_ma[20] = second*64 - 0.5;
 297
 298         if (ray->h.lon != 0.0) {
 299           degree = (int)ray->h.lon;
 300           minute = (int)((ray->h.lon - degree) * 60);
 301           second = (ray->h.lon - degree - minute/60.0) * 3600.0;
 302         } else {
 303           degree = r->h.lond;
 304           minute = r->h.lonm;
 305           second = r->h.lons;
 306         }
 307         uf_ma[21] = degree;
 308         uf_ma[22] = minute;
 309         if (second > 0.0) uf_ma[23] = second*64 + 0.5;
 310         else uf_ma[23] = second*64 - 0.5;
 311         if (ray->h.alt != 0)
 312           uf_ma[24] = ray->h.alt;
 313         else
 314           uf_ma[24] = r->h.height;
 315
 316         uf_ma[25] = ray->h.year % 100; /* By definition: not year 2000 compliant. */
 317         uf_ma[26] = ray->h.month;
 318         uf_ma[27] = ray->h.day;
 319         uf_ma[28] = ray->h.hour;
 320         uf_ma[29] = ray->h.minute;
 321         uf_ma[30] = ray->h.sec;
 322         memcpy(&uf_ma[31], "UT", 2);
 323         if (little_endian()) memcpy(&uf_ma[31], "TU", 2);
 324         if (ray->h.azimuth > 0) uf_ma[32] = ray->h.azimuth*64 + 0.5;
 325         else uf_ma[32] = ray->h.azimuth*64 - 0.5;
 326         uf_ma[33] = ray->h.elev*64 + 0.5;
 327         uf_ma[34] = uf_sweep_mode;
 328         if (ray->h.fix_angle != 0.)
 329              uf_ma[35] = ray->h.fix_angle*64.0 + 0.5;
 330         else uf_ma[35] = sweep[k]->h.elev*64.0 + 0.5;
 331         uf_ma[36] = ray->h.sweep_rate*(360.0/60.0)*64.0 + 0.5;
 332
 333         the_time = time(NULL);
 334         tm = gmtime(&the_time);
 335
 336         uf_ma[37] = tm->tm_year % 100; /* Same format as data year */
 337         uf_ma[38] = tm->tm_mon+1;
 338         uf_ma[39] = tm->tm_mday;
 339         memcpy(&uf_ma[40], "RSL" RSL_VERSION_STR, 8);
 340         if (little_endian()) swap2(&uf_ma[40], 8/2);
 341         uf_ma[44] = (signed short)UF_NO_DATA;
 342         len_ma = 45;
 343         uf_ma[2] = len_ma+1;
 344         /*
 345          * ---- End of MANDITORY HEADER BLOCK.
 346          */
 347
 348         /* ---- Begining of OPTIONAL HEADER BLOCK. */
 349         len_op = 0;
 350         if (q_op) {
 351           q_op = 0;  /* Only once. */
 352           uf_op = uf+len_ma;
 353           memcpy(&uf_op[0], "TRMMGVUF", 8);
 354           if (little_endian()) swap2(&uf_op[0], 8/2);
 355           uf_op[4] = (signed short)UF_NO_DATA;
 356           uf_op[5] = (signed short)UF_NO_DATA;
 357           uf_op[6] = r->h.hour;
 358           uf_op[7] = r->h.minute;
 359           uf_op[8] = r->h.sec;
 360           memcpy(&uf_op[9], "RADAR_UF", 8);
 361           if (little_endian()) swap2(&uf_op[9], 8/2);
 362           uf_op[13] = 2;
 363           len_op = 14;
 364         }
 365         /* ---- End of OPTIONAL HEADER BLOCK. */
 366
 367         /* ---- Begining of LOCAL USE HEADER BLOCK. */
 368         q_lu = 0;
 369         len_lu = 0;
 370         if (q_lu) {
 371           uf_lu = uf+len_ma+len_op;
 372           q_lu = 0; /* Only once. */
 373           len_lu = 0;  /* I don't have anything yet. */
 374         }
 375         /* ---- End  of LOCAL USE HEADER BLOCK. */
 376
 377
 378        /* Here is where we loop on each field type.  We need to keep
 379         * track of how many FIELD HEADER and FIELD DATA sections, one
 380         * for each field type, we fill.  The variable that tracks this
 381         * index into 'uf' is 'current_fh_index'.  It is bumped by
 382         * the length of the FIELD HEADER and FIELD DATA for each field
 383         * type encountered.  Field types expected are: Reflectivity,
 384         * Velocity, and Spectrum width; this is a typicial list but it
 385         * is not restricted to it.
 386         */
 387
 388          for (k=0; k<nvolumes; k++) {
 389           if (sweep[k])
 390             if (j < sweep[k]->h.nrays && sweep[k]->ray[j])
 391               ray = sweep[k]->ray[j];
 392             else
 393               ray = NULL;
 394           else ray = NULL;
 395
 396           if (ray) {
 397             /* ---- Begining of DATA HEADER. */
 398             nfield++;
 399             if (q_dh) {
 400               len_dh = 2*true_nvolumes + 3;
 401               uf_dh = uf+len_ma+len_op+len_lu;
 402               uf_dh[0] = nfield;
 403               uf_dh[1] = 1;
 404               uf_dh[2] = nfield;
 405               /* 'nfield' indexes the field number.
 406                * 'k' indexes the particular field from the volume.
 407                *  RSL_ftype contains field names and is defined in rsl.h.
 408                */
 409               if (k > max_field_names-1) {
 410                 fprintf(stderr,
 411                   "RSL_uf_to_radar: No field name for volume index %d\n", k);
 412                 fprintf(stderr,"RSL_ftype must be updated in rsl.h for new field.\n");
 413                 fprintf(stderr,"Quitting now.\n");
 414                 return;
 415               }
 416               memcpy(&uf_dh[3+2*(nfield-1)], RSL_ftype[k], 2);
 417               if (little_endian()) swap2(&uf_dh[3+2*(nfield-1)], 2/2);
 418               if (current_fh_index == 0) current_fh_index = len_ma+len_op+len_lu+len_dh;
 419               uf_dh[4+2*(nfield-1)] = current_fh_index + 1;
 420             }
 421             /* ---- End of DATA HEADER. */
 422
 423             /* ---- Begining of FIELD HEADER. */
 424             if (q_fh) {
 425               uf_fh = uf+current_fh_index;
 426               if (k != PH_INDEX) scale_factor = 100;
 427               else scale_factor = 10;
 428               uf_fh[1] = scale_factor;
 429               uf_fh[2] = ray->h.range_bin1/1000.0;
 430               uf_fh[3] = ray->h.range_bin1 - (1000*uf_fh[2]);
 431               uf_fh[4] = ray->h.gate_size;
 432               uf_fh[5] = ray->h.nbins;
 433               uf_fh[6] = ray->h.pulse_width*(RSL_SPEED_OF_LIGHT/1.0e6);
 434               uf_fh[7] = sweep[k]->h.beam_width*64.0 + 0.5;
 435               uf_fh[8] = sweep[k]->h.beam_width*64.0 + 0.5;
 436               uf_fh[9] = ray->h.frequency * 1000.; /* Bandwidth (mHz). */
 437               uf_fh[10] = 0; /* Horizontal polarization. */
 438               uf_fh[11] = ray->h.wavelength*64.0*100.0; /* m to cm. */
 439               uf_fh[12] = ray->h.pulse_count;
 440               memcpy(&uf_fh[13], "  ", 2);
 441               uf_fh[14] = (signed short)UF_NO_DATA;
 442               uf_fh[15] = (signed short)UF_NO_DATA;
 443               if (k == DZ_INDEX || k == ZT_INDEX) {
 444                 uf_fh[16] = volume[k]->h.calibr_const*100.0 + 0.5;
 445               }
 446               else {
 447                 memcpy(&uf_fh[16], "  ", 2);
 448               }
 449               if (ray->h.prf != 0)
 450                 uf_fh[17] = 1.0/ray->h.prf*1000000.0; /* Pulse repetition time(msec)  = 1/prf */
 451               else
 452                 uf_fh[17] = (signed short)UF_NO_DATA; /* Pulse repetition time  = 1/prf */
 453               uf_fh[18] = 16;
 454               if (VR_INDEX == k || VE_INDEX == k) {
 455                 uf_fh[19] = scale_factor*ray->h.nyq_vel;
 456                 uf_fh[20] = 1;
 457                 len_fh = 21;
 458               } else {
 459                 len_fh = 19;
 460               }
 461
 462               uf_fh[0] = current_fh_index + len_fh + 1;
 463               /* ---- End of FIELD HEADER. */
 464
 465               /* ---- Begining of FIELD DATA. */
 466               uf_data = uf+len_fh+current_fh_index;
 467               len_data = ray->h.nbins;
 468               for (m=0; m<len_data; m++) {
 469                 x = ray->h.f(ray->range[m]);
 470                 if (x == BADVAL || x == RFVAL || x == APFLAG || x == NOECHO)
 471                   uf_data[m] = (signed short)UF_NO_DATA;
 472                 else
 473                   uf_data[m] = scale_factor * x;
 474               }
 475
 476               current_fh_index += (len_fh+len_data);
 477             }
 478           }
 479         /* ---- End of FIELD DATA. */
 480         }
 481         /* Fill in some infomation we didn't know.  Like, buffer length,
 482          * record number, etc.
 483          */
 484         rec_num++;
 485         uf_ma[1] = current_fh_index;
 486         uf_ma[3] = len_ma + len_op + 1;
 487         uf_ma[4] = len_ma + len_op + len_lu + 1;
 488         uf_ma[5] = rec_num;
 489
 490         /* WRITE the UF buffer. */
 491         rec_len =(int)uf_ma[1]*2;
 492         save_rec_len = rec_len;  /* We destroy 'rec_len' when making it
 493                         big endian on a little endian machine. */
 494         if (little_endian()) swap_4_bytes(&rec_len);
 495         (void)fwrite(&rec_len, sizeof(int), 1, fp);
 496         if (little_endian()) swap_uf_buffer(uf);
 497         (void)fwrite(uf, sizeof(char), save_rec_len, fp);
 498         (void)fwrite(&rec_len, sizeof(int), 1, fp);
 499       } /* if (ray) */
 500     }
 501   }
 502
 503   /* If Radar argument "r" was modified (WSR-88D), restore the saved copy. */
 504   if (r_save != NULL) r = r_save;
 505 }
 506
 507 /**********************************************************************/
 508 /*                                                                    */
 509 /*                     RSL_radar_to_uf                                */
 510 /*                                                                    */
 511 /**********************************************************************/
 512 void RSL_radar_to_uf(Radar *r, char *outfile)
 513 {
 514   FILE *fp;
 515   if (r == NULL) {
 516     fprintf(stderr, "radar_to_uf: radar pointer NULL\n");
 517     return;
 518   }
 519
 520   if ((fp = fopen(outfile, "w")) == NULL) {
 521     perror(outfile);
 522     return;
 523   }
 524
 525   RSL_radar_to_uf_fp(r, fp);
 526   fclose(fp);
 527 }
 528
 529 /**********************************************************************/
 530 /*                                                                    */
 531 /*                     RSL_radar_to_uf_gzip                           */
 532 /*                                                                    */
 533 /**********************************************************************/
 534 void RSL_radar_to_uf_gzip(Radar *r, char *outfile)
 535 {
 536   FILE *fp;
 537   if (r == NULL) {
 538     fprintf(stderr, "radar_to_uf_gzip: radar pointer NULL\n");
 539     return;
 540   }
 541
 542   if ((fp = fopen(outfile, "w")) == NULL) {
 543     perror(outfile);
 544     return;
 545   }
 546
 547   fp = compress_pipe(fp);
 548   RSL_radar_to_uf_fp(r, fp);
 549   rsl_pclose(fp);
 550 }