RSL v1.44
[~andy/rsl] / wsr88d_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 #include <stdio.h>
24 #include <string.h>
25
26 #include <stdlib.h>
27
28 #include "rsl.h"
29 #include "wsr88d.h"
30
31 extern int radar_verbose_flag;
32 /*
33  * These externals can be found in the wsr88d library; secret code.
34  */
35 void print_head(Wsr88d_file_header wsr88d_file_header);
36 void clear_sweep(Wsr88d_sweep *wsr88d_sweep, int x, int n);
37 void free_and_clear_sweep(Wsr88d_sweep *wsr88d_sweep, int x, int n);
38 /*
39  * Secretly in uf_to_radar.c
40  */
41 Volume *copy_sweeps_into_volume(Volume *new_volume, Volume *old_volume);
42
43  
44 void float_to_range(float *x, Range *c, int n, Range (*function)(float x) )
45 {
46   while (n--) {
47     if (*x == WSR88D_BADVAL) *c = function(BADVAL);
48     else if (*x == WSR88D_RFVAL) *c = function(RFVAL);
49     else *c = function(*x);
50     c++; x++;
51   }
52 }
53
54 /**********************************************************************/
55 /*                                                                    */
56 /* done 3/23         wsr88d_load_sweep_into_volume                    */
57 /*                                                                    */
58 /*  By: John Merritt                                                  */
59 /*      Space Applications Corporation                                */
60 /*      March 3, 1994                                                 */
61 /**********************************************************************/
62 int wsr88d_load_sweep_into_volume(Wsr88d_sweep ws,
63                        Volume *v, int nsweep, unsigned int vmask)
64 {
65   int i;
66   int iray;
67   float v_data[1000];
68   Range  c_data[1000];
69   int n;
70
71   int mon, day, year;
72   int hh, mm, ss;
73   float fsec;
74   int vol_cpat;
75
76   Ray *ray_ptr;
77   Range (*invf)(float x);
78   float (*f)(Range x);
79
80   /* Allocate memory for MAX_RAYS_IN_SWEEP rays. */
81   v->sweep[nsweep] = RSL_new_sweep(MAX_RAYS_IN_SWEEP);
82   if (v->sweep[nsweep] == NULL) {
83     perror("wsr88d_load_sweep_into_volume: RSL_new_sweep");
84     return -1;
85   }
86
87   v->sweep[nsweep]->h.elev = 0;
88   v->sweep[nsweep]->h.nrays = 0;
89   f = (float (*)(Range x))NULL;
90   invf = (Range (*)(float x))NULL;
91   if (vmask & WSR88D_DZ) { invf = DZ_INVF; f = DZ_F; }
92   if (vmask & WSR88D_VR) { invf = VR_INVF; f = VR_F; }
93   if (vmask & WSR88D_SW) { invf = SW_INVF; f = SW_F; }
94   
95   v->h.invf = invf;
96   v->h.f    = f;
97   v->sweep[nsweep]->h.invf = invf;
98   v->sweep[nsweep]->h.f    = f;
99
100   for (i=0,iray=0; i<MAX_RAYS_IN_SWEEP; i++) {
101     if (ws.ray[i] != NULL) {
102       wsr88d_ray_to_float(ws.ray[i], vmask, v_data, &n);
103       float_to_range(v_data, c_data, n, invf);
104       if (n > 0) {
105         wsr88d_get_date(ws.ray[i], &mon, &day, &year);
106         wsr88d_get_time(ws.ray[i], &hh, &mm, &ss, &fsec);
107         /*
108       fprintf(stderr,"n %d, mon %d, day %d, year %d,  hour %d, min %d, sec %d, fsec %f\n",
109             n, mon, day, year, hh, mm, ss, fsec);
110             */
111         /*
112          * Load the sweep/ray headar information.
113          */
114         
115         v->sweep[nsweep]->ray[iray] = RSL_new_ray(n);
116         /*(Range *)calloc(n, sizeof(Range)); */
117         
118         ray_ptr = v->sweep[nsweep]->ray[iray]; /* Make code below readable. */
119         ray_ptr->h.f        = f;
120         ray_ptr->h.invf     = invf;
121         ray_ptr->h.month    = mon;
122         ray_ptr->h.day      = day;
123         ray_ptr->h.year     = year + 1900; /* Yes 1900 makes this year 2000 compliant, due to wsr88d using unix time(). */
124         ray_ptr->h.hour     = hh;
125         ray_ptr->h.minute   = mm;
126         ray_ptr->h.sec      = ss + fsec;
127         ray_ptr->h.unam_rng = wsr88d_get_range   (ws.ray[i]);
128         ray_ptr->h.azimuth  = wsr88d_get_azimuth (ws.ray[i]);
129 /* -180 to +180 is converted to 0 to 360 */
130         if (ray_ptr->h.azimuth < 0) ray_ptr->h.azimuth += 360;
131         ray_ptr->h.ray_num  = ws.ray[i]->ray_num;
132         ray_ptr->h.elev       = wsr88d_get_elevation_angle(ws.ray[i]);
133         ray_ptr->h.elev_num   = ws.ray[i]->elev_num;
134         if (vmask & WSR88D_DZ) {
135           ray_ptr->h.range_bin1 = ws.ray[i]->refl_rng;
136           ray_ptr->h.gate_size  = ws.ray[i]->refl_size;
137         } else {
138           ray_ptr->h.range_bin1 = ws.ray[i]->dop_rng;
139           ray_ptr->h.gate_size  = ws.ray[i]->dop_size;
140         }
141         ray_ptr->h.vel_res  = wsr88d_get_velocity_resolution(ws.ray[i]);
142         vol_cpat = wsr88d_get_volume_coverage(ws.ray[i]);
143         switch (vol_cpat) {
144         case 11: ray_ptr->h.sweep_rate = 16.0/5.0;  break;
145         case 12: ray_ptr->h.sweep_rate = 17.0/4.2;  break;
146         case 21: ray_ptr->h.sweep_rate = 11.0/6.0;  break;
147         case 31: ray_ptr->h.sweep_rate =  8.0/10.0; break;
148         case 32: ray_ptr->h.sweep_rate =  7.0/10.0; break;
149         case 121:ray_ptr->h.sweep_rate = 20.0/5.5;  break;
150         default: ray_ptr->h.sweep_rate =  0.0; break;
151         }
152
153         ray_ptr->h.nyq_vel  = wsr88d_get_nyquist(ws.ray[i]);
154         ray_ptr->h.azim_rate   = wsr88d_get_azimuth_rate(ws.ray[i]);
155         ray_ptr->h.fix_angle   = wsr88d_get_fix_angle(ws.ray[i]);
156         ray_ptr->h.pulse_count = wsr88d_get_pulse_count(ws.ray[i]);
157         ray_ptr->h.pulse_width = wsr88d_get_pulse_width(ws.ray[i]);
158         ray_ptr->h.beam_width  = .95;
159         ray_ptr->h.prf         = wsr88d_get_prf(ws.ray[i]);
160         ray_ptr->h.frequency   = wsr88d_get_frequency(ws.ray[i]);
161         ray_ptr->h.wavelength = 0.1071; /* Previously called
162                                      * wsr88d_get_wavelength(ws.ray[i]).
163                                      * See wsr88d.c for explanation.
164                          */
165
166         /* It is no coincidence that the 'vmask' and wsr88d datatype
167          * values are the same.  We expect 'vmask' to be one of
168          * REFL_MASK, VEL_MASK, or SW_MASK.  These match WSR88D_DZ,
169          * WSR88D_VR, and WSR88D_SW in the wsr88d library.
170          */
171         ray_ptr->h.nbins = n;
172         memcpy(ray_ptr->range, c_data, n*sizeof(Range));
173         v->sweep[nsweep]->h.nrays = iray+1;
174         v->sweep[nsweep]->h.elev += ray_ptr->h.elev;
175         v->sweep[nsweep]->h.sweep_num = ray_ptr->h.elev_num;
176         iray++;
177       }
178     }
179   }
180   v->sweep[nsweep]->h.beam_width = .95;
181   v->sweep[nsweep]->h.vert_half_bw = .475;
182   v->sweep[nsweep]->h.horz_half_bw = .475;
183   /* Now calculate the mean elevation angle for this sweep. */
184   if (v->sweep[nsweep]->h.nrays > 0)
185     v->sweep[nsweep]->h.elev /= v->sweep[nsweep]->h.nrays;
186   else {
187     RSL_free_sweep(v->sweep[nsweep]); /* No rays loaded, free this sweep. */
188     v->sweep[nsweep] = NULL;
189   }
190   
191   return 0;
192 }
193
194 /**********************************************************************/
195 /*                                                                    */
196 /*                     RSL_wsr88d_to_radar                            */
197 /*                                                                    */
198 /*  By: John Merritt                                                  */
199 /*      Space Applications Corporation                                */
200 /*      March 3, 1994                                                 */
201 /**********************************************************************/
202
203 Radar *RSL_wsr88d_to_radar(char *infile, char *call_or_first_tape_file)
204 /*
205  * Gets all volumes from the nexrad file.  Input file is 'infile'.
206  * Site information is extracted from 'call_or_first_tape_file'; this
207  * is typically a disk file called 'nex.file.1'.
208  *
209  *  -or-
210  *
211  * Uses the string in 'call_or_first_tape_file' as the 4 character call sign
212  * for the sight.  All UPPERCASE characters.  Normally, this call sign
213  * is extracted from the file 'nex.file.1'.
214  *
215  * Returns a pointer to a Radar structure; that contains the different
216  * Volumes of data.
217  */
218 {
219   Radar *radar;
220   Volume *new_volume;
221   Wsr88d_file *wf;
222   Wsr88d_sweep wsr88d_sweep;
223   Wsr88d_file_header wsr88d_file_header;
224   Wsr88d_tape_header wsr88d_tape_header;
225   int n;
226   int nsweep;
227   int i;
228   int iv;
229   int nvolumes;
230   int volume_mask[] = {WSR88D_DZ, WSR88D_VR, WSR88D_SW};
231   char *field_str[] = {"Reflectivity", "Velocity", "Spectrum width"};
232   Wsr88d_site_info *sitep;
233   char site_id_str[5];
234   char *the_file;
235   int expected_msgtype = 0;
236   char version[8];
237   int vnum;
238
239   extern int rsl_qfield[]; /* See RSL_select_fields in volume.c */
240   extern int *rsl_qsweep; /* See RSL_read_these_sweeps in volume.c */
241   extern int rsl_qsweep_max;
242
243   Radar *wsr88d_load_m31_into_radar(Wsr88d_file *wf);
244
245   sitep = NULL;
246 /* Determine the site quasi automatically.  Here is the procedure:
247  *    1. Determine if we have a call sign.
248  *    2. Try reading 'call_or_first_tape_file' from disk.  This is done via
249  *       wsr88d_read_tape_header.
250  *    3. If no valid site info, abort.
251  */
252   if (call_or_first_tape_file == NULL) {
253     fprintf(stderr, "wsr88d_to_radar: No valid site ID info provided.\n");
254     return(NULL);
255   } else if (strlen(call_or_first_tape_file) == 4)
256     sitep =  wsr88d_get_site(call_or_first_tape_file);
257   else if (strlen(call_or_first_tape_file) == 0) {
258     fprintf(stderr, "wsr88d_to_radar: No valid site ID info provided.\n");
259     return(NULL);
260   }  
261
262   if (sitep == NULL)
263     if (wsr88d_read_tape_header(call_or_first_tape_file, &wsr88d_tape_header) > 0) {
264       memcpy(site_id_str, wsr88d_tape_header.site_id, 4);
265       sitep  = wsr88d_get_site(site_id_str);
266     }
267   if (sitep == NULL) {
268       fprintf(stderr,"wsr88d_to_radar: No valid site ID info found.\n");
269         return(NULL);
270   }
271     if (radar_verbose_flag)
272       fprintf(stderr,"SITE: %c%c%c%c\n", sitep->name[0], sitep->name[1],
273              sitep->name[2], sitep->name[3]);
274
275   
276   memset(&wsr88d_sweep, 0, sizeof(Wsr88d_sweep)); /* Initialize to 0 a 
277                                                    * heavily used variable.
278                                                    */
279
280 /* 1. Open the input wsr88d file. */
281   if (infile == NULL) the_file = "stdin";  /* wsr88d.c understands this to
282                                             * mean read from stdin.
283                                             */
284   else the_file = infile;
285
286   if ((wf = wsr88d_open(the_file)) == NULL) {
287     wsr88d_perror(the_file);
288     return NULL;
289   }
290
291
292
293 /* 2. Read wsr88d headers. */
294   /* Return # bytes, 0 or neg. on fail. */
295   n = wsr88d_read_file_header(wf, &wsr88d_file_header);
296   /*
297    * Get the expected digital radar message type based on version string
298    * from the Archive II header.  The message type is 31 for Build 10, and 1
299    * for prior builds.
300    */
301   if (n > 0) {
302       strncpy(version, wsr88d_file_header.title.filename, 8);
303       if (strncmp(version,"AR2V",4) == 0) {
304           sscanf(version, "AR2V%4d", &vnum);
305           if (vnum > 1) expected_msgtype = 31;
306           else expected_msgtype = 1;
307       }
308       else if (strncmp(version,"ARCHIVE2",8) == 0) expected_msgtype = 1;
309   }
310
311   if (n <= 0 || expected_msgtype == 0) {
312       fprintf(stderr,"RSL_wsr88d_to_radar: ");
313       if (n <= 0)
314       fprintf(stderr,"wsr88d_read_file_header failed\n");
315       else
316       fprintf(stderr,"Archive II header contains unknown version "
317           ": '%s'\n", version);
318       wsr88d_close(wf);
319       return NULL;
320   }
321
322   if (radar_verbose_flag)
323     print_head(wsr88d_file_header);
324
325
326   if (expected_msgtype == 31) {
327       /* Get radar for message type 31. */
328       radar = wsr88d_load_m31_into_radar(wf);
329       if (radar == NULL) return NULL;
330   }
331   else {
332       /* Get radar for message type 1. */
333       nvolumes = 3;
334       /* Allocate all Volume pointers. */
335       radar = RSL_new_radar(MAX_RADAR_VOLUMES);
336       if (radar == NULL) return NULL;
337
338     /* Clear the sweep pointers. */
339       clear_sweep(&wsr88d_sweep, 0, MAX_RAYS_IN_SWEEP);
340
341     /* Allocate a maximum of 30 sweeps for the volume. */
342     /* Order is important.  WSR88D_DZ, WSR88D_VR, WSR88D_SW, is 
343      * assigned to the indexes DZ_INDEX, VR_INDEX and SW_INDEX respectively.
344      */ 
345
346       for (iv=0; iv<nvolumes; iv++)
347         if (rsl_qfield[iv]) radar->v[iv] = RSL_new_volume(20);
348
349
350     /* LOOP until EOF */
351       nsweep = 0;
352       for (;(n = wsr88d_read_sweep(wf, &wsr88d_sweep)) > 0; nsweep++) {
353         if (rsl_qsweep != NULL) {
354           if (nsweep > rsl_qsweep_max) break;
355           if (rsl_qsweep[nsweep] == 0) continue;
356         }
357         if (radar_verbose_flag)  
358         fprintf(stderr,"Processing for SWEEP # %d\n", nsweep);
359
360           /*  wsr88d_print_sweep_info(&wsr88d_sweep); */
361         
362         for (iv=0; iv<nvolumes; iv++) {
363           if (rsl_qfield[iv]) {
364             /* Exceeded sweep limit.
365              * Allocate more sweeps.
366              * Copy all previous sweeps.
367              */
368             if (nsweep >= radar->v[iv]->h.nsweeps) {
369               if (radar_verbose_flag)
370                 fprintf(stderr,"Exceeded sweep allocation of %d. "
371                     "Adding 20 more.\n", nsweep);
372               new_volume = RSL_new_volume(radar->v[iv]->h.nsweeps+20);
373               new_volume = copy_sweeps_into_volume(new_volume, radar->v[iv]);
374               radar->v[iv] = new_volume;
375             }
376             if (wsr88d_load_sweep_into_volume(wsr88d_sweep,
377                radar->v[iv], nsweep, volume_mask[iv]) != 0) {
378               RSL_free_radar(radar);
379               return NULL;
380             }
381           }
382         }
383         if (nsweep == 0) {
384           /* Get Volume Coverage Pattern number for radar header. */
385           i=0;
386           while (i < MAX_RAYS_IN_SWEEP && wsr88d_sweep.ray[i] == NULL) i++;
387           if (i < MAX_RAYS_IN_SWEEP) radar->h.vcp = wsr88d_get_volume_coverage(
388             wsr88d_sweep.ray[i]);
389         }
390
391         free_and_clear_sweep(&wsr88d_sweep, 0, MAX_RAYS_IN_SWEEP);
392       }
393
394       for (iv=0; iv<nvolumes; iv++) {
395         if (rsl_qfield[iv]) {
396           radar->v[iv]->h.type_str = strdup(field_str[iv]);
397           radar->v[iv]->h.nsweeps = nsweep;
398         }
399       }
400   }
401   wsr88d_close(wf);
402
403 /*
404  * Here we will assign the Radar_header information.  Take most of it
405  * from an existing volume's header.  
406  */
407   radar_load_date_time(radar);  /* Magic :-) */
408
409     radar->h.number = sitep->number;
410     memcpy(&radar->h.name, sitep->name, sizeof(sitep->name));
411     memcpy(&radar->h.radar_name, sitep->name, sizeof(sitep->name)); /* Redundant */
412     memcpy(&radar->h.city, sitep->city, sizeof(sitep->city));
413     memcpy(&radar->h.state, sitep->state, sizeof(sitep->state));
414     strcpy(radar->h.radar_type, "wsr88d");
415     radar->h.latd = sitep->latd;
416     radar->h.latm = sitep->latm;
417     radar->h.lats = sitep->lats;
418     if (radar->h.latd < 0) { /* Degree/min/sec  all the same sign */
419       radar->h.latm *= -1;
420       radar->h.lats *= -1;
421     }
422     radar->h.lond = sitep->lond;
423     radar->h.lonm = sitep->lonm;
424     radar->h.lons = sitep->lons;
425     if (radar->h.lond < 0) { /* Degree/min/sec  all the same sign */
426       radar->h.lonm *= -1;
427       radar->h.lons *= -1;
428     }
429     radar->h.height = sitep->height;
430     radar->h.spulse = sitep->spulse;
431     radar->h.lpulse = sitep->lpulse;
432
433     free(sitep);
434
435   if (wsr88d_merge_split_cuts_is_set()) radar = wsr88d_merge_split_cuts(radar);
436   return radar;
437 }