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