/* NASA/TRMM, Code 910.1. This is the TRMM Office Radar Software Library. Copyright (C) 1996-1999 John H. Merritt Space Applications Corporation Vienna, Virginia This library is free software; you can redistribute it and/or modify it under the terms of the GNU Library General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public License for more details. You should have received a copy of the GNU Library General Public License along with this library; if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */ #include #include #include #include #define USE_RSL_VARS #include "rsl.h" #include "dorade.h" extern int radar_verbose_flag; /********************************************************************/ /* */ /* find_rsl_field_index */ /* */ /********************************************************************/ int find_rsl_field_index(char *dorade_field_name) { /* * Dorade: VE, DM, SW, DBZ, ZDR, PHI, RHO, LDR, DX, CH, AH, CV, AV * RSL: VR, DM, SW, DZ, ZD, PH, RH, LR, *DX, *CH, *AH, *CV, *AV. */ if (strncasecmp(dorade_field_name, "ve", 2) == 0) return VR_INDEX; if (strncasecmp(dorade_field_name, "vr", 2) == 0) return VR_INDEX; if (strncasecmp(dorade_field_name, "dm", 2) == 0) return DM_INDEX; if (strncasecmp(dorade_field_name, "sw", 2) == 0) return SW_INDEX; if (strncasecmp(dorade_field_name, "dbz", 3) == 0) return DZ_INDEX; if (strncasecmp(dorade_field_name, "zdr", 3) == 0) return ZD_INDEX; if (strncasecmp(dorade_field_name, "phi", 3) == 0) return PH_INDEX; if (strncasecmp(dorade_field_name, "rho", 3) == 0) return RH_INDEX; if (strncasecmp(dorade_field_name, "ldr", 3) == 0) return LR_INDEX; if (strncasecmp(dorade_field_name, "dx", 2) == 0) return DX_INDEX; if (strncasecmp(dorade_field_name, "ch", 2) == 0) return CH_INDEX; if (strncasecmp(dorade_field_name, "ah", 2) == 0) return AH_INDEX; if (strncasecmp(dorade_field_name, "cv", 2) == 0) return CV_INDEX; if (strncasecmp(dorade_field_name, "av", 2) == 0) return AV_INDEX; if (strncasecmp(dorade_field_name, "vs", 2) == 0) return VS_INDEX; if (strncasecmp(dorade_field_name, "vl", 2) == 0) return VL_INDEX; if (strncasecmp(dorade_field_name, "vg", 2) == 0) return VG_INDEX; if (strncasecmp(dorade_field_name, "vt", 2) == 0) return VT_INDEX; if (strncasecmp(dorade_field_name, "ncp", 2) == 0) return NP_INDEX; return -1; } void prt_skipped_field_msg(char *dorade_field_name) { char prtname[9]; int i, already_printed; #define MAXFIELDS 20 static int nskipped = 0; static char skipped_list[MAXFIELDS][9]; /* Make sure name is a properly formed string. */ strncpy(prtname, dorade_field_name, 8); prtname[8] = '\0'; /* We don't want to repeat messages for the same field, so keep a list of * fields already printed. */ already_printed = 0; i = 0; while (!already_printed && i < nskipped) { if (strncmp(prtname, skipped_list[i], 2) == 0) already_printed = 1; i++; } if (!already_printed) { fprintf(stderr, "Unknown DORADE parameter type <%s> -- skipping.\n", prtname); strcpy(skipped_list[nskipped], prtname); nskipped++; } } /* For date conversion routines. */ static int daytab[2][13] = { {0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334, 365}, {0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335, 366} }; /*************************************************************/ /* */ /* julian */ /* */ /*************************************************************/ static int julian(int year, int mo, int day) { /* Converts a calendar date (month, day, year) to a Julian date. Returns: Julian day. */ int leap; leap = (year%4 == 0 && year%100 != 0) || year%400 == 0; return(day + daytab[leap][mo-1]); } /*************************************************************/ /* */ /* ymd */ /* */ /*************************************************************/ static void ymd(int jday, int year, int *mm, int *dd) { /* Input: jday, yyyy */ /* Output: mm, dd */ /* Copied from hdf_to_radar.c, written by Mike Kolander. */ int leap; int i; leap = (year%4 == 0 && year%100 != 0) || year%400 == 0; for (i=0; daytab[leap][i]nsensors, sizeof(Sensor_desc *)); for (i=0; insensors; i++) { sd[i] = dorade_read_sensor(fp); } /* P R I N T */ if (radar_verbose_flag) { for (i=0; insensors; i++) { fprintf(stderr, "============ S E N S O R # %d =====================\n", i); dorade_print_sensor(sd[i]); } } /* R E A D sweeps. */ if (vd->nsensors > 1) { fprintf(stderr, "RSL_dorade_to_radar: Unable to process for more than 1 sensor.\n"); fprintf(stderr, "RSL_dorade_to_radar: Number of sensors is %d\n", vd->nsensors); return NULL; } /* Use sensor 0 for vitals. */ rd = sd[0]->radar_desc; range_bin1 = sd[0]->cell_range_vector->range_cell[0]; gate_size = sd[0]->cell_range_vector->range_cell[1] - range_bin1; radar = RSL_new_radar(MAX_RADAR_VOLUMES); radar->h.month = vd->month; radar->h.day = vd->day; radar->h.year = vd->year; radar->h.hour = vd->hour; radar->h.minute = vd->minute; radar->h.sec = vd->second; sprintf(radar->h.radar_type, "dorade"); radar->h.number = 0; strncpy(radar->h.name, vd->flight_num, sizeof(radar->h.name)); strncpy(radar->h.radar_name, rd->radar_name, sizeof(radar->h.radar_name)); strncpy(radar->h.project, vd->project_name, sizeof(radar->h.project)); sprintf(radar->h.city, "Unknown"); strncpy(radar->h.state, "UKN", 3); sprintf(radar->h.country, "Unknown"); /* Convert lat to d:m:s */ degree = (int)rd->latitude; minute = (int)((rd->latitude - degree) * 60); second = (rd->latitude - degree - minute/60.0) * 3600.0; radar->h.latd = degree; radar->h.latm = minute; radar->h.lats = second; /* Convert lat to d:m:s */ degree = (int)rd->longitude; minute = (int)((rd->longitude - degree) * 60); second = (rd->longitude - degree - minute/60.0) * 3600.0; radar->h.lond = degree; radar->h.lonm = minute; radar->h.lons = second; radar->h.height = rd->altitude * 1000.0; radar->h.spulse = 0; /* FIXME */ radar->h.lpulse = 0; /* FIXME */ year = vd->year; month = vd->month; day = vd->day; jday_vol = julian(year, month, day); /* TODO: Get any threshold values present in parameter descriptors. Note: The parameter descriptor contains an 8-character string which may contain the, "Name of parameter upon which this parameter is thresholded". According to documentation, if there is no thresholding parameter, the string will be "NONE". In practice, this string has occasionally been seen to contain only spaces for this condition. The corresponding missing-data-flag for threshold may also vary, the nominal value being -999, but has also been -32768. Pseudo code: nparam = rd->nparam_desc; [<--useable code, i.e., not pseudo] create string array for threshold parameter names, size nparams. for each param: strcpy thresh param name from param desc to threshold param array. if string is all blanks, replace with 'NONE'. endfor */ /* Begin volume code. */ /* We don't know how many sweeps per volume exist, until we read * the file. So allocate a large number of pointers, hope we don't * exceed it, and adjust the pointer array at the end. This is * efficient because we'll be manipulating pointers to the sweeps and * not the sweeps themselves. */ if (radar_verbose_flag) fprintf(stderr, "Number of parameters: %d\n", rd->nparam_desc); /* All the parameters are together, however, their order within * the ray is not guarenteed. For instance, VE could appear after * DM. For this we'll keep a list of parameter names and perform * a (linear) search. The result will be an index into the RSL * volume array (radar->v[i]). It is likely that the order will be * consistant within a file, therefore, we'll keep track of the index of * our previous parameter type and begin the search from there; the next * index should be a match. * * The dorade parameter names and the rsl mapping is: * * Dorade: VE, DM, SW, DBZ, ZDR, PHI, RHO, LDR, DX, CH, AH, CV, AV * RSL: VR, DM, SW, DZ, ZD, PH, RH, LR, *DX, *CH, *AH, *CV, *AV. * * * means this is a new RSL name. */ #define DORADE_MAX_SWEEP 20 nsweep = 0; while((sr = dorade_read_sweep(fp, sd))) { for(iray = 0; iray < sr->nrays; iray++) { dray = sr->data_ray[iray]; /* Now, loop through the parameters and fill the rsl structures. */ for (iparam = 0; iparam < dray->nparam; iparam++) { pd = dray->parameter_data[iparam]; iv = find_rsl_field_index(pd->name); if (iv < 0) { prt_skipped_field_msg(pd->name); continue; } if (radar->v[iv] == NULL) { radar->v[iv] = RSL_new_volume(DORADE_MAX_SWEEP); /* Expandable */ } else if (nsweep >= radar->v[iv]->h.nsweeps) { /* Must expand the number of sweeps. */ /* Expand by another DORADE_MAX_SWEEP. */ if (radar_verbose_flag) { fprintf(stderr, "nsweeps (%d) exceeds radar->v[%d]->h.nsweeps (%d)." "\n", nsweep, iv, radar->v[iv]->h.nsweeps); fprintf(stderr, "Increasing it to %d sweeps\n", radar->v[iv]->h.nsweeps+DORADE_MAX_SWEEP); } new_volume = RSL_new_volume(radar->v[iv]->h.nsweeps+DORADE_MAX_SWEEP); /* Look in uf_to_radar.c for 'copy_sweeps_into_volume' */ new_volume = copy_sweeps_into_volume(new_volume, radar->v[iv]); radar->v[iv] = new_volume; } /* Assign f and invf switch (iv) . . . * Or just use RSL_f_list[iv] and RSL_invf_list[iv] as in sweep.h below. */ /* Allocate the ray and load the parameter data. */ if ((sweep = radar->v[iv]->sweep[nsweep]) == NULL) { sweep = radar->v[iv]->sweep[nsweep] = RSL_new_sweep(sr->s_info->nrays); sweep->h.sweep_num = sr->s_info->sweep_num; sweep->h.elev = sr->s_info->fixed_angle; sweep->h.beam_width = rd->horizontal_beam_width; sweep->h.vert_half_bw = rd->vertical_beam_width / 2.0; /* sweep->h.vert_half_bw = radar->v[iv]->sweep[nsweep]->h.beam_width / 2.0; */ sweep->h.horz_half_bw = rd->horizontal_beam_width / 2.0; sweep->h.f = RSL_f_list[iv]; sweep->h.invf = RSL_invf_list[iv]; } data_len = dray->data_len[iparam]; word_size = dray->word_size[iparam]; if (word_size != 2 && word_size != 4) { fprintf(stderr,"RSL_dorade_to_radar: dray->word_size[%d] = %d\n", iparam, dray->word_size[iparam]); fprintf(stderr,"Expected value is 2 or 4.\n"); return NULL; } nbins = data_len / word_size; if ((ray = sweep->ray[iray]) == NULL) { if (radar_verbose_flag) fprintf(stderr, "Allocating %d bins for ray %d\n", dray->data_len[iparam], iray); ray = sweep->ray[iray] = RSL_new_ray(nbins); } /* Load ray header. */ ray_info = dray->ray_info; jday = ray_info->jday; if (jday != jday_vol) { if (jday > 0 && jday < 367) { /* Recompute year month day */ if (jday < jday_vol) year++; ymd(jday, year, &month, &day); jday_vol = jday; } else { /* Invalid jday */ } } ray->h.year = year; ray->h.month = month; ray->h.day = day; ray->h.hour = ray_info->hour; ray->h.minute = ray_info->minute; ray->h.sec = ray_info->second + ray_info->msec / 1000.; ray->h.azimuth = ray_info->azimuth; ray->h.ray_num = iray + 1; ray->h.elev = ray_info->elevation; ray->h.elev_num = ray_info->sweep_num; ray->h.unam_rng = rd->unambiguous_range; ray->h.nyq_vel = rd->unambiguous_velocity; ray->h.azim_rate = ray_info->scan_rate; /* TODO ray->h.fix_angle = ; */ ray->h.range_bin1 = range_bin1; ray->h.gate_size = gate_size; platform_info = dray->platform_info; ray->h.pitch = platform_info->pitch; ray->h.roll = platform_info->roll; ray->h.heading = platform_info->heading; ray->h.pitch_rate = platform_info->pitch_rate; ray->h.heading_rate = platform_info->heading_rate; ray->h.lat = platform_info->latitude; ray->h.lon = platform_info->longitude; ray->h.alt = platform_info->altitude; ray->h.vel_east = platform_info->ew_speed; ray->h.vel_north = platform_info->ns_speed; ray->h.vel_up = platform_info->v_speed; /* TODO: Does DORADE have Doppler velocity res? A: No. ray->h.vel_res = N/A */ ray->h.beam_width = rd->horizontal_beam_width; /* TODO ray->h.frequency = @Radar Descriptor Get parm_desc.xmit_freq Find first set-bit, use frequency from corresponding index in rad_desc frequenies. ray->h.pulse_count = Is it number samples used? @Param desc. Probably not. * The following are DONE * ray->h.pulse_width = @Parameter Descriptor--not same--it's in meters-- we use micro-sec. They're using pulse length. pulse width = pulse length/c ray->h.wavelength = Can compute using Nyquist velocity and PRF or PRT: wl = 4*vmax/PRF or wl = 4*vmax*PRT ray->h.prf = Can compute if nothing else: prf = c/(2*Rmax) */ /* DORADE is actually using pulse length. Convert to pulse width. */ ray->h.pulse_width = (sd[0]->p_desc[iparam]->pulse_width/ RSL_SPEED_OF_LIGHT) * 1e6; prf = RSL_SPEED_OF_LIGHT/(2.*rd->unambiguous_range*1000.); ray->h.prf = prf; ray->h.wavelength = (4.*rd->unambiguous_velocity)/prf; ray->h.nbins = nbins; ray->h.f = RSL_f_list[iv]; ray->h.invf = RSL_invf_list[iv]; /* Copy the ray data into the RSL ray. */ /* .... fill here .... */ /* Assign pointer to data. * Get datum using word size and proper cast. * Convert and store in rsl ray->range. * Increment pointer to data based on word size. */ ptr2bytes = (short *) pd->data; ptr4bytes = (int *) pd->data; scale = sd[0]->p_desc[iparam]->scale_factor; offset = sd[0]->p_desc[iparam]->offset_factor; missing_data_flag = sd[0]->p_desc[iparam]->missing_data_flag; for (i=0; irange[i] = ray->h.invf(value); } if (iray == 0) { radar->v[iv]->h.nsweeps = nsweep + 1; radar->v[iv]->h.f = RSL_f_list[iv]; radar->v[iv]->h.invf = RSL_invf_list[iv]; } } } nsweep++; if (radar_verbose_flag) fprintf(stderr, "______NEW SWEEP__<%d>____\n", nsweep); /* Save for loading into volume structure. */ dorade_free_sweep(sr); } /* The following avoids a broken pipe message, since a VOLD at the end * is not read yet. */ while(fread(buf, sizeof(buf), 1, fp)) continue; /* Read til EOF */ rsl_pclose(fp); return radar; }