3 This is the TRMM Office Radar Software Library.
4 Copyright (C) 1996-1999
6 Space Applications Corporation
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.
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.
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.
32 extern int radar_verbose_flag;
34 /********************************************************************/
36 /* find_rsl_field_index */
38 /********************************************************************/
39 int find_rsl_field_index(char *dorade_field_name)
42 * Dorade: VE, DM, SW, DBZ, ZDR, PHI, RHO, LDR, DX, CH, AH, CV, AV
43 * RSL: VR, DM, SW, DZ, ZD, PH, RH, LR, *DX, *CH, *AH, *CV, *AV.
45 if (strncasecmp(dorade_field_name, "ve", 2) == 0) return VR_INDEX;
46 if (strncasecmp(dorade_field_name, "vr", 2) == 0) return VR_INDEX;
47 if (strncasecmp(dorade_field_name, "dm", 2) == 0) return DM_INDEX;
48 if (strncasecmp(dorade_field_name, "sw", 2) == 0) return SW_INDEX;
49 if (strncasecmp(dorade_field_name, "dbz", 3) == 0) return DZ_INDEX;
50 if (strncasecmp(dorade_field_name, "zdr", 3) == 0) return ZD_INDEX;
51 if (strncasecmp(dorade_field_name, "phi", 3) == 0) return PH_INDEX;
52 if (strncasecmp(dorade_field_name, "rho", 3) == 0) return RH_INDEX;
53 if (strncasecmp(dorade_field_name, "ldr", 3) == 0) return LR_INDEX;
54 if (strncasecmp(dorade_field_name, "dx", 2) == 0) return DX_INDEX;
55 if (strncasecmp(dorade_field_name, "ch", 2) == 0) return CH_INDEX;
56 if (strncasecmp(dorade_field_name, "ah", 2) == 0) return AH_INDEX;
57 if (strncasecmp(dorade_field_name, "cv", 2) == 0) return CV_INDEX;
58 if (strncasecmp(dorade_field_name, "av", 2) == 0) return AV_INDEX;
59 if (strncasecmp(dorade_field_name, "vs", 2) == 0) return VS_INDEX;
60 if (strncasecmp(dorade_field_name, "vl", 2) == 0) return VL_INDEX;
61 if (strncasecmp(dorade_field_name, "vg", 2) == 0) return VG_INDEX;
62 if (strncasecmp(dorade_field_name, "vt", 2) == 0) return VT_INDEX;
63 if (strncasecmp(dorade_field_name, "ncp", 2) == 0) return NP_INDEX;
68 void prt_skipped_field_msg(char *dorade_field_name)
71 int i, already_printed;
73 static int nskipped = 0;
74 static char skipped_list[MAXFIELDS][9];
76 /* Make sure name is a properly formed string. */
77 strncpy(prtname, dorade_field_name, 8);
80 /* We don't want to repeat messages for the same field, so keep a list of
81 * fields already printed.
85 while (!already_printed && i < nskipped) {
86 if (strncmp(prtname, skipped_list[i], 2) == 0) already_printed = 1;
89 if (!already_printed) {
90 fprintf(stderr, "Unknown DORADE parameter type <%s> -- skipping.\n", prtname);
91 strcpy(skipped_list[nskipped], prtname);
96 /* For date conversion routines. */
97 static int daytab[2][13] = {
98 {0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334, 365},
99 {0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335, 366}
102 /*************************************************************/
106 /*************************************************************/
107 static int julian(int year, int mo, int day)
109 /* Converts a calendar date (month, day, year) to a Julian date.
115 leap = (year%4 == 0 && year%100 != 0) || year%400 == 0;
116 return(day + daytab[leap][mo-1]);
119 /*************************************************************/
123 /*************************************************************/
124 static void ymd(int jday, int year, int *mm, int *dd)
126 /* Input: jday, yyyy */
128 /* Copied from hdf_to_radar.c, written by Mike Kolander. */
133 leap = (year%4 == 0 && year%100 != 0) || year%400 == 0;
134 for (i=0; daytab[leap][i]<jday; i++) continue;
137 *dd = jday - daytab[leap][i];
140 /* Secretly defined in uf_to_radar.c */
141 Volume *copy_sweeps_into_volume(Volume *new_volume, Volume *old_volume);
143 /**********************************************************************/
145 /* RSL_dorade_to_radar */
147 /**********************************************************************/
148 Radar *RSL_dorade_to_radar(char *infile)
154 int iv, iray, iparam;
155 int range_bin1, gate_size;
157 int nbins, data_len, word_size;
160 float scale, offset, value;
161 int datum, missing_data_flag;
165 Range (*invf)(float x);
176 Platform_info *platform_info;
185 int year, month, day, jday, jday_vol;
188 if (infile == NULL) {
191 fp = fdopen(save_fd, "r");
193 if((fp=fopen(infile, "r"))==(FILE *)NULL) {
198 fp = uncompress_pipe(fp); /* Transparently, use gunzip. */
200 cb = dorade_read_comment_block(fp);
202 /**********************************************************************/
204 vd = dorade_read_volume_desc(fp); /* R E A D */
205 if (radar_verbose_flag) dorade_print_volume_desc(vd); /* P R I N T */
208 sd = (Sensor_desc **) calloc(vd->nsensors, sizeof(Sensor_desc *));
209 for (i=0; i<vd->nsensors; i++) {
210 sd[i] = dorade_read_sensor(fp);
214 if (radar_verbose_flag) {
215 for (i=0; i<vd->nsensors; i++) {
216 fprintf(stderr, "============ S E N S O R # %d =====================\n", i);
217 dorade_print_sensor(sd[i]);
220 /* R E A D sweeps. */
221 if (vd->nsensors > 1) {
222 fprintf(stderr, "RSL_dorade_to_radar: Unable to process for more than 1 sensor.\n");
223 fprintf(stderr, "RSL_dorade_to_radar: Number of sensors is %d\n", vd->nsensors);
227 /* Use sensor 0 for vitals. */
228 rd = sd[0]->radar_desc;
229 range_bin1 = sd[0]->cell_range_vector->range_cell[0];
230 gate_size = sd[0]->cell_range_vector->range_cell[1] - range_bin1;
232 radar = RSL_new_radar(MAX_RADAR_VOLUMES);
233 radar->h.month = vd->month;
234 radar->h.day = vd->day;
235 radar->h.year = vd->year;
236 radar->h.hour = vd->hour;
237 radar->h.minute = vd->minute;
238 radar->h.sec = vd->second;
239 sprintf(radar->h.radar_type, "dorade");
241 strncpy(radar->h.name, vd->flight_num, sizeof(radar->h.name));
242 strncpy(radar->h.radar_name, rd->radar_name, sizeof(radar->h.radar_name));
243 strncpy(radar->h.project, vd->project_name, sizeof(radar->h.project));
244 sprintf(radar->h.city, "Unknown");
245 strncpy(radar->h.state, "UKN", 3);
246 sprintf(radar->h.country, "Unknown");
247 /* Convert lat to d:m:s */
248 degree = (int)rd->latitude;
249 minute = (int)((rd->latitude - degree) * 60);
250 second = (rd->latitude - degree - minute/60.0) * 3600.0;
251 radar->h.latd = degree;
252 radar->h.latm = minute;
253 radar->h.lats = second;
254 /* Convert lat to d:m:s */
255 degree = (int)rd->longitude;
256 minute = (int)((rd->longitude - degree) * 60);
257 second = (rd->longitude - degree - minute/60.0) * 3600.0;
258 radar->h.lond = degree;
259 radar->h.lonm = minute;
260 radar->h.lons = second;
261 radar->h.height = rd->altitude * 1000.0;
262 radar->h.spulse = 0; /* FIXME */
263 radar->h.lpulse = 0; /* FIXME */
268 jday_vol = julian(year, month, day);
271 Get any threshold values present in parameter descriptors.
273 Note: The parameter descriptor contains an 8-character string which may
274 contain the, "Name of parameter upon which this parameter is thresholded".
275 According to documentation, if there is no thresholding parameter, the
276 string will be "NONE". In practice, this string has occasionally been
277 seen to contain only spaces for this condition. The corresponding
278 missing-data-flag for threshold may also vary, the nominal value being
279 -999, but has also been -32768.
282 nparam = rd->nparam_desc; [<--useable code, i.e., not pseudo]
283 create string array for threshold parameter names, size nparams.
285 strcpy thresh param name from param desc to threshold param array.
286 if string is all blanks, replace with 'NONE'.
290 /* Begin volume code. */
291 /* We don't know how many sweeps per volume exist, until we read
292 * the file. So allocate a large number of pointers, hope we don't
293 * exceed it, and adjust the pointer array at the end. This is
294 * efficient because we'll be manipulating pointers to the sweeps and
295 * not the sweeps themselves.
298 if (radar_verbose_flag)
299 fprintf(stderr, "Number of parameters: %d\n", rd->nparam_desc);
301 /* All the parameters are together, however, their order within
302 * the ray is not guarenteed. For instance, VE could appear after
303 * DM. For this we'll keep a list of parameter names and perform
304 * a (linear) search. The result will be an index into the RSL
305 * volume array (radar->v[i]). It is likely that the order will be
306 * consistant within a file, therefore, we'll keep track of the index of
307 * our previous parameter type and begin the search from there; the next
308 * index should be a match.
310 * The dorade parameter names and the rsl mapping is:
312 * Dorade: VE, DM, SW, DBZ, ZDR, PHI, RHO, LDR, DX, CH, AH, CV, AV
313 * RSL: VR, DM, SW, DZ, ZD, PH, RH, LR, *DX, *CH, *AH, *CV, *AV.
315 * * means this is a new RSL name.
318 #define DORADE_MAX_SWEEP 20
320 while((sr = dorade_read_sweep(fp, sd))) {
321 for(iray = 0; iray < sr->nrays; iray++) {
322 dray = sr->data_ray[iray];
324 /* Now, loop through the parameters and fill the rsl structures. */
325 for (iparam = 0; iparam < dray->nparam; iparam++) {
326 pd = dray->parameter_data[iparam];
327 iv = find_rsl_field_index(pd->name);
329 prt_skipped_field_msg(pd->name);
332 if (radar->v[iv] == NULL) {
333 radar->v[iv] = RSL_new_volume(DORADE_MAX_SWEEP); /* Expandable */
334 } else if (nsweep >= radar->v[iv]->h.nsweeps) {
335 /* Must expand the number of sweeps. */
336 /* Expand by another DORADE_MAX_SWEEP. */
337 if (radar_verbose_flag) {
338 fprintf(stderr, "nsweeps (%d) exceeds radar->v[%d]->h.nsweeps (%d)."
339 "\n", nsweep, iv, radar->v[iv]->h.nsweeps);
340 fprintf(stderr, "Increasing it to %d sweeps\n",
341 radar->v[iv]->h.nsweeps+DORADE_MAX_SWEEP);
343 new_volume = RSL_new_volume(radar->v[iv]->h.nsweeps+DORADE_MAX_SWEEP);
344 /* Look in uf_to_radar.c for 'copy_sweeps_into_volume' */
345 new_volume = copy_sweeps_into_volume(new_volume, radar->v[iv]);
346 radar->v[iv] = new_volume;
350 * Or just use RSL_f_list[iv] and RSL_invf_list[iv] as in sweep.h below.
353 /* Allocate the ray and load the parameter data. */
354 if ((sweep = radar->v[iv]->sweep[nsweep]) == NULL) {
355 sweep = radar->v[iv]->sweep[nsweep] = RSL_new_sweep(sr->s_info->nrays);
356 sweep->h.sweep_num = sr->s_info->sweep_num;
357 sweep->h.elev = sr->s_info->fixed_angle;
358 sweep->h.beam_width = rd->horizontal_beam_width;
359 sweep->h.vert_half_bw = rd->vertical_beam_width / 2.0;
361 sweep->h.vert_half_bw = radar->v[iv]->sweep[nsweep]->h.beam_width / 2.0;
363 sweep->h.horz_half_bw = rd->horizontal_beam_width / 2.0;
364 sweep->h.f = RSL_f_list[iv];
365 sweep->h.invf = RSL_invf_list[iv];
369 data_len = dray->data_len[iparam];
370 word_size = dray->word_size[iparam];
371 if (word_size != 2 && word_size != 4) {
372 fprintf(stderr,"RSL_dorade_to_radar: dray->word_size[%d] = %d\n",
373 iparam, dray->word_size[iparam]);
374 fprintf(stderr,"Expected value is 2 or 4.\n");
377 nbins = data_len / word_size;
379 if ((ray = sweep->ray[iray]) == NULL) {
380 if (radar_verbose_flag)
381 fprintf(stderr, "Allocating %d bins for ray %d\n",
382 dray->data_len[iparam], iray);
383 ray = sweep->ray[iray] = RSL_new_ray(nbins);
386 /* Load ray header. */
388 ray_info = dray->ray_info;
389 jday = ray_info->jday;
390 if (jday != jday_vol) {
391 if (jday > 0 && jday < 367) {
392 /* Recompute year month day */
393 if (jday < jday_vol) year++;
394 ymd(jday, year, &month, &day);
397 else { /* Invalid jday */
401 ray->h.month = month;
403 ray->h.hour = ray_info->hour;
404 ray->h.minute = ray_info->minute;
405 ray->h.sec = ray_info->second + ray_info->msec / 1000.;
406 ray->h.azimuth = ray_info->azimuth;
407 ray->h.ray_num = iray + 1;
408 ray->h.elev = ray_info->elevation;
409 ray->h.elev_num = ray_info->sweep_num;
410 ray->h.unam_rng = rd->unambiguous_range;
411 ray->h.nyq_vel = rd->unambiguous_velocity;
412 ray->h.azim_rate = ray_info->scan_rate;
416 ray->h.range_bin1 = range_bin1;
417 ray->h.gate_size = gate_size;
418 platform_info = dray->platform_info;
419 ray->h.pitch = platform_info->pitch;
420 ray->h.roll = platform_info->roll;
421 ray->h.heading = platform_info->heading;
422 ray->h.pitch_rate = platform_info->pitch_rate;
423 ray->h.heading_rate = platform_info->heading_rate;
424 ray->h.lat = platform_info->latitude;
425 ray->h.lon = platform_info->longitude;
426 ray->h.alt = platform_info->altitude;
427 ray->h.vel_east = platform_info->ew_speed;
428 ray->h.vel_north = platform_info->ns_speed;
429 ray->h.vel_up = platform_info->v_speed;
430 /* TODO: Does DORADE have Doppler velocity res?
434 ray->h.beam_width = rd->horizontal_beam_width;
436 ray->h.frequency = @Radar Descriptor
437 Get parm_desc.xmit_freq
438 Find first set-bit, use frequency from corresponding index in rad_desc frequenies.
440 ray->h.pulse_count = Is it number samples used? @Param desc. Probably not.
442 * The following are DONE *
443 ray->h.pulse_width = @Parameter Descriptor--not same--it's in meters--
444 we use micro-sec. They're using pulse length.
445 pulse width = pulse length/c
446 ray->h.wavelength = Can compute using Nyquist velocity and PRF or PRT:
447 wl = 4*vmax/PRF or wl = 4*vmax*PRT
448 ray->h.prf = Can compute if nothing else: prf = c/(2*Rmax)
450 /* DORADE is actually using pulse length. Convert to pulse width. */
451 ray->h.pulse_width = (sd[0]->p_desc[iparam]->pulse_width/
452 RSL_SPEED_OF_LIGHT) * 1e6;
453 prf = RSL_SPEED_OF_LIGHT/(2.*rd->unambiguous_range*1000.);
455 ray->h.wavelength = (4.*rd->unambiguous_velocity)/prf;
456 ray->h.nbins = nbins;
457 ray->h.f = RSL_f_list[iv];
458 ray->h.invf = RSL_invf_list[iv];
460 /* Copy the ray data into the RSL ray. */
462 /* .... fill here .... */
464 /* Assign pointer to data.
465 * Get datum using word size and proper cast.
466 * Convert and store in rsl ray->range.
467 * Increment pointer to data based on word size.
470 ptr2bytes = (short *) pd->data;
471 ptr4bytes = (int *) pd->data;
472 scale = sd[0]->p_desc[iparam]->scale_factor;
473 offset = sd[0]->p_desc[iparam]->offset_factor;
474 missing_data_flag = sd[0]->p_desc[iparam]->missing_data_flag;
476 for (i=0; i<nbins; i++) {
477 if (word_size == 2) datum = *ptr2bytes++;
478 else datum = *ptr4bytes++;
480 TODO: If there is a threshold parameter for this parameter then
481 apply threshold value. I think threshold works like this: If there's a
482 threshold parameter, as with VT (threshold parm = NCP), then use
483 threshold value from that parameter unless it is the missing data value.
485 if (datum != missing_data_flag) {
486 value = ((float)datum - offset) / scale;
490 ray->range[i] = ray->h.invf(value);
494 radar->v[iv]->h.nsweeps = nsweep + 1;
495 radar->v[iv]->h.f = RSL_f_list[iv];
496 radar->v[iv]->h.invf = RSL_invf_list[iv];
501 if (radar_verbose_flag) fprintf(stderr, "______NEW SWEEP__<%d>____\n", nsweep);
502 /* Save for loading into volume structure. */
503 dorade_free_sweep(sr);
506 /* The following avoids a broken pipe message, since a VOLD at the end
509 while(fread(buf, sizeof(buf), 1, fp)) continue; /* Read til EOF */