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

   1 #ifdef HAVE_CONFIG_H
   2 #include "config.h"
   3 #endif
   4 #ifdef HAVE_LASSEN
   5 /*
   6     NASA/TRMM, Code 910.1.
   7     This is the TRMM Office Radar Software Library.
   8     Copyright (C) 1996, 1997
   9             John H. Merritt
  10             Space Applications Corporation
  11             Vienna, Virginia
  12
  13     This library is free software; you can redistribute it and/or
  14     modify it under the terms of the GNU Library General Public
  15     License as published by the Free Software Foundation; either
  16     version 2 of the License, or (at your option) any later version.
  17
  18     This library is distributed in the hope that it will be useful,
  19     but WITHOUT ANY WARRANTY; without even the implied warranty of
  20     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  21     Library General Public License for more details.
  22
  23     You should have received a copy of the GNU Library General Public
  24     License along with this library; if not, write to the Free
  25     Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  26 */
  27 /*----------------------------------------------------------------------
  28  * PURPOSE:
  29  *
  30  * Ingest a Lassen file and fill the Lassen_volume data structure.
  31  * It can handle version 1.3 and 1.4 files (MCTEX data).
  32  *
  33  * The algorthm is pretty simple, call an 'xdr' routine each time for
  34  * each and every member of Lassen_volume and all substructures.
  35  *
  36  * The order and type of xdr routine called defines the file organization.
  37  *----------------------------------------------------------------------
  38  */
  39 #include <stdio.h>
  40 #include <stdlib.h>
  41 #include <string.h>
  42 #include <rpc/rpc.h>
  43 #include <netinet/in.h>
  44 #include "lassen.h"
  45
  46 /* xdr_destroy is broken on HPUX, SGI, and SUN; Linux, the only working one? */
  47 /* This broken behavior could be from old xdr versions, too. */
  48 #undef BROKEN_XDR_DESTROY
  49 #if defined(__hpux) || defined(__sgi) || defined(__sun)
  50 #define BROKEN_XDR_DESTROY
  51 #endif
  52 /*************************************************************/
  53 /*                                                           */
  54 /*                     read_lassen_head                      */
  55 /*                                                           */
  56 /*************************************************************/
  57 int read_lassen_head(FILE *f, Lassen_head *head)
  58 {
  59   int rc=1, i;
  60   char *tmp;
  61   XDR xdrs;
  62
  63   memset(head->magic, 0, sizeof(head->magic));
  64
  65   xdrstdio_create(&xdrs, f, XDR_DECODE);
  66
  67   /*
  68    * Can I read the first string? If not, then this is probably
  69    * not a valid header.  If this happens, then return.
  70    */
  71   tmp=head->magic;
  72   if((rc &= xdr_string(&xdrs, &tmp, 8))==0) {
  73 #ifndef BROKEN_XDR_DESTROY
  74         xdr_destroy(&xdrs);
  75 #endif
  76         return(0);
  77   }
  78   rc &= xdr_u_char(&xdrs, &head->mdate.year);
  79   rc &= xdr_u_char(&xdrs, &head->mdate.month);
  80   rc &= xdr_u_char(&xdrs, &head->mdate.day);
  81   rc &= xdr_u_char(&xdrs, &head->mdate.hour);
  82   rc &= xdr_u_char(&xdrs, &head->mdate.minute);
  83   rc &= xdr_u_char(&xdrs, &head->mdate.second);
  84   rc &= xdr_u_char(&xdrs, &head->cdate.year);
  85   rc &= xdr_u_char(&xdrs, &head->cdate.month);
  86   rc &= xdr_u_char(&xdrs, &head->cdate.day);
  87   rc &= xdr_u_char(&xdrs, &head->cdate.hour);
  88   rc &= xdr_u_char(&xdrs, &head->cdate.minute);
  89   rc &= xdr_u_char(&xdrs, &head->cdate.second);
  90   rc &= xdr_int(&xdrs, &head->type);
  91   tmp=head->mwho;
  92   rc &= xdr_string(&xdrs, &tmp, 16);
  93   tmp=head->cwho;
  94   rc &= xdr_string(&xdrs, &tmp, 16);
  95   rc &= xdr_int(&xdrs, &head->protection);
  96   rc &= xdr_int(&xdrs, &head->checksum);
  97   tmp=head->description;
  98   rc &= xdr_string(&xdrs, &tmp, 40);
  99   rc &= xdr_int(&xdrs, &head->id);
 100   for(i=0;i<12;i++)
 101         rc &= xdr_int(&xdrs, &head->spare[i]);
 102
 103 #ifndef BROKEN_XDR_DESTROY
 104   xdr_destroy(&xdrs);
 105 #endif
 106
 107   return rc;
 108 }
 109
 110 /*************************************************************/
 111 /*                                                           */
 112 /*                   free_lassen_volume                      */
 113 /*                                                           */
 114 /*************************************************************/
 115
 116 void free_lassen_volume(Lassen_volume *vol)
 117 {
 118   int i,j;
 119
 120   if (vol == NULL) return;
 121   for(i=0; i<vol->numsweeps; i++ ) {
 122         if(vol->index[i] == NULL) continue;
 123         for(j=0; j<vol->index[i]->numrays; j++ )
 124           if (vol->index[i]->ray[j] != NULL)
 125                 free( vol->index[i]->ray[j] );
 126         free( vol->index[i] );
 127   }
 128   free( vol );
 129   return;
 130 }
 131
 132 static void short_shift(unsigned short *n)
 133 {
 134   *n = ntohs(*n);
 135 }
 136
 137 /*************************************************************/
 138 /*                                                           */
 139 /*                   read_lassen_volume                      */
 140 /*                                                           */
 141 /*************************************************************/
 142 int read_lassen_volume( XDR *xdr, Lassen_volume *vol)
 143 {
 144   int   i,j;
 145   int   rc = TRUE;
 146   char          *tmp;
 147   unsigned short        ushort_tmp;
 148
 149   /*  Read the version number. */
 150   rc &= xdr_u_short(xdr, &vol->version);
 151
 152   /*
 153    *  Currently, Version 1.3 and 1.4 seem to be the same on disk.
 154    *  I have no doc's.  The only thing that I know of to be different
 155    *  is how the ray azimuth angles are calculated.  For MCTEX data,
 156    *  this is time dependant.  The correction will happen in the
 157    *  application; after RSL_lassen_to_radar is called, for example.
 158    *
 159    *  Version 1.3 is coded as 13 and version 1.4 is coded as 14.
 160    */
 161   if( (int)vol->version != 13 && (int)vol->version != 14 ) {
 162         (void)fprintf(stderr, "Incompatible version of LASSEN file.\n" );
 163         (void)fprintf(stderr, "File version: %d.%d\n",
 164                                   (int)vol->version/10, (int)vol->version%10);
 165         return( 0 );
 166   }
 167
 168   /*  Read in the Lassen_volume.  */
 169   rc &= xdr_short(xdr, &vol->filled);
 170   rc &= xdr_u_int(xdr, &vol->volume);
 171   rc &= xdr_u_short(xdr, &vol->sweep);
 172   rc &= xdr_u_short(xdr, &vol->sweep_type);
 173   rc &= xdr_u_short(xdr, &vol->max_height);
 174   rc &= xdr_u_short(xdr, &vol->status);
 175   rc &= xdr_u_short(xdr, &vol->min_fangle);
 176   rc &= xdr_u_short(xdr, &vol->max_fangle);
 177   rc &= xdr_u_short(xdr, &vol->min_var);
 178   rc &= xdr_u_short(xdr, &vol->max_var);
 179   rc &= xdr_u_short(xdr, &vol->a_start);
 180   rc &= xdr_u_short(xdr, &vol->a_stop);
 181   rc &= xdr_u_short(xdr, &vol->numsweeps);
 182   for(i=0;i<LASSEN_MAX_SWEEPS;i++)
 183         rc &= xdr_u_short(xdr, &vol->fangles[i]);
 184   rc &= xdr_u_short(xdr, &vol->gatewid);
 185   rc &= xdr_u_short(xdr, &vol->rangeg1);
 186   for(i=0;i<LASSEN_MAX_SWEEPS;i++)
 187         rc &= xdr_u_short(xdr, &vol->numgates[i]);
 188   rc &= xdr_u_short(xdr, &vol->maxgates);
 189   rc &= xdr_u_short(xdr, &ushort_tmp);
 190   rc &= xdr_u_short(xdr, &ushort_tmp);
 191   rc &= xdr_u_short(xdr, &ushort_tmp);
 192   rc &= xdr_u_short(xdr, &ushort_tmp);
 193   rc &= xdr_u_short(xdr, &vol->prf);
 194   rc &= xdr_u_short(xdr, &vol->prflow);
 195   rc &= xdr_u_int(xdr, &vol->freq);
 196   rc &= xdr_u_short(xdr, &vol->n_pulses);
 197
 198   for(i=0;i<LASSEN_MAX_SWEEPS;i++)
 199         for(j=0;j<NUMOFFSETS;j++)
 200           rc &= xdr_u_short(xdr, &vol->offset[i][j]);
 201
 202   rc &= xdr_u_char(xdr, &vol->year);
 203   rc &= xdr_u_char(xdr, &vol->month);
 204   rc &= xdr_u_char(xdr, &vol->day);
 205   rc &= xdr_u_char(xdr, &vol->shour);
 206   rc &= xdr_u_char(xdr, &vol->sminute);
 207   rc &= xdr_u_char(xdr, &vol->ssecond);
 208   rc &= xdr_u_char(xdr, &vol->ehour);
 209   rc &= xdr_u_char(xdr, &vol->eminute);
 210   rc &= xdr_u_char(xdr, &vol->esecond);
 211   /* The flags use bits.  The data type length is short; I expect. */
 212   rc &= xdr_u_short(xdr, (unsigned short *)&vol->volflags);
 213   short_shift((unsigned short *)&vol->volflags);
 214   /* This is radar info. */
 215   tmp=vol->radinfo.radar_name;
 216   rc &= xdr_string(xdr, &tmp, 8);
 217   tmp=vol->radinfo.site_name;
 218   rc &= xdr_string(xdr, &tmp, 8);
 219   rc &= xdr_u_short(xdr, &vol->radinfo.antenna_height);
 220   rc &= xdr_short(xdr, &vol->radinfo.latitude.degree);
 221   rc &= xdr_short(xdr, &vol->radinfo.latitude.minute);
 222   rc &= xdr_short(xdr, &vol->radinfo.latitude.second);
 223   rc &= xdr_short(xdr, &vol->radinfo.longitude.degree);
 224   rc &= xdr_short(xdr, &vol->radinfo.longitude.minute);
 225   rc &= xdr_short(xdr, &vol->radinfo.longitude.second);
 226
 227   /*
 228    *  The return value is the logical and of all the xdr reads.
 229    */
 230   return( rc );
 231 }
 232
 233 /*************************************************************/
 234 /*                                                           */
 235 /*                   read_lassen_ray                         */
 236 /*                                                           */
 237 /*************************************************************/
 238 int read_lassen_ray(XDR *xdr, Lassen_ray *ray )
 239 {
 240         int     i, rc = 1;
 241
 242         /*  Read in Lassen_ray */
 243         rc &= xdr_u_short( xdr, &ray->vangle);
 244         rc &= xdr_u_short( xdr, &ray->fanglet);
 245         rc &= xdr_u_short( xdr, &ray->fanglea);
 246         rc &= xdr_u_short( xdr, &ray->a_start);
 247         rc &= xdr_u_short( xdr, &ray->a_stop);
 248         rc &= xdr_u_char( xdr, &ray->max_height);
 249         rc &= xdr_u_char( xdr, &ray->volume);
 250         rc &= xdr_u_char( xdr, &ray->sweep);
 251         rc &= xdr_u_char( xdr, &ray->sweep_type);
 252         rc &= xdr_u_short( xdr, &ray->gatewid);
 253         rc &= xdr_u_short( xdr, &ray->rangeg1);
 254         rc &= xdr_u_short( xdr, &ray->numgates);
 255         rc &= xdr_u_short( xdr, &ray->prf);
 256         rc &= xdr_u_short( xdr, &ray->prflow);
 257         rc &= xdr_u_short( xdr, &ray->n_pulses);
 258         rc &= xdr_u_char( xdr, &ray->p_width);
 259         rc &= xdr_u_char( xdr, &ray->cfilter);
 260         rc &= xdr_u_short( xdr, &ray->status);
 261         memset(&ray->flags, 0, 4);
 262         /* I expect the length of the bit field to be a 'short', never 'long' */
 263         rc &= xdr_u_short( xdr, (unsigned short *)&ray->flags); /* a bit field */
 264         short_shift((unsigned short *)&ray->flags);
 265         /*
 266         fprintf(stderr, "short_shift: flags = %x\n", ray->flags);
 267         fprintf(stderr, "packed = %d", (int)ray->flags.packed);
 268         fprintf(stderr, "good_data = %d", (int)ray->flags.good_data);
 269         fprintf(stderr, "uz = %d", (int)ray->flags.uz);
 270         fprintf(stderr, "cz = %d", (int)ray->flags.cz);
 271         fprintf(stderr, "vel= %d", (int)ray->flags.vel);
 272         fprintf(stderr, "wid= %d", (int)ray->flags.wid);
 273         fprintf(stderr, "zdr= %d", (int)ray->flags.zdr);
 274         fprintf(stderr, "phi= %d", (int)ray->flags.phi);
 275         fprintf(stderr, "rho= %d", (int)ray->flags.rho);
 276         fprintf(stderr, "ldr= %d", (int)ray->flags.ldr);
 277         fprintf(stderr, "kdp= %d", (int)ray->flags.kdp);
 278         fprintf(stderr, "tim= %d", (int)ray->flags.time);
 279         fprintf(stderr, "spa= %x", (int)ray->flags.spares);
 280         */
 281         for(i=0;i<NUMOFFSETS;i++)
 282                 rc &= xdr_u_short( xdr, &ray->offset[i]);
 283
 284         rc &= xdr_u_char( xdr, &ray->year);
 285         rc &= xdr_u_char( xdr, &ray->month);
 286         rc &= xdr_u_char( xdr, &ray->day);
 287         rc &= xdr_u_char( xdr, &ray->hour);
 288         rc &= xdr_u_char( xdr, &ray->minute);
 289         rc &= xdr_u_char( xdr, &ray->second);
 290
 291         return rc;
 292 }
 293
 294
 295 /*************************************************************/
 296 /*                                                           */
 297 /*                  read_entire_lassen_file                  */
 298 /*                                                           */
 299 /*************************************************************/
 300 int read_entire_lassen_file(FILE *f, Lassen_volume *vol)
 301 {
 302   Lassen_sweep  *sweep;
 303   int i;
 304   XDR                           xdr;
 305   Lassen_ray            ray;
 306   Lassen_head           head;
 307   int                           isweep = (-99);
 308   int                           lastray=0,numrays=0, rc;
 309   unsigned int      size;
 310   int tbytes = 0;
 311   unsigned char *p;
 312
 313
 314   /* Skip the header. */
 315   read_lassen_head(f, &head);
 316
 317   /*  Connect to the XDR stream. */
 318   xdrstdio_create(&xdr, f, XDR_DECODE);
 319
 320   /*  Check the volume header. Is the version is correct? */
 321   if( read_lassen_volume( &xdr, vol ) == 0) {
 322 #ifndef BROKEN_XDR_DESTROY
 323         xdr_destroy(&xdr);
 324 #endif
 325         return( 0 );
 326   }
 327
 328   rc = TRUE;
 329
 330   /*  Read until XDR error. */
 331   while(rc == TRUE) {
 332
 333         /*
 334          *  Read in the ray header.  If there is an error in reading
 335          *      in the ray header, then make sure the rc value is
 336          *      set to 0.
 337          */
 338         if( read_lassen_ray( &xdr, &ray ) == 0 )
 339           break;
 340
 341         if( (short)ANGLE_CONVERT(ray.vangle) < 0 ||
 342                 (short)ANGLE_CONVERT(ray.vangle) > 360 ) {
 343           fprintf(stderr, "read_entire_lassen_file: angle out of range (%d)\n",
 344                           ANGLE_CONVERT(ray.vangle));
 345           rc = 0;
 346           break;
 347         }
 348
 349         if(ray.volume != (vol->volume&0x00ff)) {
 350           fprintf(stderr, "read_entire_lassen_file: Volume serial number out of sync:\n");
 351           fprintf(stderr, "         Ray->volume = %u, Should be %u\n",
 352                           ray.volume, (unsigned char)vol->volume&0x00ff );
 353           fprintf(stderr, "         (Ray->volume == Vol->volume(%d) & 0x00ff)\n",
 354                           vol->volume );
 355           if(isweep == -99) rc = 0;
 356           break;
 357         }
 358
 359         /*
 360          * If we have a new sweep, then create
 361          * the new sweep data structure.  The 'isweep'
 362          * contains the offset number in the volume. 'ray.sweep'
 363      * has the sweep number.
 364          */
 365         if((ray.sweep-1) != isweep) {
 366           if(ray.sweep == 0) {
 367                 fprintf(stderr,"Sweep number in a ray is 0!\n");
 368                 fprintf(stderr,"Aborting read_entire_lassen_file()\n");
 369 #ifndef BROKEN_XDR_DESTROY
 370                 xdr_destroy(&xdr);
 371 #endif
 372                 rc = 0;
 373                 break;
 374           }
 375
 376
 377           isweep = ray.sweep-1;
 378
 379           if(isweep > 0) {
 380                 /* Reset the last ray when passing 2 pi. */
 381                 lastray--;
 382                 if( lastray < 0 ) lastray = 359;
 383
 384                 sweep = vol->index[isweep-1];
 385                 if (sweep == NULL) {
 386                   rc = 0;
 387                   fprintf(stderr, "vol->index[%d] is NULL.\n", isweep-1);
 388                   break;
 389                 }
 390                 sweep->ehour   = sweep->ray[lastray]->hour;
 391                 sweep->eminute = sweep->ray[lastray]->minute;
 392                 sweep->esecond = sweep->ray[lastray]->second;
 393                 sweep->numrays = numrays;
 394                 sweep->max_var = sweep->ray[lastray]->vangle;
 395
 396           }
 397
 398           /* We have a new sweep. */
 399           if((vol->index[isweep] = (Lassen_sweep *)calloc
 400                   (1,sizeof(Lassen_sweep))) == 0) {
 401                 perror("read_entire_lassen_file");
 402                 rc = 0;
 403                 break;
 404           }
 405
 406           lastray = 0;
 407
 408           /* Copy to the Sweep. */
 409           sweep          = vol->index[isweep];
 410           sweep->sweep   = isweep;
 411           sweep->shour   = ray.hour;
 412           sweep->sminute = ray.minute;
 413           sweep->ssecond = ray.second;
 414           sweep->volume  = ray.volume;
 415           sweep->sweep   = ray.sweep;
 416           sweep->sweep_type = ray.sweep_type;
 417           sweep->max_height = ray.max_height;
 418           sweep->fangle     = ray.fanglet;
 419           sweep->a_start    = ray.a_start;
 420           sweep->a_stop     = ray.a_stop;
 421           sweep->gatewid    = ray.gatewid;
 422           sweep->rangeg1    = ray.rangeg1;
 423           sweep->numgates   = ray.numgates;
 424           sweep->prf        = ray.prf;
 425           sweep->prflow     = ray.prflow;
 426           sweep->n_pulses   = ray.n_pulses;
 427           sweep->p_width    = ray.p_width;
 428           sweep->cfilter    = ray.cfilter;
 429           for(i=0; i<NUMOFFSETS; i++)
 430                 sweep->offset[i] = ray.offset[i];
 431           sweep->year    = ray.year;
 432           sweep->month   = ray.month;
 433           sweep->day     = ray.day;
 434           sweep->min_var = ray.vangle;
 435
 436           numrays = 0;
 437         }
 438
 439         /*
 440          * If everything is ok, then add the ray to the
 441          * volume.
 442          */
 443         numrays++;
 444         vol->index[isweep]->status |= ray.status;
 445
 446         /* Determine the 'size' of the ray and allocate memory. */
 447         size = 0;
 448         /*                      fprintf(stderr,"\nOFFSETS (%d):", NUMOFFSETS); */
 449         for(i=0; i<NUMOFFSETS; i++) {
 450           size += (ray.numgates*(ray.offset[i] !=0 ));
 451           /*                            fprintf(stderr, " %d", ray.numgates*(ray.offset[i])); */
 452         }
 453
 454         if((vol->index[isweep]->ray[lastray] =
 455                 (Lassen_ray *) calloc(1,sizeof(Lassen_ray)+size))==0) {
 456           perror("read_entire_lassen_file");
 457           rc = 0;
 458 #ifndef BROKEN_XDR_DESTROY
 459           xdr_destroy(&xdr);
 460 #endif
 461           break;
 462         }
 463
 464
 465         p = (unsigned char *)vol->index[isweep]->ray[lastray];
 466
 467         /* Copy the ray header. */
 468         memcpy(p, &ray, sizeof(Lassen_ray));
 469         p = (unsigned char *)p + sizeof(Lassen_ray);
 470
 471         /* Read the data. */
 472         /*                      printf("Previous rc=%d ... ", rc); */
 473         rc &= xdr_bytes(&xdr, (char **)&p, &size,size);
 474         tbytes += size;
 475         /*                      printf("after rc=%d.  xdr_bytes %d bytes, tbytes %d.\n", rc, size, tbytes);
 476          */
 477         lastray++;
 478         if(lastray > 359) lastray = 0;
 479   } /* End while. */
 480
 481   if(rc == 0) {
 482         /* Clean up and get out. */
 483         free_lassen_volume(vol);
 484 #ifndef BROKEN_XDR_DESTROY
 485         xdr_destroy(&xdr);
 486 #endif
 487         fprintf(stderr,"\tAborting read_entire_lassen_file.\n\n");
 488         return(0);
 489   }
 490
 491   /*
 492    * Check to see if we read a sweep!
 493    */
 494   if( isweep == -99 )
 495         vol->filled = -1;
 496   else {
 497         lastray--;
 498         if( lastray < 0 ) lastray = 359;
 499
 500         /* Update the last sweep. */
 501         sweep          = vol->index[isweep];
 502         sweep->ehour   = sweep->ray[lastray]->hour;
 503         sweep->eminute = sweep->ray[lastray]->minute;
 504         sweep->esecond = sweep->ray[lastray]->second;
 505         sweep->max_var = ray.vangle;
 506         sweep->numrays = numrays;
 507         vol->filled = 1;
 508   }
 509
 510   /* Check the number of tilts read vs. the number in the file. */
 511   if(vol->numsweeps != isweep+1) {
 512         fprintf(stderr, "read_entire_lassen_file: Warning, read error.  The number of sweeps ");
 513         fprintf(stderr, "detected is different from the number written.");
 514   }
 515
 516   vol->numsweeps = isweep+1;
 517
 518 #ifndef BROKEN_XDR_DESTROY
 519   xdr_destroy(&xdr);
 520 #endif
 521   return( 1 );
 522 }
 523 #endif