#!/usr/bin/perl
#======================================================================
# L I S T B I N S
# doc: Fri Aug 25 15:57:05 2006
# dlm: Fri Jun 5 13:45:23 2020
# (c) 2006 A.M. Thurnherr
# uE-Info: 367 103 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# Split data file into per-bin time series.
# HISTORY:
# Aug 25, 2006: - created from [listEns]
# Aug 26, 2006: - added -M)agdecl
# - changed -b to -f
# Aug 27, 2006: - added %bin
# Aug 28, 2006: - BUG: options were confused
# Jan 4, 2007: - improved usage message for -a
# - added %mag_decl
# - BUG: roundoff error in %pct_good_vels
# Sep 19, 2007: - adapted to new [RDI_BB_Read.pl]
# Feb 7, 2008: - added sound-speed correction
# - enabled 3-beam solutions
# Feb 8, 2008: - added -d)iscard <beam>
# - added -b)eam coordinate output
# Feb 12, 2008: - modified 3-beam output
# - added -p)ct_good <min>
# Feb 13, 2008: - various improvements
# Feb 19, 2008: - BUG: division by zero
# BUG: min() did not work with 1st elt undef
# Feb 21, 2008: - BUG: had forgotten to undo debugging changes
# - removed missing magdecl warning on -b
# Feb 22, 2008: - moved ssCorr() to [RDI_Utils.pl]
# - BUG: %d complete ensembles was written to STDOUT
# - BUG: %N_ensembles was total number of ensembles, not
# only the ones used
# - BUG: 0 errvel was output for 3-beam solutions =>
# wrong time-average statistics
# May 19, 2009: - added -w to calculate vertical velocities from
# two beam pairs separately
# May 21, 2009: - added horizontal beampair velocities on -w
# - -P)itchRoll <bias/bias>
# May 22, 2009: - added -B) <bias/bias/bias/bias>
# May 23, 2009: - adapted to changed beampair-velocity fun name
# Aug 22, 2010: - added -R
# Apr 29, 2013: - cosmetics
# - added warning on missing -S
# Aug 7, 2013: - BUG: -w did not respect -d
# Mar 4, 2014: - added support for missing PITCH/ROLL/HEADING
# Apr 22, 2014: - BUG: extraneous semicolon
# Nov 24, 2014: - enabled -w always
# Mar 22, 2015: - replaced -f by -o (allowing for pipes)
# Jun 16, 2015: - BUG: velocity bias code did not respect bad velocities
# Jan 5, 2016: - adapted to [ANTS_tools_lib.pl]
# - adapted to calculation of w12, w34 from earth-coordinate data
# - several other changes to the code that should not affect the results
# Jan 6, 2016: - -b removed (always output beamvels)
# Jan 31, 2016: - started debugging the obviously wrong Earth2Beam() transformation
# Feb 29, 2016: - continued debugging; removed debugging code
# Mar 17, 2016: - adapted to new Getopt library
# Jun 9, 2016: - minor improvements
# - BUG: velBeamToEarth() has new interface
# Jan 27, 2017: - BUG: earth-coordinate label was wrong
# - changed semantics to include all %-good fields of earth-coordinate data
# and also report total pcg (3+4 beam) in min_pcg
# - BUG: 3-beam %ages were incorrect: 1) they were based on goodvels instead of the
# entire ensemble range; 2) pings-per-ensemble were not considered
# - BUG: output layout was all messed up for non-valid velocities
# Feb 6, 2018: - added support for PD0_IO first_ens, last_ens
# Apr 10, 2018: - added day number to output
# - added -l)ast bin
# - activate output files
# Aug 29, 2018: - added error message on -r decoding failures
# Jun 13, 2018: - adpated to RTI files (disabled BIT error check)
# - BUG: dn did not have sufficient digits
# Feb 13, 2020: - added -z
# May 11, 2020: - removed -z, added -t -m
# General Notes:
# - everything (e.g. beams) is numbered from 1
# - no support for BT data
# Post-Processing:
# - e.g. in a Makefile,
# BFILT = | list -Lconv ensemble elapsed dn12='str2dec_time(\$$date,\$$time,2012)' heading pitch roll temp u v w err_vel w12 w34 corr1 corr2 corr3 corr4 amp1 amp2 amp3 amp4 > `basename $@ .b1`.b%d
# - used as follows
# listBins -o "${BFILT}"
# - this writes only the interesting fields and calculates a 2012-referenced decimal day field
# Soundspeed Correction:
# - applied as described in the RDI coord-trans manual
# - sound-speed variation over range is ignored (valid for small gradients)
# => - same simple correction for all velocity components
# - simple correction for cell depths
# Min %-good (min_pcg):
# - nan for records w/o valid velocities
# - min(%-good) of the beams used for the velocity solution
# - min_pcg does not have to decrease monotonically with distance,
# at least when 3-beam solutions are allowed and when -p is used to
# edit the data
# - non-monotonic min_pcg is particularly obvious with the DYNAMUCK BM_ADCP
# data, where one of the beams performed much worse than the others
# %-good of earth-coordinate data:
# pcg1 3-beam (correlation threshold)
# pcg2 bad (errvel threshold)
# pcg3 bad (>1 bad beam)
# pcg4 4-beam
use Getopt::Std;
$ADCP_tools_minVersion = 2.2;
($ADCP_TOOLS) = ($0 =~ m{(.*/)[^/]+});
require "$ADCP_TOOLS/ADCP_tools_lib.pl";
$antsMinLibVersion = 7.0;
($ANTS) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
require "$ANTS/ants.pl";
require "$ANTS/libconv.pl";
die("Usage: $0 [-r)ange <first_ens,last_ens>] [-l)ast <bin>] [-R)enumber ensembles from 1] " .
"[-o)utput <redirection[>bin%d.raw]>] " .
"[output -a)ll ens (not just those with good vels)] " .
"[-M)agnetic <declination>] " .
"[-S)oundspeed correction <salin|*,temp|*,depth|*> " .
"[Instrument -T)ransformation Matrix <file>] " .
'[disable bin -m)apping] [use TRDI beam-to-earth -t)ransformation] ' .
"[-P)itch/Roll <bias/bias>] [-B)eamvel <bias/bias/bias/bias>] " .
"[require -4)-beam solutions] [-d)iscard <beam#>] " .
"[-p)ct-good <min>] " .
"<RDI file>\n")
unless (&getopts("4aB:d:l:mM:o:p:r:P:RS:tT:") && @ARGV == 1);
($P{pitch_bias},$P{roll_bias}) = split('[,/]',$opt_P);
($P{velbias_b1},$P{velbias_b2},$P{velbias_b3},$P{velbias_b4}) = split('[,/]',$opt_B);
die("$0: -4 and -d are mutually exclusive\n")
if ($opt_4 && defined($opt_d));
$opt_p = 0 unless defined($opt_p);
$RDI_Coords::minValidVels = 4 if ($opt_4); # no 3-beam solutions
$RDI_Coords::binMapping = 'none' if ($opt_m); # 'linterp' is default
$RDI_Coords::beamTransformation = 'RDI' if ($opt_t); # 'LHR90' is default
print(STDERR "WARNING: magnetic declination not set!\n")
unless defined($opt_M);
$opt_o = '>bin%d.raw' unless defined($opt_o);
$ifn = $ARGV[0];
if (defined($opt_r)) {
($first_ens,$last_ens) = split(',',$opt_r);
die("Cannot decode -r $opt_r\n")
unless ($first_ens>0 && $last_ens>$first_ens);
}
if (defined($opt_S)) {
($SS_salin,$SS_temp,$SS_depth) = split(',',$opt_S);
$variable_ssCorr = ($SS_salin eq '*' || $SS_temp eq '*' || $SS_depth eq '*');
} else {
print(STDERR "WARNING: no soundspeed correction applied!\n");
}
loadInstrumentTransformation($opt_T) # load instrument-transformation matrix
if (defined($opt_T));
#----------------------------------------------------------------------
sub min(@) # return minimum
{
my($min) = 99e99;
for (my($i)=0; $i<=$#_; $i++) {
$min = $_[$i] if defined($_[$i]) && ($_[$i] < $min);
}
return ($min == 99e99) ? nan : $min;
}
sub dumpBin($$$) # write time series of single bin
{
my($b,$fe,$le) = @_;
my($out) = sprintf($opt_o,$b+1);
open(P,"$out") || die("$out: $!\n");
print(P "#!/usr/bin/perl -S list\n");
chmod(0777&~umask,*P); # activate output
print(P "#ANTS#PARAMS# ");
foreach my $k (keys(%P)) {
print(P "$k\{$P{$k}\} ");
}
my($pct3b) = 100*$three_beam[$b]/(($le-$fe+1) * $dta{PINGS_PER_ENSEMBLE});
# die("assertion failed ($pct3b > 100): gv[$b] = $good_vels[$b]; 3b = $three_beam[$b]; ens=$fe..$le; ")
# if ($pct3b > 100);
printf(STDERR "%02d:%.0f%%/%.0f%% ",$b+1,100*$good_vels[$b]/($le-$fe+1),$pct3b);
printf(P "pct_3_beam{%.0f} ",$pct3b);
printf(P "pct_good_vels{%.0f} ",100*$good_vels[$b]/($le-$fe+1));
printf(P "bin{%d}",$b+1);
printf(P " soundspeed_correction{%s}",defined($opt_S) ? $opt_S : 'NONE!');
printf(P " dz{%g}",$dz[$b] *
(defined($opt_S) ? ssCorr($dta{ENSEMBLE}[$fe],$SS_salin,$SS_temp,$SS_depth) : 1)
) unless ($variable_ssCorr);
print( P "\n");
print(P "#ANTS#FIELDS# " .
"{ensemble} {date} {time} {elapsed} {dn} " .
"{heading} {pitch} {roll} " .
"{sig_heading} {sig_pitch} {sig_roll} " .
"{xmit_current} {xmit_voltage} " .
"{temp} " .
"{bv1} {bv2} {bv3} {bv4} {u} {v} {w} {err_vel} " .
"{v12} {w12} {v34} {w34} " .
"{corr1} {corr2} {corr3} {corr4} " .
"{amp1} {amp2} {amp3} {amp4} " .
"{pcg1} {pcg2} {pcg3} {pcg4} {3_beam} {min_pcg}"
);
print(P " {dz}") if ($variable_ssCorr);
print(P "\n");
my($t0) = $dta{ENSEMBLE}[$fe]->{UNIX_TIME};
for (my($e)=$fe; $e<=$le; $e++) {
next unless ($opt_a || $dta{ENSEMBLE}[$e]->{GOOD_VEL}[$b]);
my($ssCorr) = defined($opt_S) ? ssCorr($dta{ENSEMBLE}[$e],$SS_salin,$SS_temp,$SS_depth) : 1;
print(P "$dta{ENSEMBLE}[$e]->{NUMBER} ");
print(P "$dta{ENSEMBLE}[$e]->{DATE} ");
print(P "$dta{ENSEMBLE}[$e]->{TIME} ");
printf(P "%d ",$dta{ENSEMBLE}[$e]->{UNIX_TIME}-$t0); # elapsed time
printf(P "%.5f ",str2dec_time($dta{ENSEMBLE}[$e]->{DATE},$dta{ENSEMBLE}[$e]->{TIME})); # decimal day
print(P defined($dta{ENSEMBLE}[$e]->{HEADING}) ? "$dta{ENSEMBLE}[$e]->{HEADING} " : 'nan ');
print(P defined($dta{ENSEMBLE}[$e]->{PITCH}) ? "$dta{ENSEMBLE}[$e]->{PITCH} " : 'nan ');
print(P defined($dta{ENSEMBLE}[$e]->{ROLL}) ? "$dta{ENSEMBLE}[$e]->{ROLL} " : 'nan ');
print(P "$dta{ENSEMBLE}[$e]->{HEADING_STDDEV} ");
print(P "$dta{ENSEMBLE}[$e]->{PITCH_STDDEV} ");
print(P "$dta{ENSEMBLE}[$e]->{ROLL_STDDEV} ");
print(P "$dta{ENSEMBLE}[$e]->{ADC_XMIT_CURRENT} ");
print(P "$dta{ENSEMBLE}[$e]->{ADC_XMIT_VOLTAGE} ");
print(P "$dta{ENSEMBLE}[$e]->{TEMPERATURE} ");
if ($dta{ENSEMBLE}[$e]->{GOOD_VEL}[$b]) {
printf(P "%g %g %g %g ",@{$dta{ENSEMBLE}[$e]->{BEAM_VELOCITY}[$b]});
printf(P "%g ",$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0] * $ssCorr);
printf(P "%g ",$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][1] * $ssCorr);
printf(P "%g ",$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][2] * $ssCorr);
if ($dta{ENSEMBLE}[$e]->{THREE_BEAM}[$b]) {
print(P "nan ");
} else {
printf(P "%g ",$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3] * $ssCorr);
}
printf(P defined($dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][0]) ? "%g " : "nan ",
$dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][0]);
printf(P defined($dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][1]) ? "%g " : "nan ",
$dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][1]);
printf(P defined($dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][2]) ? "%g " : "nan ",
$dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][2]);
printf(P defined($dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][3]) ? "%g " : "nan ",
$dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][3]);
} else {
print(P "nan nan nan nan ");
print(P "nan nan nan nan ");
print(P "nan nan nan nan ");
}
print(P "@{$dta{ENSEMBLE}[$e]->{CORRELATION}[$b]} ");
print(P "@{$dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b]} ");
if ($dta{ENSEMBLE}[$e]->{GOOD_VEL}[$b]) {
print(P "@{$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b]} ");
} else {
print(P "nan nan nan nan ");
}
if ($beamCoords) {
printf(P "%d ",$dta{ENSEMBLE}[$e]->{THREE_BEAM}[$b]);
printf(P "%s ",min(@{$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b]}));
} else {
if ($dta{ENSEMBLE}[$e]->{GOOD_VEL}[$b]) {
print(P "$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][0] ");
} else {
print(P "nan ");
}
printf(P "%g ",$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][0]+$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][3]);
}
printf(P "%g ",$dz[$b]*$ssCorr) if ($variable_ssCorr);
print(P "\n");
}
close(P);
}
#----------------------------------------------------------------------
# MAIN
#----------------------------------------------------------------------
$P{RDI_file} = $ifn;
$P{mag_decl} = $opt_M if defined($opt_M);
readData($ifn,\%dta,$first_ens,$last_ens,$opt_l);
printf(STDERR "%d complete ensembles...\n",scalar(@{$dta{ENSEMBLE}}));
$dta{HEADING_BIAS} = -$opt_M; # magnetic declination
if ($dta{BEAM_COORDINATES}) { # coords
$beamCoords = 1;
} else {
die("$ifn: only beam and earth coordinates supported\n")
if (!$dta{EARTH_COORDINATES});
die("$ifn: -p requires beam-coordinate data\n")
if ($opt_p > 0);
die("$ifn: -d requires beam-coordinate data\n")
if defined($opt_d);
die("$ifn: -B requires beam-coordinate data\n")
if defined($opt_B);
}
for (my($b)=0; $b<$dta{N_BINS}; $b++) { # calc dz
$dz[$b] = $dta{DISTANCE_TO_BIN1_CENTER} + $b*$dta{BIN_LENGTH};
}
$lastGoodBin = 0;
for ($e=0; $e<=$#{$dta{ENSEMBLE}}; $e++) { # check/transform velocities
$dta{ENSEMBLE}[$e]->{NUMBER} = $e+1 if ($opt_R); # renumber ensembles
next if (defined($first_ens) &&
$dta{ENSEMBLE}[$e]->{NUMBER} < $first_ens);
$dta{ENSEMBLE}[$e]->{PITCH} -= $P{pitch_bias}
if defined($dta{ENSEMBLE}[$e]->{PITCH});
$dta{ENSEMBLE}[$e]->{ROLL} -= $P{roll_bias}
if defined($dta{ENSEMBLE}[$e]->{ROLL});
$P{first_ens} = $dta{ENSEMBLE}[$e]->{NUMBER},$fe = $e
unless defined($P{first_ens});
last if (defined($last_ens) &&
$dta{ENSEMBLE}[$e]->{NUMBER} > $last_ens);
$P{last_ens} = $dta{ENSEMBLE}[$e]->{NUMBER};
$le = $e;
die("3-beams used in ensemble #$dta{ENSEMBLE}[$e]->{NUMBER}\n")
if ($dta{ENSEMBLE}[$e]->{N_BEAMS_USED} < 4);
# die("BIT error in ensemble $dta{ENSEMBLE}[$e]->{NUMBER}\n")
# if defined($dta{ENSEMBLE}[$e]->{BUILT_IN_TEST_ERROR});
die("Low gain in ensemble #$dta{ENSEMBLE}[$e]->{NUMBER}\n")
if ($dta{ENSEMBLE}[$e]->{LOW_GAIN});
for (my($b)=0; $b<$dta{N_BINS}; $b++) {
if ($beamCoords) {
for (my($i)=0; $i<4; $i++) { # percent-good editing (-p)
if ($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][$i] < $opt_p) {
undef($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][$i]);
undef($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][$i]);
}
}
$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0] -= $P{velbias_b1} # beam-velocity biases (-B)
if defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0]);
$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][1] -= $P{velbias_b2}
if defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][1]);
$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][2] -= $P{velbias_b3}
if defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][2]);
$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3] -= $P{velbias_b4}
if defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3]);
if (defined($opt_d)) { # discard data from given beam (-d)
undef($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][$opt_d-1]);
undef($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][$opt_d-1]);
}
@{$dta{ENSEMBLE}[$e]->{BEAM_VELOCITY}[$b]} = # save beam velocities
@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]};
@{$dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b]} = # calculate v12, w12, v34, w34
velBeamToBPEarth(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{BEAM_VELOCITY}[$b]});
@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]} = # calculate earth velocities
velBeamToEarth(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{BEAM_VELOCITY}[$b]});
$dta{ENSEMBLE}[$e]->{THREE_BEAM}[$b] = $RDI_Coords::threeBeamFlag;
$three_beam[$b] += $RDI_Coords::threeBeamFlag;
unless (defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0])) { # not a valid transformation
undef(@{$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b]});
next;
}
} else { # Earth coordinates
@{$dta{ENSEMBLE}[$e]->{BEAM_VELOCITY}[$b]} = # calculate beam velocities
velEarthToBeam(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]});
@{$dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b]} = # calculate w12, w34
velBeamToBPEarth(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{BEAM_VELOCITY}[$b]});
@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]} = # correct for heading bias
velApplyHdgBias(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]});
$three_beam[$b] += ($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][0]/100 * $dta{PINGS_PER_ENSEMBLE});
unless (defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0])) { # no valid velocity
undef(@{$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b]});
next;
}
}
$dta{ENSEMBLE}[$e]->{GOOD_VEL}[$b] = 1;
$good_vels[$b]++;
$lastGoodBin = $b if ($b > $lastGoodBin);
$firstGoodEns = $e unless defined($firstGoodEns);
$lastGoodEns = $e;
}
}
unless (defined($opt_r)) {
$fe = $firstGoodEns;
$le = $lastGoodEns;
}
$P{N_ensembles} = $le - $fe + 1;
$firstBin = 0;
$lastBin = $lastGoodBin;
print( STDERR "Start : $dta{ENSEMBLE}[$fe]->{DATE} $dta{ENSEMBLE}[$fe]->{TIME}\n");
print( STDERR "End : $dta{ENSEMBLE}[$le]->{DATE} $dta{ENSEMBLE}[$le]->{TIME}\n");
printf(STDERR "Bins : %d-%d\n",$firstBin+1,$lastBin+1);
if ($beamCoords) {
printf(STDERR "3-Beam : %d %d %d %d\n",$RDI_Coords::threeBeam_1,
$RDI_Coords::threeBeam_2,
$RDI_Coords::threeBeam_3,
$RDI_Coords::threeBeam_4);
}
print(STDERR "Good/3-Beam: ");
for ($b=$firstBin; $b<=$lastBin; $b++) { # generate output
dumpBin($b,$fe,$le);
}
print(STDERR "\n");
exit(0);