X-Git-Url: http://pileus.org/git/?p=~andy%2Frsl;a=blobdiff_plain;f=dorade_to_radar.c;h=b66b9a8beb85adf627ec5383e71bd03116f0b0d0;hp=ff963a8a85f335b66ace0191abe821cc2a7f238f;hb=d08a7f8a699a044bc4ac5f93917aa7f6c463923b;hpb=012916676d26251849e408aaf574458e196df2c4 diff --git a/dorade_to_radar.c b/dorade_to_radar.c index ff963a8..b66b9a8 100644 --- a/dorade_to_radar.c +++ b/dorade_to_radar.c @@ -25,6 +25,7 @@ #include #include #include +#include #define USE_RSL_VARS #include "rsl.h" #include "dorade.h" @@ -39,30 +40,112 @@ extern int radar_verbose_flag; int find_rsl_field_index(char *dorade_field_name) { /* - * Dorade: VE, DM, SW, DZ, ZDR, PHI, RHO, LDR, DX, CH, AH, CV, AV - * RSL: VR, DM, SW, DZ, ZD, PH, RH, LR, *DX, *CH, *AH, *CV, *AV. + * 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, "dm", 2) == 0) return DM_INDEX; - if (strncasecmp(dorade_field_name, "sw", 2) == 0) return SW_INDEX; - if (strncasecmp(dorade_field_name, "dz", 2) == 0) return DZ_INDEX; + 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, "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; - fprintf(stderr, "Unknown DORADE type <%s>\n", dorade_field_name); 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); + 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]); - } + 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; + 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; @@ -159,7 +262,32 @@ Radar *RSL_dorade_to_radar(char *infile) 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 @@ -169,21 +297,21 @@ Radar *RSL_dorade_to_radar(char *infile) */ if (radar_verbose_flag) - fprintf(stderr, "Number of parameters: %d\n", rd->nparam_desc); + 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 + * 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, DZ, ZDR, PHI, RHO, LDR, DX, CH, AH, CV, AV - * RSL: VR, DM, SW, DZ, ZD, PH, RH, LR, *DX, *CH, *AH, *CV, *AV. + * 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. */ @@ -191,56 +319,189 @@ Radar *RSL_dorade_to_radar(char *infile) #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 (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, "Increading 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; - } - - /* 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; + 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 = 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]; - } - - - 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(dray->data_len[iparam]); - } - - /* Copy the ray data into the RSL ray. */ + 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); } - nsweep++; - if (radar_verbose_flag) fprintf(stderr, "______NEW SWEEP__<%d>____\n", nsweep); - /* Save for loading into volume structure. */ - dorade_free_sweep(sr); + + 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 @@ -252,4 +513,3 @@ Radar *RSL_dorade_to_radar(char *infile) return radar; } -