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