Pileus Git - ~andy/rsl/blob - src/rainbow_to_radar.c

   1 /*
   2     NASA/TRMM, Code 912
   3     This is the TRMM Office Radar Software Library.
   4     Copyright (C) 2004
   5             Bart Kelley
   6             George Mason University
   7             Fairfax, 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 Software
  21     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  22 */
  23
  24 #include <stdio.h>
  25 #include <stdlib.h>
  26 #include <string.h>
  27 #include <unistd.h>
  28 #include "rsl.h"
  29 #include "rainbow.h"
  30
  31 struct dms {
  32     int deg;
  33     int minute;
  34     int sec;
  35 };
  36
  37 struct dms deg_to_dms(float deg)
  38 {
  39     /* Convert value in degrees to degrees-minutes-seconds. */
  40
  41     int sign, ideg, iminute, sec;
  42     float fminute;  /* minute with fractional part */
  43     struct dms dms;
  44
  45     sign = (deg < 0.) ? -1 : 1;
  46     deg = sign * deg;  /* absolute value */
  47     ideg = deg;        /* truncated degree */
  48     iminute = fminute = (deg - ideg) * 60.;
  49     sec = (fminute - iminute) * 60. + .5;
  50     dms.deg = sign * ideg;
  51     dms.minute = sign * iminute;
  52     dms.sec = sign * sec;
  53     return dms;
  54 }
  55
  56 static float (*f)(Range x);
  57 static Range (*invf)(float x);
  58
  59 /**********************************************************/
  60 /*                                                        */
  61 /*                  RSL_rainbow_to_radar                  */
  62 /*                                                        */
  63 /**********************************************************/
  64
  65 Radar *RSL_rainbow_to_radar(char *infile)
  66 {
  67     /* This function reads the Rainbow format scan data file and returns a
  68      * radar structure.
  69      */
  70
  71     Radar *radar;
  72     FILE *fp;
  73     int c;
  74     int nvolumes;
  75     Rainbow_hdr rainbow_hdr;
  76     struct dms latdms, londms;
  77
  78     void read_rainbow_header(Rainbow_hdr *, FILE *);
  79
  80     /* These next lines allow program to read from a regular file, a
  81      * compressed file, or standard input.  I lifted them from RSL_uf_to_radar
  82      * by John Merritt.
  83      */
  84
  85     if (infile == NULL) {
  86         int save_fd;
  87         save_fd = dup(0);
  88         fp = fdopen(save_fd, "rb");
  89     }
  90     else if ((fp = fopen(infile, "rb")) == NULL) {
  91         perror(infile);
  92         return NULL;
  93     }
  94     fp = uncompress_pipe(fp); /* Transparently gunzip. */
  95
  96     /* Read first character and verify file format. */
  97
  98     if ((c = fgetc(fp)) != SOH) {
  99         fprintf(stderr,"%s is not a valid Rainbow format file.\n",infile);
 100         return NULL;
 101     }
 102
 103     /* Read Rainbow file header and check for correct product. */
 104
 105     read_rainbow_header(&rainbow_hdr, fp);
 106
 107     if (rainbow_hdr.filetype != SCAN_DATA ) {
 108         fprintf(stderr,"ERROR: File is not a scan data file.\n");
 109         fprintf(stderr,"File type number (header label H3) is %d\n",
 110                 rainbow_hdr.filetype);
 111         fprintf(stderr,"See Rainbow File Format Document for details on "
 112                 "header labels.\n");
 113         return NULL;
 114     }
 115     if (rainbow_hdr.product != VOLUME_SCAN) {
 116         fprintf(stderr,"WARNING: Product is not volume scan as expected.\n");
 117         fprintf(stderr,"Header label N is %d, expected %d\n",
 118                 rainbow_hdr.product, VOLUME_SCAN);
 119         fprintf(stderr,"See Rainbow File Format Document for details on "
 120                 "header labels.\n");
 121     }
 122     if (rainbow_hdr.compressed) {
 123         fprintf(stderr,"RSL_rainbow_to_radar: Label F3 indicates data "
 124                 "compression.\n");
 125         fprintf(stderr,"This routine can not handle compressed data.\n");
 126         fprintf(stderr,"See Rainbow File Format Document for details on "
 127                 "header labels.\n");
 128         return NULL;
 129     }
 130
 131     /* Make radar structure and assign header information.  */
 132
 133     nvolumes = 1; /* There is only one raw data type in a volume scan file. */
 134     radar = RSL_new_radar(MAX_RADAR_VOLUMES);
 135     if (radar == NULL) {
 136         perror("RSL_new_radar returned NUL to RSL_rainbow_to_radar");
 137         return NULL;
 138     }
 139     radar->h.nvolumes = nvolumes;
 140     strcpy(radar->h.radar_name, rainbow_hdr.radarname);
 141     strcpy(radar->h.radar_type, "rainbow");
 142     radar->h.month = rainbow_hdr.month;
 143     radar->h.day = rainbow_hdr.day;
 144     radar->h.year = rainbow_hdr.year;
 145     radar->h.hour = rainbow_hdr.hour;
 146     radar->h.minute = rainbow_hdr.minute;
 147     radar->h.sec = (float) rainbow_hdr.sec;
 148     latdms = deg_to_dms(rainbow_hdr.lat);
 149     londms = deg_to_dms(rainbow_hdr.lon);
 150     radar->h.latd = latdms.deg;
 151     radar->h.latm = latdms.minute;
 152     radar->h.lats = latdms.sec;
 153     radar->h.lond = londms.deg;
 154     radar->h.lonm = londms.minute;
 155     radar->h.lons = londms.sec;
 156
 157     /* Read the data portion of the file into the radar structure. */
 158
 159     if (rainbow_data_to_radar(radar, rainbow_hdr, fp) < 1)
 160         radar = NULL;
 161
 162     return radar;
 163 }
 164
 165 /**********************************************************/
 166 /*                                                        */
 167 /*                 rainbow_data_to_radar                  */
 168 /*                                                        */
 169 /**********************************************************/
 170
 171 int rainbow_data_to_radar(Radar *radar, Rainbow_hdr rainbow_hdr, FILE *fp)
 172 {
 173     /* Read Rainbow data into the radar structure.  Data in the file is stored
 174      * as bytes, where each byte contains the value for one range bin.  The
 175      * data is ordered as follows: The first byte corresponds to the first bin
 176      * of the first ray, followed by the remaining bins of that ray, followed
 177      * by the remaining rays of the first sweep.  This sequence of bins, rays,
 178      * and sweeps continues for the remainder of sweeps in the volume scan.
 179      */
 180
 181     int iray, isweep, nread, nsweeps, nrays, nbins;
 182     unsigned char *rainbow_ray;
 183     Volume *v;
 184     Sweep *sweep;
 185     Ray *ray;
 186     int i;
 187     float azim_rate, beam_width, dz, elev_angle, prf, unam_rng;
 188     int vol_index = VR_INDEX;
 189
 190     beam_width = 1.0; /* for now */
 191
 192     switch (rainbow_hdr.datatype) {
 193         /* TODO: Add f and invf for each field */
 194         case 0:
 195             vol_index = VR_INDEX;
 196             f = VR_F;
 197             invf = VR_INVF;
 198             break;
 199         case 1:
 200             vol_index = DZ_INDEX;
 201             f = DZ_F;
 202             invf = DZ_INVF;
 203             break;
 204         case 2:
 205             vol_index = SW_INDEX;
 206             break;
 207         case 3:
 208             vol_index = ZT_INDEX;
 209             break;
 210         case 5:
 211             vol_index = DR_INDEX;
 212             break;
 213     }
 214     /* TODO: remove the following if-statement once the code for other data
 215      * types has been implemented.  Currently only handles reflectivity (DZ).
 216      */
 217     if (vol_index != DZ_INDEX) {
 218         fprintf(stderr, "RSL_rainbow_to_radar: currently only handles "
 219                 "field type DZ\n");
 220         fprintf(stderr,"Rainbow data type value (label F9) is %d\n",
 221                 rainbow_hdr.datatype);
 222         fprintf(stderr,"Corresponding vol_INDEX number is %d\n", vol_index);
 223         return 0;
 224     }
 225
 226     nsweeps = rainbow_hdr.nsweeps;
 227     nrays = (rainbow_hdr.az_stop - rainbow_hdr.az_start + 1) /
 228         rainbow_hdr.az_step + .5;
 229     nbins = rainbow_hdr.nbins;
 230     if (nrays != 360) {
 231         fprintf(stderr,"WARNING: number of rays computed is not the number "
 232                 "expected.\n");
 233         fprintf(stderr,"Computed = nrays: "
 234             "azstart = %d, az_stop = %d, az_step = %f\n",
 235             rainbow_hdr.az_start, rainbow_hdr.az_stop, rainbow_hdr.az_step);
 236         fprintf(stderr,"Expected 360\n");
 237     }
 238     radar->v[vol_index] = RSL_new_volume(nsweeps);
 239     v = radar->v[vol_index];
 240     v->h.nsweeps = nsweeps;
 241     v->h.f = f;
 242     v->h.invf = invf;
 243     if (vol_index == DZ_INDEX) v->h.type_str = strdup("Reflectivity");
 244     rainbow_ray = (unsigned char *) malloc(nbins);
 245
 246     /* Load sweeps. */
 247
 248     for (isweep = 0; isweep < nsweeps; isweep++) {
 249         sweep = RSL_new_sweep(nrays);
 250         prf = rainbow_hdr.elev_params[isweep]->prf_high;
 251         unam_rng = RSL_SPEED_OF_LIGHT / (2. * prf * 1000.);
 252         azim_rate = rainbow_hdr.elev_params[isweep]->az_rate;
 253         elev_angle = rainbow_hdr.elev_params[isweep]->elev_angle;
 254
 255         /* Load rays. */
 256
 257         for (iray = 0; iray < nrays; iray++) {
 258             nread = fread(rainbow_ray, 1, nbins, fp);
 259             if (nread != nbins) {
 260                 fprintf(stderr, "ERROR: Could not read enough bytes to fill "
 261                         "ray.\n");
 262                 fprintf(stderr, "Sweep = %d, ray = %d, number read = %d\n",
 263                         isweep, iray, nread);
 264                 return 0;
 265             }
 266             ray = RSL_new_ray(nbins);
 267
 268             /* TODO: Add code for other fields. */
 269
 270             if (vol_index == DZ_INDEX) {
 271                 for (i=0; i < ray->h.nbins; i++) {
 272                     dz = -31.5 + (rainbow_ray[i] - 1) * 0.5;
 273                     if (dz < -31.5)
 274                         ray->range[i] = invf(BADVAL);
 275                     else
 276                         ray->range[i] = invf(dz);
 277                 }
 278             }
 279             /* Load ray header. Note: the rainbow data file has only one time
 280              * stamp, which is the data acquisition time at the start of the
 281              * volume scan. For now, that is the time assigned to all rays.
 282              */
 283             ray->h.month = rainbow_hdr.month;
 284             ray->h.day = rainbow_hdr.day;
 285             ray->h.year = rainbow_hdr.year;
 286             ray->h.hour = rainbow_hdr.hour;
 287             ray->h.minute = rainbow_hdr.minute;
 288             ray->h.sec = rainbow_hdr.sec;
 289             ray->h.f = f;
 290             ray->h.invf = invf;
 291             ray->h.unam_rng = unam_rng;
 292             ray->h.prf = prf;
 293             ray->h.azim_rate = azim_rate;
 294             ray->h.elev = elev_angle;
 295             ray->h.elev_num = isweep + 1;
 296             ray->h.fix_angle = elev_angle;
 297             ray->h.azimuth = rainbow_hdr.az_start + iray * rainbow_hdr.az_step;
 298             ray->h.ray_num = iray + 1;
 299             ray->h.range_bin1 = rainbow_hdr.range_start;
 300             ray->h.gate_size = rainbow_hdr.range_step * 1000.;
 301             ray->h.beam_width = beam_width;
 302             sweep->ray[iray] = ray;
 303         } /* iray */
 304         sweep->h.sweep_num = isweep + 1;
 305         sweep->h.beam_width = beam_width;
 306         sweep->h.vert_half_bw = sweep->h.beam_width * 0.5;
 307         sweep->h.horz_half_bw = sweep->h.beam_width * 0.5;
 308         sweep->h.nrays = nrays;
 309         sweep->h.f = f;
 310         sweep->h.invf = invf;
 311         v->sweep[isweep] = sweep;
 312     } /* isweep */
 313     return 1;
 314 }