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