Pileus Git - ~andy/rsl/blob - lassen_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
  24 19 Jan 1998
  25    Michael Whimpey from BMRC, Australia changed code, for the
  26    different callibrations, between, Pre_mctex, mctex, Gunn_Pt
  27    periods.
  28 22 Apr 1999
  29    Michael Whimpey added more callibration periods over berrimah and scsmex.
  30    please see code where there is m.whimpey
  31
  32 */
  33 #include <stdio.h>
  34 #include <stdlib.h>
  35 #include <unistd.h>
  36 #include <string.h>
  37 #include <math.h>
  38
  39 #define USE_RSL_VARS
  40 #include "rsl.h"
  41
  42 #ifdef HAVE_LASSEN
  43 #include "lassen.h"
  44 extern int radar_verbose_flag;
  45
  46 /* Some parameter headaches are prevented when vol is declared global. */
  47 Lassen_volume vol;
  48 extern int read_entire_lassen_file(FILE *f, Lassen_volume *vol);
  49
  50 /**********************************************************************/
  51 /*                                                                    */
  52 /*                     lassen_load_sweep                              */
  53 /*                                                                    */
  54 /*  By: John Merritt                                                  */
  55 /*      Space Applications Corporation                                */
  56 /*      May  26, 1994                                                 */
  57 /**********************************************************************/
  58 void lassen_load_sweep(Sweep *s, int isweep_num, int ifield, int period, Lassen_sweep *ptr)
  59 {
  60   float c = RSL_SPEED_OF_LIGHT;
  61   float elev;
  62   Lassen_ray *aray;
  63   unsigned char *ray_data;
  64   int i,j,m,kk=0;
  65   float x[2000];
  66   double Vu;
  67   Range (*invf)(float x);
  68   float (*f)(Range x);
  69
  70 /*  calibration period flags   m.whimpey  */
  71 #define BERRIMAH1 0
  72 #define MCTEX_EARLY 1
  73 #define MCTEX 2
  74 #define BERRIMAH2 3
  75 #define GUNN_PT 4
  76 #define GUNN_PT1 5
  77 #define SCSMEX 6
  78 #define SCSMEX1 7
  79
  80   if (s == NULL) return;
  81   f    = (float (*)(Range x))NULL;  /* This quiets the pedantic warning. */
  82   invf = (Range (*)(float x))NULL;
  83   if (ifield == ZT_INDEX) {kk = OFF_UZ;  invf = ZT_INVF; f = ZT_F;}
  84   if (ifield == DZ_INDEX) {kk = OFF_CZ;  invf = DZ_INVF; f = DZ_F;}
  85   if (ifield == VR_INDEX) {kk = OFF_VEL; invf = VR_INVF; f = VR_F;}
  86   if (ifield == SW_INDEX) {kk = OFF_WID; invf = SW_INVF; f = SW_F;}
  87   if (ifield == ZD_INDEX) {kk = OFF_ZDR; invf = ZD_INVF; f = ZD_F;}
  88   if (ifield == PH_INDEX) {kk = OFF_PHI; invf = PH_INVF; f = PH_F;}
  89   if (ifield == RH_INDEX) {kk = OFF_RHO; invf = RH_INVF; f = RH_F;}
  90   if (ifield == LR_INDEX) {kk = OFF_LDR; invf = LR_INVF; f = LR_F;}
  91   if (ifield == KD_INDEX) {kk = OFF_KDP; invf = KD_INVF; f = KD_F;}
  92   if (ifield == TI_INDEX) {kk = OFF_TIME; invf = TI_INVF; f = TI_F;}
  93
  94   elev = (float)ptr->fangle*360.0/16384.0;
  95   Vu = c*((float)vol.prf/10.)/(4.*(float)vol.freq*100000.0);
  96
  97   s->h.sweep_num = ptr->sweep;
  98   s->h.elev = elev;
  99   s->h.nrays = ptr->numrays;
 100   s->h.beam_width = 1.0;  /* What is it really? */
 101   s->h.horz_half_bw = .5;
 102   s->h.vert_half_bw = .5;
 103   s->h.f = f;
 104   s->h.invf = invf;
 105
 106   for(i=0;i<(int)ptr->numrays;i++) {
 107         aray = ptr->ray[i];
 108         if (aray == NULL) continue;
 109         s->ray[i] = RSL_new_ray((int)aray->numgates);
 110         s->ray[i]->h.month    = aray->month;
 111         s->ray[i]->h.day      = aray->day;
 112         s->ray[i]->h.year     = (int)aray->year + 1900;
 113         if (s->ray[i]->h.year < 1980) s->ray[i]->h.year += 100; /* Year > 2000. */
 114         s->ray[i]->h.hour     = aray->hour;
 115         s->ray[i]->h.minute   = aray->minute;
 116         s->ray[i]->h.sec      = aray->second;
 117         s->ray[i]->h.azimuth  =  (float)aray->vangle*360.0/16384.0;
 118
 119         s->ray[i]->h.ray_num    = i;
 120         s->ray[i]->h.elev       = elev;
 121         s->ray[i]->h.elev_num   = s->h.sweep_num;
 122         s->ray[i]->h.range_bin1 = aray->rangeg1;
 123         s->ray[i]->h.gate_size  = aray->gatewid;
 124         s->ray[i]->h.vel_res  = 0.5; /* What is this really? */
 125
 126         s->ray[i]->h.fix_angle = s->h.elev;
 127         s->ray[i]->h.frequency =        (float)vol.freq*1.0e-4;  /* GHz */
 128         s->ray[i]->h.wavelength = c / s->ray[i]->h.frequency * 1.0e-9;
 129         s->ray[i]->h.prf      = (int)aray->prf/10;
 130         s->ray[i]->h.nyq_vel = s->ray[i]->h.prf * s->ray[i]->h.wavelength / 4.0;
 131         if (s->ray[i]->h.prf != 0)
 132           s->ray[i]->h.unam_rng = c / (2.0 * s->ray[i]->h.prf * 1000.0); /* km */
 133         else
 134           s->ray[i]->h.unam_rng = 0.0;
 135         s->ray[i]->h.pulse_width = (float)(aray->p_width * 0.05);
 136         s->ray[i]->h.pulse_count = aray->n_pulses;
 137
 138         s->ray[i]->h.beam_width = 1.0; /* What is it really? */
 139         s->ray[i]->h.f = f;
 140         s->ray[i]->h.invf = invf;
 141
 142         ray_data = (unsigned char *)aray;
 143
 144         m = aray->offset[kk];
 145         if(m==0) continue;
 146
 147 /*  conversion changes by m.whimpey   */
 148         for(j=0; j<s->ray[i]->h.nbins; j++) {
 149           switch (kk) {
 150           case OFF_UZ: /* UZ field. */
 151           case OFF_CZ: /* CZ field. */
 152                 /* Apply a 1.4 dB correction. Ken Glasson did this. */
 153                 /* Removed 1.4 dB correction 09/27/2006--BMRC doesn't use it. */
 154                 if (period == BERRIMAH1)
 155                     x[j] = ((float)ray_data[m+j] - 56.0)/2.0; /* Removed +1.4 */
 156                 else
 157                     x[j] = ((float)ray_data[m+j] - 64.0)/2.0; /* Removed +1.4 */
 158                 break;
 159           case OFF_VEL: /* VR field */
 160                 if (period == BERRIMAH1)
 161                    x[j] =  (float)(Vu*((double)ray_data[m+j]-128.0)/128.);
 162                 else if (period == SCSMEX)
 163                   x[j] =  (float)(Vu*((double)(ray_data[m+j]^0x80)-128.0)/127.);
 164                 else
 165                    x[j] =  (float)(Vu*((double)ray_data[m+j]-128.0)/127.);
 166 /*              fprintf(stderr,"Velocity for ray[%d] at x[%d] = %f, nyquist = %f\n",i,j,x[j],s->ray[i]->h.nyq_vel); */
 167                 break;
 168           case OFF_WID: /* SW field */
 169                 if (period == BERRIMAH1)
 170                     x[j] = (float)(Vu*(double)ray_data[m+j]/100.);
 171                 else
 172                     x[j] = (float)(Vu*(double)ray_data[m+j]/256.);
 173                 break;
 174           case OFF_ZDR: /* ZD field */
 175                 if (period < MCTEX_EARLY) break;
 176                 if (period <= BERRIMAH2)
 177                     x[j] = ((float)ray_data[m+j] - 64.0)/21.25;
 178                 else
 179                     x[j] = ((float)ray_data[m+j] - 128.)*18./254.;
 180                 break;
 181           case OFF_PHI: /* PH field */
 182                 if (period < MCTEX_EARLY) break;
 183                 if (period <= BERRIMAH2)
 184                     x[j] = ((float)ray_data[m+j] - 128.0)*32.0/127.0;
 185                 else if (period == GUNN_PT)
 186                     x[j] = ((float)ray_data[m+j] - 64.5)*360.0/254.0;
 187                 else if (period > GUNN_PT)
 188                     x[j] = 90.0+((float)ray_data[m+j] - 128.0)*180.0/254.0;
 189                     /*Extra 90 degrees added by Scott Collis, CAWCR/BMRC 2008*/
 190                 break;
 191           case OFF_RHO: /* RH field */
 192                 if (period < MCTEX_EARLY) break;
 193                 if (period <= BERRIMAH2)
 194                     x[j] = sqrt((float)ray_data[m+j]/256.822 + 0.3108);
 195                 else
 196                     x[j] = (((float)ray_data[m+j]-1.)*1.14/254.) + 0.01;
 197                 break;
 198           case OFF_LDR: /* LR field -- no 'period' here. */
 199                 x[j] = ((float)ray_data[m+j] - 250.0)/6;
 200                 break;
 201           case OFF_KDP: /* KP field -- no 'period' here. */
 202                 x[j] = ((float)ray_data[m+j]); /* This conversion is taken
 203                                                                                   from SIGMET.  Is this right? */
 204                 break;
 205           case OFF_TIME: /* TIME field -- no 'period' here. */
 206                 x[j] = ((float)ray_data[m+j]); /* This conversion is taken
 207                                                                                   from SIGMET.  Is this right? */
 208                 break;
 209           }
 210           if (ray_data[m+j] == 0) x[j] = BADVAL;
 211         }
 212         for (j=0; j<s->ray[i]->h.nbins; j++)
 213           s->ray[i]->range[j] = invf(x[j]);
 214   }
 215 }
 216
 217
 218
 219 /**********************************************************************/
 220 /*                                                                    */
 221 /*                     RSL_lassen_to_radar                            */
 222 /*                                                                    */
 223 /*  By: John Merritt                                                  */
 224 /*      Space Applications Corporation                                */
 225 /*      May  26, 1994                                                 */
 226 /**********************************************************************/
 227 Radar *RSL_lassen_to_radar(char *infile)
 228 {
 229 /* Lassen specific. */
 230   Lassen_sweep *ptr;
 231   Lassen_ray *aray;
 232   int period;  /*   m.whimpey changed early variable to period  */
 233   FILE *f;
 234   int q[MAX_RADAR_VOLUMES];
 235   extern int rsl_qfield[];
 236   extern int *rsl_qsweep; /* See RSL_read_these_sweeps in volume.c */
 237   extern int rsl_qsweep_max;
 238
 239 /* Radar specific */
 240   Radar *radar;
 241   int i, j, k;
 242   /* Listing of the RSL field type indexes, in the order that
 243    * LASSEN stores them.
 244    */
 245   int rsl_index[] = {ZT_INDEX, DZ_INDEX, VR_INDEX, SW_INDEX, ZD_INDEX,
 246                                          PH_INDEX, RH_INDEX, LR_INDEX, KD_INDEX, TI_INDEX};
 247
 248   char *ltype[] = {"UZ", "CZ", "Vel", "Wid", "Zdr", "Phi",
 249                                    "Rho", "Ldr", "Kdp", "Time"};
 250
 251
 252   struct compare_date {   /*   introduced by m.whimpey   */
 253         int year;
 254         int month;
 255         int day;
 256         int hour;
 257         int minute;
 258         int second;
 259   };
 260 /*  M.Whimpey changes made 19990422 for extra periods */
 261 #define NUM_DATES 8
 262   struct compare_date  cvrt_date[NUM_DATES] =
 263   { {1992,  1,  1,  0,  0,  0},
 264         {1995, 11,  1,  0,  0,  0},
 265         {1995, 11, 25, 20,  5,  0},
 266         {1996,  1,  1,  0,  0,  0},
 267         {1997, 10,  1,  0,  0,  0},
 268         {1997, 11, 20,  3, 40,  0},
 269         {1998, 05, 04,  0,  0,  0},
 270         {1998, 05, 17, 03, 47,  0} };
 271
 272   int d;        /* date counter */
 273
 274   unsigned long vt;     /* vol time */
 275   unsigned long dt;     /* date time */
 276
 277
 278     /*   open Lassen file  */
 279   if (infile == NULL) {
 280         int save_fd;
 281         save_fd = dup(0);
 282         f = fdopen(save_fd, "r");
 283   }  else
 284     if((f=fopen(infile, "r"))==(FILE *)NULL) {
 285           perror(infile);
 286           return NULL;
 287     }
 288   f = uncompress_pipe(f); /* Transparently, use gunzip. */
 289 #define NEW_BUFSIZ 16384
 290   setvbuf(f,NULL,_IOFBF,(size_t)NEW_BUFSIZ); /* Faster i/o? */
 291
 292     if((read_entire_lassen_file(f, &vol)) == 0)
 293     {
 294         perror("RSL_lassen_to_radar ... read_entire_lassen_file");
 295         exit(1);
 296     }
 297
 298         rsl_pclose(f);
 299
 300         if (radar_verbose_flag) {
 301           fprintf(stderr,"\n Version   = %d",vol.version);
 302           fprintf(stderr,"\n Volume    = %d",vol.volume);
 303           fprintf(stderr,"\n Numsweeps = %d",vol.numsweeps);
 304           fprintf(stderr,"\n Time      = %2.2d/%2.2d/%2.2d %2.2d:%2.2d:%2.2d - %2.2d:%2.2d:%2.2d",
 305                  (int)vol.month, (int)vol.day, (int)vol.year,
 306                  (int)vol.shour, (int)vol.sminute, (int)vol.ssecond,
 307                  (int)vol.ehour, (int)vol.eminute, (int)vol.esecond);
 308           fprintf(stderr,"\n Angle: start %d, stop %d", (int)vol.a_start, (int)vol.a_stop);
 309           fprintf(stderr,"\n");
 310         }
 311
 312 /*  determine which period the lassen volume belongs   m.whimpey   */
 313         vt = (vol.year-90) * 32140800;
 314         vt += vol.month * 2678400;
 315         vt += vol.day * 86400;
 316         vt += vol.shour * 3600;
 317         vt += vol.sminute * 60;
 318         vt += vol.ssecond;
 319
 320         for(d=0; d<NUM_DATES; d++) {
 321             dt = (cvrt_date[d].year-1990) * 32140800;
 322             dt += cvrt_date[d].month * 2678400;
 323             dt += cvrt_date[d].day * 86400;
 324             dt += cvrt_date[d].hour * 3600;
 325             dt += cvrt_date[d].minute * 60;
 326             dt += cvrt_date[d].second;
 327
 328             if(vt<=dt) {
 329                 break;
 330             }
 331         }
 332
 333         if(d==0) {
 334             fprintf(stderr, "%s: Error Vol date before first known!\n",
 335                      infile);
 336             exit(3);
 337         }
 338         period = d-1;
 339
 340 /* Max. expected volumes. */
 341   radar = RSL_new_radar(MAX_RADAR_VOLUMES);
 342
 343   radar->h.month      = vol.month;
 344   radar->h.day        = vol.day;
 345   radar->h.year       = vol.year + 1900;
 346   radar->h.hour       = vol.shour;
 347   radar->h.minute     = vol.sminute;
 348   radar->h.sec        = vol.ssecond;
 349   strcpy(radar->h.radar_type, "lassen");
 350   radar->h.nvolumes   = MAX_RADAR_VOLUMES;
 351   memcpy(&radar->h.radar_name, vol.radinfo.radar_name, 8);
 352   memcpy(&radar->h.name, vol.radinfo.site_name, 8);
 353   memcpy(&radar->h.city, "????", 4);
 354   memcpy(&radar->h.state,"AU", 2);
 355   radar->h.latd       = vol.radinfo.latitude.degree;
 356   radar->h.latm       = vol.radinfo.latitude.minute;
 357   radar->h.lats       = vol.radinfo.latitude.second;
 358   /* Is there a problem with the minutes/seconds when negative?
 359    * The degree/minute/sec all should have the same sign.
 360    */
 361   if (radar->h.latd < 0) {
 362         if (radar->h.latm > 0) radar->h.latm *= -1;
 363         if (radar->h.lats > 0) radar->h.lats *= -1;
 364   }
 365   radar->h.lond       = vol.radinfo.longitude.degree;
 366   radar->h.lonm       = vol.radinfo.longitude.minute;
 367   radar->h.lons       = vol.radinfo.longitude.second;
 368   if (radar->h.lond < 0) {
 369         if (radar->h.lonm > 0) radar->h.lonm *= -1;
 370         if (radar->h.lons > 0) radar->h.lons *= -1;
 371   }
 372   radar->h.height     = vol.radinfo.antenna_height;
 373   radar->h.spulse     = 0;
 374   radar->h.lpulse     = 0;
 375
 376
 377   /*   iterate for each sweep in radar volume   */
 378
 379 /* Determine which field types exist. The array 'q' is a boolean to
 380  * force one print messages, if requested.
 381  *
 382  * Well, I tried looping through all the rays and examining aray->flags.<?>,
 383  * but, it turns out that these flags bounce around, ie. toggle, for
 384  * fields that don't really exist.  Therefore, I cannot just logically OR
 385  * the field flags together to determine which fields exist -- is this
 386  * a bug when the file was created?
 387  *
 388  * What I've seen in lass2uf is to examine the first ray of the first sweep,
 389  * but, I think that is unreliable too, due to the above finding.
 390  *
 391  * The solution, now, is to examine the existance of OFFSET values.
 392  * This seems to be consistant, throughout the volume.
 393  * The OFFSET is really what is used, anyway, to extract the data.
 394  */
 395   memset(q, 0, sizeof(q));
 396   for (i=0; i<vol.numsweeps; i++) {
 397         ptr = vol.index[i];
 398         for (j=0; j<ptr->numrays; j++) {
 399           aray = ptr->ray[j];
 400           for (k=0; k<NUMOFFSETS; k++) {
 401                 if (aray->offset[k] != 0 && !q[rsl_index[k]]) {
 402                   /* From RSL_select_fields */
 403                   if (rsl_qfield[rsl_index[k]] == 1)
 404                         q[rsl_index[k]]=1;
 405                 }
 406           }
 407         }
 408   }
 409   if (radar_verbose_flag)
 410         fprintf(stderr,"\n Fields are (Lassen nomenclature):");
 411   for (k=0; k<NUMOFFSETS; k++) {
 412         i = rsl_index[k]; /* Lassen index order to RSL index order translation. */
 413         if (q[i])  {
 414           if (radar_verbose_flag) fprintf(stderr," %s", ltype[k]);
 415         }
 416   }
 417   if (radar_verbose_flag)   fprintf(stderr,"\n");
 418
 419   if (radar_verbose_flag) fprintf(stderr," Fields are    (RSL nomenclature):");
 420   for (k=0; k<NUMOFFSETS; k++) {
 421
 422         /* BTW, it doesn't matter if we allocate volumes, sweeps, or rays that
 423          * have no data.  The radar library will check for NULL pointers.
 424          */
 425         i = rsl_index[k]; /* Lassen index order to RSL index order translation. */
 426         if (q[i])  {
 427           radar->v[i] = RSL_new_volume(vol.numsweeps);
 428           radar->v[i]->h.f    = RSL_f_list[i];
 429           radar->v[i]->h.invf = RSL_invf_list[i];
 430           if (radar_verbose_flag)  fprintf(stderr," %s", RSL_ftype[i]);
 431           if (k >= 2 && radar_verbose_flag) fprintf(stderr," "); /* Alignment. */
 432         }
 433   }
 434   if (radar_verbose_flag)   fprintf(stderr,"\n");
 435
 436   for(j=0;j<radar->h.nvolumes; j++) {
 437         for(i=0;i<(int)vol.numsweeps;i++) {
 438           if (rsl_qsweep != NULL) {
 439                 if (i > rsl_qsweep_max) break;
 440                 if (rsl_qsweep[i] == 0) continue;
 441           }
 442           ptr = vol.index[i];
 443           if (radar->v[j]) {
 444                 radar->v[j]->sweep[i] = RSL_new_sweep(ptr->numrays);
 445
 446                 /* 'period' is a flag for different calibrations */
 447                 lassen_load_sweep(radar->v[j]->sweep[i], i, j, period, ptr);
 448           }
 449         }
 450   }
 451   radar = RSL_prune_radar(radar);
 452   return radar;
 453 }
 454 #else
 455 Radar *RSL_lassen_to_radar(char *infile)
 456 {
 457   fprintf(stderr, "LASSEN is not installed in this version of RSL.\n");
 458   fprintf(stderr, "Reinstall RSL w/ -DHAVE_LASSEN in the Makefile.\n");
 459   return NULL;
 460 }
 461 #endif