.
#!/usr/bin/perl
#======================================================================
# L A D C P _ W
# doc: Fri Dec 17 18:11:13 2010
# dlm: Wed Oct 17 12:03:12 2012
# (c) 2010 A.M. Thurnherr
# uE-Info: 121 63 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# TODO:
# detection so that editing stats make sense
# - own seabed detection (P403)
# - make upcast-flag valid for yoyo casts
# - make diagnostic output 3-beam field work for Earth coordinates
# - remove water-depth from BT code, which is not really used and a bit of an outlier
# because mean and stddev are used instead of median/mad
# - read assumed ADCP soundspeed from data file, instead of assuming 1500m/s
# ? use instrument tilt in sidelobe editing?
$antsSummary = 'calculate vertical velocities from LADCP & CTD time series';
# HISTORY:
# Dec 17, 2010: - created from [mergeCTD+LADCP]
# Dec 18, 2010: - made to work
# Dec 19, 2010: - improved considerably
# Dec 20, 2010: - onward
# - BUG: depth-binning was off by 1 bin?!
# - added binning correction for instrument tilt
# Dec 21, 2010: - added -h (seafloor depth)
# Dec 22, 2010: - BUG: had not applied soundspeed-correction to w
# - debugged opt_d
# Dec 23, 2010: - continued implementation of soundspeed corrections
# Dec 24, 2010: - added winch_w, wave_w
# - removed beampair velocities from code
# Dec 25, 2010: - adapted for surface-wave correction in terms of acceleration (CTD_w_t)
# - removed elapsed_mismatch
# - removed winch_w, wave_w
# Dec 26, 2010: - made -p output layout independent of -x to avoid Makefile problems
# Dec 30, 2010: - cleaned up some
# - folded-in backscatter calculation from shear method
# - folded-in BT calculation from shear method
# Dec 31, 2010: - added weighted mean w profile to output
# Jan 2, 2010: - BUG: BT_w.mad could bomb with division by 0
# - BUG: division by zero if no valid data
# Jan 5, 2010: - adapted to allow for gaps in CTD time series
# Feb 16, 2011: - cosmetics
# Jun 22, 2011: - cosmetics
# Jun 23, 2011: - disabled error on large rms reflr w
# - added -l
# - removed CTD headers from output
# Jun 26, 2011: - added -u
# - changed package correction from acceleration to velocity, because of
# Stan's Antarctic data set where accelerations are zero but package effects are
# there
# Jul 2, 2011: - increased tilt default to 15 degrees
# Jul 3, 2011: - replaced old package-velocity correction -x code by new beamvel correction
# - removed -p
# - replaced -d by residual (diagnostics) output
# Jul 4, 2011: - saved a snapshot in [Jul_04_2011]
# - removed output of binned profile (but not calculation code, which is required for residual)
# - BUG: firstGoodEns or lastGoodEns could end up in a reflr_w gap when valid LADCP data begin
# or end with CTD in water
# - moved rarely used option -s to -k
# - added -s)kip <ens> option
# - had to very very slightly relax an assertion (by 1e-10 seconds...)
# Aug 3, 2011: - added -z)ap to truncate range, based on Stan's 2009 data set (017-019 most clearly)
# obs that final 2 velocities are outliers; could be due to color bar, should check with
# residuals
# Sep 22, 2011: - removed (commented out) CTD acceleration
# - added basename to Parameter File Matching
# Sep 23, 2011: - cosmetics (lots)
# - added -d)
# Oct 6, 2011: - removed old commented-out code (CTD acceleration & old -x)
# - renamed time-series field names
# - added uncorrected reflr w to time-series output
# Oct 10, 2011: - BUG: LADCP_w in .tds output had not been soundspeed corrected
# - removed -x
# - added "false-positives" editing filter
# Oct 11, 2011: - adapted to the new parameter file, removing a few options
# Oct 12, 2011: - added surface-layer editing
# - made %min/max_depth/elapsed more accurate
# Oct 13, 2011: - fiddled
# Oct 14, 2011: - renamed .prof fields from .N to .nsamp
# - added %output_basename
# - removed all ">" from open to allow use in pipelines
# - made <run-label> argument optional
# - replaced stdout output by $w_out
# - renamed _out to out_; output_basename to out_basename
# Oct 15, 2011: - added editWOutliers()
# - added step to remove single-ping noise
# Oct 16, 2011: - BUG: ensemble and bin numbers in output files were off by 1
# Oct 17, 2011: - added default run label
# - added out_BR
# - BUG: closed STDOUT caused problems with tee in plotting scripts
# - added %run_label
# - added PostProcess.sh
# Oct 18, 2011: - BUG: %{min,max}_ens used ensemble indices, rather than numbers
# - disabled false-positives editing (removes too much good data on 2011_IWISE)
# Oct 19, 2011: - BUG: time-series output included ensembles without valid w
# - increased post-process sleep time to 30s
# - added corr, Sv and echo_amplitude as new fields to w output
# Oct 20, 2011: - added editFarBins()
# - added downcast flag to time-series output
# - added ensemble number to LADCP-time-series output
# Oct 21, 2011: - BUG: Sv code was not enabled for uplooker data
# Oct 24, 2011: - added code to copy %lat/%lon from CTD data
# - added %ADCP_freuquency, %ADCP_blanking_distance
# - %LADCP_bin_length => %ADCP_bin_length
# Oct 26, 2011: - added $first_guess_timelag
# Oct 27, 2011: - added w12, w34, ww (pitch/roll-weighted average) to full output
# - removed ww because it apparently does not improve the solutions
# - moved editTilt() to after time-series calculation to allow
# time matching to work even with strict tilt limit
# - added correctAttitude()
# May 22, 2012: - adapted to ANTS V5
# Oct 15, 2012: - added $edit_data_hook
#--- hg commit
# - separated dc/uc time lagging
# - removed support for TLhist
# Oct 16, 2012: - added support for dc/uc only solutions
# Oct 17, 2012: - renamed $edit_data_hook to $post_merge_hook
# CTD REQUIREMENTS
# - elapsed elapsed seconds; see note below
# - depth
# - ss sound speed
# - w ddepth/dt
# - temp OPTIONAL; used for backscatter calculation (i.e. not very important)
# - %lat/%lon OPTIONAL
# NUMERICAL OPTIONS
# - the first option in the list cannot be numerical!
# - if need be, use -v 1 as a dummy option
# ELAPSED TIMES
# - there are 2 different elapsed times used in this program:
# 1) elapsed based on firstgoodens in the LADCP time series
# 2) CTD elapsed time
# - CTD elapsed time does not have to start with zero!
# - do not use the Seabird elapsed field, which is only reported to
# 3 significant digits, causing significant jitter in dt; however,
# at least up to 2010 Seabird simply calculates elapsed time by
# assuming a 24Hz sampling rate and no record drop; therefore,
# it is best to calculate elapsed time as %RECNO/24
# - the elapsed field of the output is the elapsed time from the CTD
# file; this is required in order to be able to compare the times
# from the uplooker and downlooker-derived vertical velocity
# profiles
# - the .TL files use the LADCP elapsed time in output; as a result,
# the DL/UL .TL files have different time axes!
# - as a result, a profile only starts with elapsed==0 if the CTD
# is turned on when the LADCP is already in the water
# OUTPUT NOTES
# - residuals are calculated with respect to down-/upcast medians in binned profiles
# - equivalence, assuming -o 10:
# 1) 004DL.prof dc_w depth
# 2) bindata -Sdowncast:1 -Fw.median,depth -n 20 depth 10 004DL.w
# VELOCITY AMBIGUITY ERRORS
# - quite extensive tests with DIMES US2 station 146, which has a lot of
# ambiguity velocity errors, reveal that -m 1 catches those errors
# quite nicely
# - even when the errors are not filtered with -m 1, they do not
# affect the w profiles, as long as the median bin values are used
($WCALC) = ($0 =~ m{^(.*)/[^/]*$});
($ANTS) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
($PERL_TOOLS) = (`which mkProfile` =~ m{^(.*)/[^/]*$});
require "$ANTS/ants.pl";
require "$ANTS/libstats.pl";
require "$WCALC/edit_data.pl";
require "$WCALC/time_series.pl";
require "$WCALC/time_lag.pl";
require "$WCALC/find_seabed.pl";
require "$WCALC/svel_corrections.pl";
require "$WCALC/acoustic_backscatter.pl";
require "$WCALC/bottom_tracking.pl";
require "$PERL_TOOLS/RDI_BB_Read.pl";
require "$PERL_TOOLS/RDI_Coords.pl";
use IO::Handle;
my($TINY) = 1e-16;
#------
# Usage
#------
@ARGS = @ARGV; # save opts & arguments
$antsParseHeader = 0;
&antsUsage('3:4a:b:c:e:g:h:i:k:m:n:o:qs:t:uv:w:x:',1,
'[-v)erbosity <level[1]>]',
'[-q)uick (no single-ping denoising)]',
'[require -4)-beam solutions]',
'[-c)orrelation <min[70]>] [-t)ilt <max[10deg]> [-e)rr-vel <max[0.1m/s]>]',
'[-h water <depth>]',
'[max LADCP time-series -g)ap <length[60s]>]',
'[-m)ax vertical <velocity[1m/s]>',
'[-a)djust CTD depth <by[0m]>]',
'[-i)nitial CTD time offset <guestimate> [-u)se as final]]',
'[calculate -n) <lags,lags[10,100]>] [lag -w)indow <sz,sz[240s,20s]>]',
'[require top-3) lags to account for <frac[0.6]> of all]',
'[valid LADCP -b)ins <bin[2],bin[*]>',
'[-o)utput bin <resolution[10m]>] [-k) require <min[20]> samples]',
'[e-x)ecute <perl-expr>]',
'<station-id> [run-label]');
&antsUsageError() if ($opt_u && !defined($opt_i));
&antsCardOpt(\$opt_s,0); # skip # initial ensembles
$STN = &antsCardArg(); # station id
if (@ARGV) { # run label
$RUN = $ARGV[0];
shift;
} else {
$RUN = 'default';
}
#-----------------------------
# Handle Processing Parameters
#-----------------------------
require "$WCALC/defaults.pl"; # load default/option parameters
require "$processing_param_file"; # load processing parameters
eval($opt_x) if defined($opt_x); # eval cmd-line expression to override anything
croak("$0: \$out_basename undefined\n") # all plotting routines use this
unless defined($out_basename);
&antsAddParams('out_basename',$out_basename);
&antsAddParams('run_label',$RUN);
$out_w = ">$out_w" unless ($out_w =~ /^$|^\s*\|/); # make sure O_CREAT is set unless it's a pipe
$out_BR = ">$out_BR" unless ($out_BR =~ /^$|^\s*\|/);
$out_TL = ">$out_TL" unless ($out_TL =~ /^$|^\s*\|/);
$out_LADCPtis = ">$out_LADCPtis" unless ($out_LADCPtis =~ /^$|^\s*\|/);
$out_profile = ">$out_profile" unless ($out_profile =~ /^$|^\s*\|/);
$out_timeseries = ">$out_timeseries" unless ($out_timeseries =~ /^$|^\s*\|/);
#----------------------------------------------------------------------
# Screen Logging
# - warning levels:
# 0 probably unimportant, e.g. nonsensical parameters that probably won't affect solution
# 1 may be somewhat important
# 2 important
#----------------------------------------------------------------------
open(LOGF,">$out_log") || croak("$out_log: $!\n")
if defined($out_log);
LOGF->autoflush(1);
sub progress(@)
{
printf(LOGF @_) if defined($out_log);
printf(STDERR @_) if ($verbosity_level > 1);
}
sub info(@)
{
printf(LOGF "\t",@_) if defined($out_log);
print(STDERR ($verbosity_level > 1) ? "\t" : "$LADCP_file: ");
printf(STDERR @_) if ($verbosity_level > 0);
}
sub warning(@)
{
my($lvl,@msg) = @_;
if (defined($out_log)) {
print(LOGF "\nWARNING (level $lvl): ");
printf(LOGF @msg);
print(LOGF "\n");
}
return if ($verbosity_level == 0);
print(STDERR "\n") if ($verbosity_level > 1);
print(STDERR "WARNING (level $lvl): ");
printf(STDERR @msg);
print(STDERR "\n") if ($verbosity_level > 1);
}
sub debugmsg(@)
{ printf(STDERR @_) if ($verbosity_level > 2); }
#----------------
# Read LADCP data
#----------------
progress("Reading LADCP data from <$LADCP_file>...\n");
readData($LADCP_file,\%LADCP);
progress("\t%d ensembles\n",scalar(@{$LADCP{ENSEMBLE}}));
croak("$LADCP_file: not enough LADCP bins ($LADCP{N_BINS}) for choice of -r\n")
unless ($LADCP{N_BINS} >= $refLr_lastBin);
croak("$0: first reference-layer bin outside valid range\n")
unless ($refLr_firstBin>=1 && $refLr_firstBin<=$LADCP{N_BINS});
croak("$0: last reference-layer bin outside valid range\n")
unless ($refLr_lastBin>=1 && $refLr_lastBin<=$LADCP{N_BINS});
croak("$0: first reference-layer bin > last reference-layer bin\n")
unless ($refLr_firstBin <= $refLr_lastBin);
$LADCP_lastBin = $LADCP{N_BINS}-1
if ($LADCP_lastBin eq '*');
croak("$0: first valid LADCP bin outside valid range\n")
unless ($LADCP_firstBin>=1 && $LADCP_firstBin<=$LADCP{N_BINS});
croak("$0: last valid LADCP bin outside valid range\n")
unless ($LADCP_lastBin>=1 && $LADCP_lastBin<=$LADCP{N_BINS});
croak("$0: first valid LADCP bin > last valid LADCP bin\n")
unless ($LADCP_firstBin <= $LADCP_lastBin);
warning(0,"first reference-layer bin < first valid LADCP bin\n")
unless ($refLr_firstBin >= $LADCP_firstBin);
warning(0,"last reference-layer bin > last valid LADCP bin\n")
unless ($refLr_lastBin <= $LADCP_lastBin);
warning(1,"if at all, bin 1 should not be used for short blank-after-transmit values\n")
if ($LADCP{BLANKING_DISTANCE}<$LADCP{BIN_LENGTH} && $refLr_firstBin==1);
&antsAddParams('ADCP_bin_length',$LADCP{BIN_LENGTH},
'ADCP_frequency',$LADCP{BEAM_FREQUENCY},
'ADCP_blanking_distance',$LADCP{BLANKING_DISTANCE});
#------------------------------------------------------------
# Edit beam-velocity data
# 1) correlation threshold
#------------------------------------------------------------
if ($LADCP{BEAM_COORDINATES}) {
progress("Editing beam-velocity data...\n");
$nvv = $cte = 0;
for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
correctAttitude($ens,$pitch_bias,$roll_bias,$heading_bias);
$nvv += countValidBeamVels($ens);
$cte += editCorr($ens,$min_correlation);
# $pte += editTilt($ens,$max_tilt);
}
croak("$LADCP_file: no valid data\n") unless ($nvv > 0);
progress("\tcorrelation threshold (-c %d counts): %d velocites removed (%d%% of total)\n",$min_correlation,$cte,round(100*$cte/$nvv));
# progress("\tattitude threshold (-t %d deg): %d velocites removed (%d%% of total)\n",$max_tilt,$pte,round(100*$pte/$nvv));
} else {
progress("Editing velocity data...\n");
$nvv = $cte = 0;
for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
$nvv += countValidBeamVels($ens);
$cte += editCorr_Earthcoords($ens,$min_correlation);
# $pte += editTilt($ens,$max_tilt);
}
croak("$LADCP_file: no valid data\n") unless ($nvv > 0);
progress("\tcorrelation threshold (-c %d counts): %d velocites removed (%d%% of total)\n",$min_correlation,$cte,round(100*$cte/$nvv));
# progress("\tattitude threshold (-t %d deg): %d velocites removed (%d%% of total)\n",$max_tilt,$pte,round(100*$pte/$nvv));
}
#-------------------------------------------------------------------
# Calculate earth velocities
# - this is done for all bins (not just valid ones), to allow
# useless possibility that invalid bins are used for reflr calcs
# - also calculate separate beam-pair velocities
# - the UNEDITED velocities are saved for the BT calculations
# (W is required, U & V are only used for stats that have not
# been very useful so far)
#-------------------------------------------------------------------
if ($LADCP{BEAM_COORDINATES}) {
my($dummy);
progress("Calculating earth-coordinate velocities...\n");
progress("\tdiscarding velocities from beam $discard_velocities_from_beam\n")
if ($discard_velocities_from_beam);
$nvw = 0;
for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
undef($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$discard_velocities_from_beam-1])
if ($discard_velocities_from_beam);
($LADCP{ENSEMBLE}[$ens]->{U}[$bin],
$LADCP{ENSEMBLE}[$ens]->{V}[$bin],
$LADCP{ENSEMBLE}[$ens]->{W}[$bin],
$LADCP{ENSEMBLE}[$ens]->{ERRVEL}[$bin]) =
velInstrumentToEarth(\%LADCP,$ens,
velBeamToInstrument(\%LADCP,
@{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]}));
if (defined($LADCP{ENSEMBLE}[$ens]->{W}[$bin])) {
$per_bin_nsamp[$bin]++;
$nvw++;
}
$LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH} =
gimbal_pitch($LADCP{ENSEMBLE}[$ens]->{PITCH},$LADCP{ENSEMBLE}[$ens]->{ROLL});
$LADCP{ENSEMBLE}[$ens]->{U_UNEDITED}[$bin] = $LADCP{ENSEMBLE}[$ens]->{U}[$bin];
$LADCP{ENSEMBLE}[$ens]->{V_UNEDITED}[$bin] = $LADCP{ENSEMBLE}[$ens]->{V}[$bin];
$LADCP{ENSEMBLE}[$ens]->{W_UNEDITED}[$bin] = $LADCP{ENSEMBLE}[$ens]->{W}[$bin];
($dummy,$LADCP{ENSEMBLE}[$ens]->{W12}[$bin],
$dummy,$LADCP{ENSEMBLE}[$ens]->{W34}[$bin]) =
velBeamToBPEarth(\%LADCP,$ens,@{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]});
}
}
progress("\t$nvw valid velocities in bins $LADCP_firstBin-$LADCP_lastBin\n");
progress("\t3-beam solutions : $RDI_Coords::threeBeam_1 " .
"$RDI_Coords::threeBeam_2 " .
"$RDI_Coords::threeBeam_3 " .
"$RDI_Coords::threeBeam_4\n")
unless ($opt_4);
} else { # Earth Coordinates
progress("Counting valid vertical velocities...\n");
$nvw = 0;
for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
($LADCP{ENSEMBLE}[$ens]->{U}[$bin],
$LADCP{ENSEMBLE}[$ens]->{V}[$bin],
$LADCP{ENSEMBLE}[$ens]->{W}[$bin],
$LADCP{ENSEMBLE}[$ens]->{ERRVEL}[$bin]) = @{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]};
if (defined($LADCP{ENSEMBLE}[$ens]->{W}[$bin])) {
$per_bin_nsamp[$bin]++;
$nvw++;
}
$LADCP{ENSEMBLE}[$ens]->{U_UNEDITED}[$bin] = $LADCP{ENSEMBLE}[$ens]->{U}[$bin];
$LADCP{ENSEMBLE}[$ens]->{V_UNEDITED}[$bin] = $LADCP{ENSEMBLE}[$ens]->{V}[$bin];
$LADCP{ENSEMBLE}[$ens]->{W_UNEDITED}[$bin] = $LADCP{ENSEMBLE}[$ens]->{W}[$bin];
}
}
progress("\t$nvw valid velocities in bins $LADCP_firstBin-$LADCP_lastBin\n");
}
croak("$LADCP_file: insufficient valid velocities\n") unless ($nvw > 1000);
#----------------------------------------------
# STEP: Edit earth-coordinate -velocity data
# 1) error-velocity threshold
# 2) vertical-velocity outliers
# 3) truncate range by deleting farthest valid velocities (disbled by default)
#----------------------------------------------
progress("Editing earth-coordinate velocity data...\n");
$evrm = $trrm = $worm = 0;
for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
$evrm += editErrVel($ens,$max_allowed_errvel);
$worm += editWOutliers($ens,$per_ens_outliers_mad_limit);
$trrm += editTruncRange($ens,$truncate_farthest_valid_bins) if ($truncate_farthest_valid_bins > 0);
}
progress("\terror-velocity threshold (-e %.1f m/s): %d velocites removed (%d%% of total in bins $LADCP_firstBin-$LADCP_lastBin)\n",
$max_allowed_errvel,$evrm,round(100*$evrm/$nvw));
progress("\tvertical-velocity outliers ($per_ens_outliers_mad_limit x mad): %d velocites removed (%d%% of total in bins $LADCP_firstBin-$LADCP_lastBin)\n",
$worm,round(100*$worm/$nvw));
progress("\trange truncation (-z %d bins): %d velocites removed (%d%% of total in bins $LADCP_firstBin-$LADCP_lastBin)\n",
$truncate_farthest_valid_bins,$trrm,round(100*$trrm/$nvw))
if ($truncate_farthest_valid_bins > 0);
#----------------------------
# Calculate LADCP time series
#----------------------------
progress("Calculating LADCP time-series...\n");
($firstGoodEns,$lastGoodEns,$LADCP_atbottom,$LADCP_w_gap_time) =
calcLADCPts(\%LADCP,$opt_s,$refLr_firstBin,$refLr_lastBin,$max_LADCP_reflr_vel_gap);
croak("$LADCP_file: no good ensembles\n")
unless defined($firstGoodEns) && ($lastGoodEns-$firstGoodEns > 0);
my($cast_duration) = $LADCP{ENSEMBLE}[$lastGoodEns]->{ELAPSED} -
$LADCP{ENSEMBLE}[$firstGoodEns]->{ELAPSED};
croak("$0: implausibly short cast ($cast_duration seconds)\n")
unless ($cast_duration > 600);
$LADCP{MEAN_DT} = $cast_duration / ($lastGoodEns-$firstGoodEns-1);
progress("\tStart of cast : %s (#%5d)\n",
$LADCP{ENSEMBLE}[$firstGoodEns]->{TIME},
$LADCP{ENSEMBLE}[$firstGoodEns]->{NUMBER});
progress("\tBottom of cast : %s (#%5d) @ dz=%6dm\n",
$LADCP{ENSEMBLE}[$LADCP_atbottom]->{TIME},
$LADCP{ENSEMBLE}[$LADCP_atbottom]->{NUMBER},
$LADCP{ENSEMBLE}[$LADCP_atbottom]->{DEPTH});
progress("\tEnd of cast : %s (#%5d) @ dz=%6dm\n",
$LADCP{ENSEMBLE}[$lastGoodEns]->{TIME},
$LADCP{ENSEMBLE}[$lastGoodEns]->{NUMBER},
$LADCP{ENSEMBLE}[$lastGoodEns]->{DEPTH});
progress("\tCast duration : %.1f hours (pinging for %.1f hours)\n",
$cast_duration / 3600,
($LADCP{ENSEMBLE}[$#{$LADCP{ENSEMBLE}}]->{UNIX_TIME} -
$LADCP{ENSEMBLE}[0]->{UNIX_TIME}) / 3600);
progress("\tMean ping interval: %.1f seconds\n",$LADCP{MEAN_DT});
if (defined($out_LADCPtis)) {
progress("Writing LADCP-timeseries to <$out_LADCPtis>...\n");
@antsNewLayout = ('ens','elapsed','reflr_w','reflr_w.stddev','reflr_w.nsamp','depth');
open(STDOUT,"$out_LADCPtis") || croak("$out_LADCPtis: $!\n");
for (my($ens)=$firstGoodEns; $ens<=$lastGoodEns; $ens++) {
&antsOut($LADCP{ENSEMBLE}[$ens]->{NUMBER},
$LADCP{ENSEMBLE}[$ens]->{ELAPSED},
$LADCP{ENSEMBLE}[$ens]->{REFLR_W},
$LADCP{ENSEMBLE}[$ens]->{REFLR_W_STDDEV},
$LADCP{ENSEMBLE}[$ens]->{REFLR_W_NSAMP},
$LADCP{ENSEMBLE}[$ens]->{DEPTH});
}
&antsOut('EOF'); open(STDOUT,">&2");
}
#----------------------------------------------------------------------
# More editing
# - this requires ${first,last}GoodEns to be known
# - TILT field is set as a side-effect
#----------------------------------------------------------------------
progress("Editing additional earth-coordinate velocity data...\n");
my($first_bad_bin);
for (my($bin)=$LADCP_firstBin-1; $bin<$LADCP_lastBin-1; $bin++) {
next unless ($per_bin_nsamp[$bin]/($lastGoodEns-$firstGoodEns) < $per_bin_valid_frac_lim);
$first_bad_bin = $bin;
last;
}
if (defined($first_bad_bin)) {
$fprm = $pte = 0;
for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
$pte += editTilt($ens,$max_tilt);
$fprm += editFarBins($ens,$first_bad_bin) if defined($first_bad_bin);
}
progress("\tattitude threshold (max_tilt = %d deg): %d velocites removed (%d%% of total)\n",
$max_tilt,$pte,round(100*$pte/$nvv));
progress("\tvelocities beyond bin $first_bad_bin (<%d%% valid values): %d velocites removed (%d%% of total in bins $LADCP_firstBin-$LADCP_lastBin)\n",
round(100*$per_bin_valid_frac_lim),$fprm,round(100*$fprm/$nvw));
}
#--------------
# Read CTD data
#--------------
progress("Reading CTD data from <$CTD_file>...\n");
open(STDIN,$CTD_file) || croak("$CTD_file: $!\n");
croak("$CTD_file: no data\n") unless (&antsIn());
undef($antsOldHeaders);
($CTD_elapsed,$CTD_depth,$CTD_svel,$CTD_w) = &fnr('elapsed','depth','ss','w');
$CTD_temp = &fnrNoErr('temp');
&antsAddParams('lat',$P{lat}) if defined($P{lat});
&antsAddParams('lon',$P{lon}) if defined($P{lon});
$CTD_maxdepth = -1;
do {
croak("$0: cannot deal with non-numeric CTD elapsed time\n")
unless &antsNumbers($CTD_elapsed);
push(@{$CTD{ELAPSED}},$ants_[0][$CTD_elapsed]);
push(@{$CTD{DEPTH}}, $ants_[0][$CTD_depth]);
push(@{$CTD{SVEL}}, $ants_[0][$CTD_svel]);
push(@{$CTD{W}}, $ants_[0][$CTD_w]);
push(@{$CTD{TEMP}}, $ants_[0][$CTD_temp]) if defined($CTD_temp);
if ($ants_[0][$CTD_depth] > $CTD_maxdepth) {
$CTD_maxdepth = $ants_[0][$CTD_depth];
$CTD_atbottom = $#{$CTD{DEPTH}};
}
} while (&antsIn());
$CTD{DT} = ($CTD{ELAPSED}[$#{$CTD{ELAPSED}}] - $CTD{ELAPSED}[0]) / $#{$CTD{ELAPSED}};
progress("\t%d scans at %.1fHz\n",scalar(@{$CTD{DEPTH}}),1/$CTD{DT});
croak("$0: CTD start depth must be numeric\n")
unless numberp($CTD{DEPTH}[0]);
$CTD_neg_press_offset = -1 * round($CTD{DEPTH}[0])
if ($CTD{DEPTH}[0] < -$CTD_neg_press_offset);
if ($CTD_neg_press_offset > 0) {
progress("\tadding $CTD_neg_press_offset dbar offset to pressure data\n");
for (my($i)=0; $i<@{$CTD{DEPTH}}; $i++) { $CTD{DEPTH}[$i] += $CTD_neg_press_offset; }
$CTD_maxdepth += $CTD_neg_press_offset;
}
progress("\tstart depth = %.1fm\n",$CTD{DEPTH}[0]);
progress("\tmax depth = %dm (# $CTD_atbottom)\n",$CTD_maxdepth);
#--------------------------------------------------------------------
# Construct sound-speed correction profile from CTD 1Hz downcast data
# very simple algorithm that stores the last value found
# in each 1m bin
#--------------------------------------------------------------------
progress("Constructing sound-speed correction profile\n");
my($scans_per_sec) = int(1/$CTD{DT}+0.5);
for (my($s)=0; $s<=$CTD_atbottom; $s+=$scans_per_sec) {
next unless ($CTD{DEPTH}[$s] >= 0 && numberp($CTD{SVEL}[$s]));
$sVelProf[int($CTD{DEPTH}[$s])] = $CTD{SVEL}[$s];
}
#-------------------
# Determine time lag
#-------------------
if (defined($first_guess_timelag)) {
progress("Setting initial time lag...\n");
$CTD{TIME_LAG} = $first_guess_timelag;
progress("\t-i => elapsed(CTD) ~ elapsed(LADCP) + %.1fs\n",$CTD{TIME_LAG});
} else {
progress("Guestimating time lag...\n");
my($CTD_10pct_down) = 0;
my($LADCP_10pct_down) = $firstGoodEns;
$CTD_10pct_down++
until ($CTD{DEPTH}[$CTD_10pct_down]-$CTD{DEPTH}[0] >= 0.1*($CTD_maxdepth-$CTD{DEPTH}[0]));
$LADCP_10pct_down++
until ($LADCP{ENSEMBLE}[$LADCP_10pct_down]->{DEPTH} >= 0.1*$LADCP{ENSEMBLE}[$LADCP_atbottom]->{DEPTH});
$CTD{TIME_LAG} =
$CTD{ELAPSED}[$CTD_10pct_down] - $LADCP{ENSEMBLE}[$LADCP_10pct_down]->{ELAPSED};
progress("\telapsed(dz(CTD)=%.1fm) ~ elapsed(dz(LADCP)=%.1fm) + %.1fs\n",
$CTD{DEPTH}[$CTD_10pct_down]-$CTD{DEPTH}[0],$LADCP{ENSEMBLE}[$LADCP_10pct_down]->{DEPTH},$CTD{TIME_LAG});
}
if ($opt_u) {
progress("\tskipping time lagging (-u)\n");
} else {
$CTD{TIME_LAG} = # stage 1: dc/uc combined
calc_lag($number_of_timelag_windows[0],$length_of_timelag_windows[0],int(1/$CTD{DT}+0.5),0);
croak("$0: Cannot proceed without valid lag!\n") unless defined($CTD{TIME_LAG});
progress("\telapsed(CTD) ~ elapsed(LADCP) + %.2fs\n",$CTD{TIME_LAG});
$CTD{TIME_LAG_DC} = calc_lag($number_of_timelag_windows[1],$length_of_timelag_windows[1],1,1);
progress("\tdowncast: elapsed(CTD) = elapsed(LADCP) + %.2fs\n",$CTD{TIME_LAG_DC})
if defined($CTD{TIME_LAG_DC});
$CTD{TIME_LAG_UC} = calc_lag($number_of_timelag_windows[1],$length_of_timelag_windows[1],1,-1);
progress("\tupcast : elapsed(CTD) = elapsed(LADCP) + %.2fs\n",$CTD{TIME_LAG_UC})
if defined($CTD{TIME_LAG_UC});
croak("$0: Cannot proceed without either down- or upcast time-lag!\n")
unless defined($CTD{TIME_LAG_DC}) || defined($CTD{TIME_LAG_UC});
$CTD_time_lag_dc = $CTD{TIME_LAG_DC}; # if one lag is missing ...
$CTD_time_lag_uc = $CTD{TIME_LAG_UC}; # use other for entire ...
$CTD_time_lag_dc = $CTD_time_lag_uc unless defined($CTD_time_lag_dc); # cast
$CTD_time_lag_uc = $CTD_time_lag_dc unless defined($CTD_time_lag_uc);
}
&antsAddParams('CTD_time_lag_dc',$CTD{TIME_LAG_DC},'CTD_time_lag_uc',$CTD{TIME_LAG_UC});
undef($CTD{TIME_LAG}); # to avoid confusion
#------------------------------------------------
# Merge CTD with LADCP data
# - after this step, reflr w is sound-speed corrected!!!
#------------------------------------------------
progress("Merging CTD with LADCP data...\n");
for (my($skipped)=0,my($ens)=$firstGoodEns; $ens<=$lastGoodEns; $ens++) {
my($scan) = ($ens < $LADCP_atbottom) ?
int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED} + $CTD_time_lag_dc - $CTD{ELAPSED}[0]) / $CTD{DT} + 0.5) :
int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED} + $CTD_time_lag_uc - $CTD{ELAPSED}[0]) / $CTD{DT} + 0.5);
if ($scan <= 0) { # NB: must be <=, rather than <, or assertion below sometimes fails
$skipped++;
$firstGoodEns = $ens+1;
next;
}
if ($skipped > 0) {
$firstGoodEns++,$skipped++,next # in gap
unless defined($LADCP{ENSEMBLE}[$firstGoodEns]->{REFLR_W});
info("$skipped initial LADCP ensembles skipped because CTD data begin with LADCP in water\n");
$skipped = 0;
}
if ($scan > $#{$CTD{ELAPSED}}) {
while (!defined($LADCP{ENSEMBLE}[$ens-1]->{REFLR_W})) { $ens--; } # in gap
info(sprintf("%d final LADCP ensembles skipped because CTD data end with LADCP in water\n",
$lastGoodEns-$ens+1));
$lastGoodEns = $ens-1;
last;
}
if ($ens < $LADCP_atbottom) {
die("assertion failed!\n" .
"\ttest: abs($LADCP{ENSEMBLE}[$ens]->{ELAPSED} + $CTD_time_lag_dc - $CTD{ELAPSED}[$scan]) <= $CTD{DT}/2\n" .
"\tens = $ens, scan = $scan\n" .
sprintf("\tadjusted LADCP time = %f\n",$LADCP{ENSEMBLE}[$ens]->{ELAPSED} + $CTD_time_lag_dc) .
sprintf("\tCTD($scan) time = %f\n",$CTD{ELAPSED}[$scan]) .
"=> Did you use SeaBird elapsed time? Don't!"
) unless (abs($LADCP{ENSEMBLE}[$ens]->{ELAPSED} + $CTD_time_lag_dc - $CTD{ELAPSED}[$scan]) <= $CTD{DT}/2 + 1e-10);
} else {
die("assertion failed!\n" .
"\ttest: abs($LADCP{ENSEMBLE}[$ens]->{ELAPSED} + $CTD_time_lag_uc - $CTD{ELAPSED}[$scan]) <= $CTD{DT}/2\n" .
"\tens = $ens, scan = $scan\n" .
sprintf("\tadjusted LADCP time = %f\n",$LADCP{ENSEMBLE}[$ens]->{ELAPSED} + $CTD_time_lag_uc) .
sprintf("\tCTD($scan) time = %f\n",$CTD{ELAPSED}[$scan]) .
"=> Did you use SeaBird elapsed time? Don't!"
) unless (abs($LADCP{ENSEMBLE}[$ens]->{ELAPSED} + $CTD_time_lag_uc - $CTD{ELAPSED}[$scan]) <= $CTD{DT}/2 + 1e-10);
}
$LADCP{ENSEMBLE}[$ens]->{CTD_ELAPSED} = $CTD{ELAPSED}[$scan]; # elapsed field for output
if (defined($LADCP{ENSEMBLE}[$ens]->{REFLR_W}) # not a gap
&& numberp($CTD{DEPTH}[$scan])) {
$LADCP{ENSEMBLE}[$ens]->{REFLR_W_NOSSCORR} = $LADCP{ENSEMBLE}[$ens]->{REFLR_W};
$LADCP{ENSEMBLE}[$ens]->{REFLR_W} *= $CTD{SVEL}[$scan]/1500; # correct for sound-speed variations at source
croak(sprintf("\n$0: negative depth (%.1fm) in CTD file at elapsed(CTD) = %.1fs (use -a?)\n",
$CTD{DEPTH}[$scan],$CTD{ELAPSED}[$scan]))
unless ($CTD{DEPTH}[$scan] >= 0);
$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH} = $CTD{DEPTH}[$scan];
$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN} = $scan;
my($reflr_ocean_w) = $LADCP{ENSEMBLE}[$ens]->{REFLR_W} - $CTD{W}[$scan];
if (abs($reflr_ocean_w) <= $max_allowed_w) {
$sumWsq += &SQR($reflr_ocean_w);
$nWsq++;
if ($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH} > 100 &&
$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH} < $LADCP{ENSEMBLE}[$LADCP_atbottom]->{CTD_DEPTH}-100) {
$sumWsqI += &SQR($reflr_ocean_w);
$nWsqI++;
}
} else {
undef($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH}); # DON'T USE THIS ENSEMBLE LATER
}
} else{
undef($LADCP{ENSEMBLE}[$ens]->{REFLR_W}); # don't output in time-series file
undef($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH}); # old DEPTH from calcLADCPts() must be removed
}
}
if ($nWsq > 0 && $nWsqI > 0) {
&antsAddParams('rms_w_reflr_err',sqrt($sumWsq/$nWsq),'rms_w_reflr_err_interior',sqrt($sumWsqI/$nWsqI));
progress("\t%.2f cm/s rms reference-layer w_ocean, %.2f cm/s away from boundaries\n",
100*sqrt($sumWsq/$nWsq),100*sqrt($sumWsqI/$nWsqI));
warning(0,"%.2f cm/s reference-layer w_ocean away from boundaries\n",100*sqrt($sumWsqI/$nWsqI))
if (sqrt($sumWsqI/$nWsqI) > 0.05);
# croak("$0: rms reference-layer w_ocean is too large\n")
# unless (sqrt($sumWsqI/$nWsqI) < 0.07);
} elsif ($nWsq > 0) {
&antsAddParams('rms_w_reflr_err',sqrt($sumWsq/$nWsq),'rms_w_reflr_err_interior',nan);
progress("\t%.2f cm/s rms reference-layer w_ocean\n",100*sqrt($sumWsq/$nWsq));
} else {
croak("$0: no valid vertical velocities\n");
}
#----------------------------------------------------------------------
# Calculate Volume Scattering Coefficients
#----------------------------------------------------------------------
progress("Calculating volume-scattering coefficients...\n");
calc_backscatter_profs($firstGoodEns,$lastGoodEns);
#----------------------------------------------------------------------------
# Remove data contaminated by sidelobe reflection from seabed and sea surface
#----------------------------------------------------------------------------
if ($LADCP{ENSEMBLE}[$LADCP_atbottom]->{XDUCER_FACING_DOWN}) {
&antsAddParams('ADCP_orientation','downlooker');
if (numberp($opt_h)) {
progress("Setting water depth (-h)\n");
$water_depth = $opt_h;
$sig_water_depth = 0;
} else {
progress("Finding seabed...\n");
($water_depth,$sig_water_depth) =
find_backscatter_seabed($LADCP{ENSEMBLE}[$LADCP_atbottom]->{CTD_DEPTH});
($water_depth_BT,$sig_water_depth_BT) =
find_seabed(\%LADCP,$LADCP_atbottom,$LADCP{BEAM_COORDINATES});
if (defined($water_depth_BT)) {
my($dd) = abs($water_depth_BT - $water_depth);
warning(2,sprintf("Large RDI vs. own water-depth difference (%.1fm)\n",$dd))
if ($dd > 5);
}
}
&antsAddParams('water_depth',$water_depth,'water_depth.sig',$sig_water_depth);
if (defined($water_depth)) {
if (defined($water_depth_BT)) {
progress("\t%.1f(%.1f) m water depth (%.1f(%.1f)m from BT)\n",
$water_depth,$sig_water_depth,$water_depth_BT,$sig_water_depth_BT);
} else {
progress("\t%.1f(%.1f) m water depth\n",$water_depth,$sig_water_depth);
}
warning(1,sprintf("large uncertainty in water-depth estimation (%.1fm)\n",$sig_water_depth))
if ($sig_water_depth > $LADCP{BIN_LENGTH});
progress("Editing data to remove sidelobe interference from seabed...\n");
($nvrm,$nerm) = editSideLobes($firstGoodEns,$lastGoodEns,$water_depth);
progress("\t$nvrm velocities from $nerm ensembles removed\n");
} else {
info("no seabed found in backscatter profiles --- no sidelobe editing done\n");
}
} else {
&antsAddParams('ADCP_orientation','uplooker');
progress("Editing data to remove sidelobe interference from sea surface...\n");
($nvrm,$nerm) = editSideLobes($firstGoodEns,$lastGoodEns,undef);
progress("\t$nvrm velocities from $nerm ensembles removed\n");
}
#----------------------------------------------------------------------
# Data Editing after LADCP and CTD data have been merged
# 1) surface layer editing
# 2) user-supplied $edit_data_hook
#----------------------------------------------------------------------
progress("Removing data from instrument at surface...\n");
$nerm = editSurfLayer($firstGoodEns,$lastGoodEns,$surface_layer_depth);
progress("\t$nerm ensembles removed\n");
if (defined($post_merge_hook)) {
progress("Executing user-supplied \$post_merge_hook...\n");
&{$post_merge_hook}($firstGoodEns,$lastGoodEns);
}
#---------------------------------------------------------------------------
# Depth-bin LADCP velocity data
#
# NOTES:
# 1) ensemble and bin numbers are saved for maximum flexibility
# 2) only ensemble/bins with valid vertical velocities are saved
# 3) applying the full soundspeed correction to w is pointless in
# practice, but hey!, CPU cycles are cheap; [in a cast in the Gulf of Mexico
# which has fairly pronounce soundspeed gradients, the max value of Kn
# is 1.00004160558372, which gives rise to a correction of less than 0.2mm/s
# at a winch+wave speed of 3m/s....]
# 4) as far as I can tell, the soundspeed correction for bin length also
# has only a minute effect
#---------------------------------------------------------------------------
progress("Binning velocities...\n");
my($min_depth) = 9e99;
my($max_depth) = 0;
progress("\tdowncast...\n");
for ($ens=$firstGoodEns; $ens<$LADCP_atbottom; $ens++) { # downcast
unless (numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH})) {
$firstGoodEns++ if ($ens == $firstGoodEns); # start has been edited away
next;
}
my(@bindepth) = calc_binDepths($ens);
for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
$min_depth = $bindepth[$bin] if ($bindepth[$bin] < $min_depth);
$max_depth = $bindepth[$bin] if ($bindepth[$bin] > $max_depth);
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin] =
sscorr_w($LADCP{ENSEMBLE}[$ens]->{W}[$bin],
$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH},
$bindepth[$bin]);
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin] =
sscorr_w($LADCP{ENSEMBLE}[$ens]->{W12}[$bin],
$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH},
$bindepth[$bin])
if numberp($LADCP{ENSEMBLE}[$ens]->{W12}[$bin]);
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin] =
sscorr_w($LADCP{ENSEMBLE}[$ens]->{W34}[$bin],
$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH},
$bindepth[$bin])
if numberp($LADCP{ENSEMBLE}[$ens]->{W34}[$bin]);
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_W}[$bin] =
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin] + $CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}];
my($bi) = $bindepth[$bin]/$output_bin_size;
push(@{$DNCAST{ENSEMBLE}[$bi]},$ens);
push(@{$DNCAST{ELAPSED}[$bi]},$CTD{ELAPSED}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]);
push(@{$DNCAST{CTD_W}[$bi]},$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]);
push(@{$DNCAST{BIN}[$bi]},$bin);
push(@{$DNCAST{DEPTH}[$bi]},$bindepth[$bin]);
push(@{$DNCAST{W}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin]);
}
}
for (my($bi)=0; $bi<=$#{$DNCAST{ENSEMBLE}}; $bi++) { # bin data into profile
$DNCAST{MEAN_DEPTH}[$bi] = avg(@{$DNCAST{DEPTH}[$bi]});
$DNCAST{MEAN_ELAPSED}[$bi] = avg(@{$DNCAST{ELAPSED}[$bi]});
$DNCAST{MEDIAN_W}[$bi] = median(@{$DNCAST{W}[$bi]});
$DNCAST{MAD_W}[$bi] = mad2($DNCAST{MEDIAN_W}[$bi],@{$DNCAST{W}[$bi]});
$DNCAST{N_SAMP}[$bi] = @{$DNCAST{W}[$bi]};
}
progress("\tupcast...\n"); # upcast
my($realLastGoodEns);
for ($ens=$LADCP_atbottom; $ens<=$lastGoodEns; $ens++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
$realLastGoodEns = $ens;
my(@bindepth) = calc_binDepths($ens);
for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
$min_depth = $bindepth[$bin] if ($bindepth[$bin] < $min_depth);
$max_depth = $bindepth[$bin] if ($bindepth[$bin] > $max_depth);
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin] =
sscorr_w($LADCP{ENSEMBLE}[$ens]->{W}[$bin],
$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH},
$bindepth[$bin]);
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin] =
sscorr_w($LADCP{ENSEMBLE}[$ens]->{W12}[$bin],
$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH},
$bindepth[$bin])
if numberp($LADCP{ENSEMBLE}[$ens]->{W12}[$bin]);
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin] =
sscorr_w($LADCP{ENSEMBLE}[$ens]->{W34}[$bin],
$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH},
$bindepth[$bin])
if numberp($LADCP{ENSEMBLE}[$ens]->{W34}[$bin]);
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_W}[$bin] =
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin] + $CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}];
my($bi) = $bindepth[$bin]/$output_bin_size;
push(@{$UPCAST{ENSEMBLE}[$bi]},$ens);
push(@{$UPCAST{ELAPSED}[$bi]},$CTD{ELAPSED}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]);
push(@{$UPCAST{CTD_W}[$bi]},$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]);
push(@{$UPCAST{BIN}[$bi]},$bin);
push(@{$UPCAST{DEPTH}[$bi]},$bindepth[$bin]);
push(@{$UPCAST{W}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin]);
}
}
for (my($bi)=0; $bi<=$#{$UPCAST{ENSEMBLE}}; $bi++) {
$UPCAST{MEAN_DEPTH}[$bi] = avg(@{$UPCAST{DEPTH}[$bi]});
$UPCAST{MEAN_ELAPSED}[$bi] = avg(@{$UPCAST{ELAPSED}[$bi]});
$UPCAST{MEDIAN_W}[$bi] = median(@{$UPCAST{W}[$bi]});
$UPCAST{MAD_W}[$bi] = mad2($UPCAST{MEDIAN_W}[$bi],@{$UPCAST{W}[$bi]});
$UPCAST{N_SAMP}[$bi] = @{$UPCAST{W}[$bi]};
}
&antsAddParams('min_depth',$min_depth,'max_depth',$max_depth, # plot area
'min_ens',$LADCP{ENSEMBLE}[$firstGoodEns]->{NUMBER},
'max_ens',$LADCP{ENSEMBLE}[$realLastGoodEns]->{NUMBER},
'min_elapsed',$LADCP{ENSEMBLE}[$firstGoodEns]->{CTD_ELAPSED},
'max_elapsed',$LADCP{ENSEMBLE}[$realLastGoodEns]->{CTD_ELAPSED});
#------------------------------------------------------------------------------------------------------
# remove single-ping noise
# - based on observation that without this step raw velocities tend to be strongly vertically striped
# => residual CTD/LADCP mismatch errors & ADCP single-ping errors
# - remove ensemble-median residual from all vertical velocities; new residual calculated below
# - improves solution but can introduce errors in ensembles with strong outliers => median instead of
# mean, and outlier filter above
#------------------------------------------------------------------------------------------------------
unless ($opt_q) {
progress("Removing single-ping residual noise...\n");
my($fe) = defined($CTD{TIME_LAG_DC}) ? $firstGoodEns : $LADCP_atbottom;
my($le) = defined($CTD{TIME_LAG_UC}) ? $realLastGoodEns : $LADCP_atbottom;
for ($ens=$fe; $ens<=$le; $ens++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
my(@residuals) = ();
my(@bindepth) = calc_binDepths($ens);
for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
my($bi) = $bindepth[$bin]/$output_bin_size;
push(@residuals,
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin] -
(($ens < $LADCP_atbottom) ? $DNCAST{MEDIAN_W}[$bi]
: $UPCAST{MEDIAN_W}[$bi]));
}
$LADCP{ENSEMBLE}[$ens]->{MEDIAN_RESIDUAL_W} = median(@residuals);
for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin] -=
$LADCP{ENSEMBLE}[$ens]->{MEDIAN_RESIDUAL_W};
}
}
progress("\tre-calculating profile mad...\n");
for (my($bi)=0; $bi<=$#{$DNCAST{ENSEMBLE}}; $bi++) { # bin data
for (my($i)=0; $i<@{$DNCAST{W}[$bi]}; $i++) {
$DNCAST{W}[$bi][$i] -= $LADCP{ENSEMBLE}[$DNCAST{ENSEMBLE}[$bi][$i]]->{MEDIAN_RESIDUAL_W};
}
$DNCAST{MAD_W}[$bi] = mad2($DNCAST{MEDIAN_W}[$bi],@{$DNCAST{W}[$bi]});
}
for (my($bi)=0; $bi<=$#{$UPCAST{ENSEMBLE}}; $bi++) {
for (my($i)=0; $i<@{$UPCAST{W}[$bi]}; $i++) {
$UPCAST{W}[$bi][$i] -= $LADCP{ENSEMBLE}[$UPCAST{ENSEMBLE}[$bi][$i]]->{MEDIAN_RESIDUAL_W};
}
$UPCAST{MAD_W}[$bi] = mad2($UPCAST{MEDIAN_W}[$bi],@{$UPCAST{W}[$bi]});
}
} # unless ($opt_q);
#-------------------------
# Output Per-Bin Resiudals
#-------------------------
if (defined($out_BR)) {
my(@dc_nres,@uc_nres,@dc_bres,@uc_bres);
progress("Writing per-bin residual data to $out_BR...\n");
for ($ens=$firstGoodEns; $ens<$LADCP_atbottom; $ens++) { # downcast
next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
my(@bindepth) = calc_binDepths($ens);
for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
push(@{$dc_bres[$bin]},
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin]
- $DNCAST{MEDIAN_W}[$bindepth[$bin]/$output_bin_size]);
}
}
for ($ens=$LADCP_atbottom; $ens<=$realLastGoodEns; $ens++) { # upcast
next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
my(@bindepth) = calc_binDepths($ens);
for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
push(@{$uc_bres[$bin]},
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin]
- $UPCAST{MEDIAN_W}[$bindepth[$bin]/$output_bin_size]);
}
}
my($saveParams) = $antsCurParams;
@antsNewLayout = ('bin','dc_residual','dc_residual.sig','dc_residual.nsamp',
,'uc_residual','uc_residual.sig','uc_residual.nsamp');
open(STDOUT,"$out_BR") || croak("$out_BR: $!\n");
&antsAddParams('BR_max_bin',max(scalar(@dc_bres),scalar(@uc_bres)));
for (my($bin)=0; $bin<max(scalar(@dc_bres),scalar(@uc_bres)); $bin++) {
my($dc_avg) = avg(@{$dc_bres[$bin]});
my($uc_avg) = avg(@{$uc_bres[$bin]});
&antsOut($bin+1,$dc_avg,stddev2($dc_avg,@{$dc_bres[$bin]}),scalar(@{$dc_bres[$bin]}),
$uc_avg,stddev2($uc_avg,@{$uc_bres[$bin]}),scalar(@{$uc_bres[$bin]}));
}
&antsOut('EOF'); open(STDOUT,">&2");
$antsCurParams = $saveParams;
}
#--------------------------------------------------
# Calculate BT-referenced vertical-velocity profile
#--------------------------------------------------
if (defined($water_depth)) {
progress("Calculating BT-referenced vertical velocities\n");
calc_BTprof($firstGoodEns,$realLastGoodEns,$water_depth,$sig_water_depth);
my($sumSq) = my($n) = 0;
for (my($bi)=0; $bi<=$#{$BT{MEDIAN_W}}; $bi++) {
next unless defined($BT{MEDIAN_W}[$bi]);
next unless ($BT{N_SAMP}[$bi]>=$min_w_nsamp && $DNCAST{N_SAMP}[$bi]>=$min_w_nsamp && $UPCAST{N_SAMP}[$bi]>=$min_w_nsamp);
$sumSq += ($BT{MEDIAN_W}[$bi] - $DNCAST{MEDIAN_W}[$bi]/2 - $UPCAST{MEDIAN_W}[$bi]/2)**2;
$n++;
}
if ($n > 0) {
my($rms) = round(sqrt($sumSq/$n),0.001);
&antsAddParams('BT_rms_w_discrepancy',$rms);
progress("\t$rms m/s rms vertical-velocity discrepancy\n");
}
}
#------------
# full output
#------------
# NB: residual field is calculated with respect to down-/upcast medians in -o-size bins
if (defined($out_w)) {
progress("Writing vertical-velocity data to $out_w...\n");
@antsNewLayout = ('ensemble','bin','elapsed','depth','CTD_depth','downcast',
'w','w12','w34','residual','CTD_w','LADCP_w','errvel',
'correlation','echo_amplitude','volume_scattering_coeff',
'pitch','roll','tilt','heading','3_beam','svel');
open(STDOUT,"$out_w") || croak("$out_w: $!\n");
if (defined($CTD{TIME_LAG_DC})) {
for ($ens=$firstGoodEns; $ens<$LADCP_atbottom; $ens++) { # downcast
next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
my(@bindepth) = calc_binDepths($ens);
for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
my($bi) = $bindepth[$bin]/$output_bin_size;
&antsOut(
$LADCP{ENSEMBLE}[$ens]->{NUMBER},$bin+1,
$CTD{ELAPSED}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$bindepth[$bin],$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH},1,
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin],
$LADCP{BEAM_COORDINATES} ? $LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin] : nan,
$LADCP{BEAM_COORDINATES} ? $LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin] : nan,
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin] - $DNCAST{MEDIAN_W}[$bi],
$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_W}[$bin],
$LADCP{ENSEMBLE}[$ens]->{ERRVEL}[$bin],
median(@{$LADCP{ENSEMBLE}[$ens]->{CORRELATION}[$bin]}),
median(@{$LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin]}),
$LADCP{ENSEMBLE}[$ens]->{SV}[$bin],
$LADCP{ENSEMBLE}[$ens]->{PITCH},
$LADCP{ENSEMBLE}[$ens]->{ROLL},
$LADCP{ENSEMBLE}[$ens]->{TILT},
$LADCP{ENSEMBLE}[$ens]->{HEADING},
(defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][0]) + # only works for beam coords
defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][1]) +
defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][2]) +
defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][3])) < 4 ? 1 : 0,
$CTD{SVEL}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
);
} # for $bin
} # for $ens
} # if defined dc time lag
if (defined($CTD{TIME_LAG_UC})) {
for ($ens=$LADCP_atbottom; $ens<=$realLastGoodEns; $ens++) { # upcast
next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
my(@bindepth) = calc_binDepths($ens);
for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
my($bi) = $bindepth[$bin]/$output_bin_size;
&antsOut(
$LADCP{ENSEMBLE}[$ens]->{NUMBER},$bin+1,
$CTD{ELAPSED}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$bindepth[$bin],$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH},0,
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin],
$LADCP{BEAM_COORDINATES} ? $LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin] : nan,
$LADCP{BEAM_COORDINATES} ? $LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin] : nan,
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin] - $UPCAST{MEDIAN_W}[$bi],
$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_W}[$bin],
$LADCP{ENSEMBLE}[$ens]->{ERRVEL}[$bin],
median(@{$LADCP{ENSEMBLE}[$ens]->{CORRELATION}[$bin]}),
median(@{$LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin]}),
$LADCP{ENSEMBLE}[$ens]->{SV}[$bin],
$LADCP{ENSEMBLE}[$ens]->{PITCH},
$LADCP{ENSEMBLE}[$ens]->{ROLL},
$LADCP{ENSEMBLE}[$ens]->{TILT},
$LADCP{ENSEMBLE}[$ens]->{HEADING},
(defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][0]) + # only works for beam coords
defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][1]) +
defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][2]) +
defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][3])) < 4 ? 1 : 0,
$CTD{SVEL}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
);
} # for $bin
} # for $ens
} # if defined uc time-lag
&antsOut('EOF'); open(STDOUT,">&2");
}
#---------------
# Output profile
#---------------
if (defined($out_profile)) {
progress("Writing vertical-velocity profiles to <$out_profile>...\n");
@antsNewLayout = ('depth','dc_depth','dc_elapsed','dc_w','dc_w.mad','dc_w.nsamp',
'uc_depth','uc_elapsed','uc_w','uc_w.mad','uc_w.nsamp',
'elapsed','w','w.mad','w.nsamp',
'BT_w','BT_w.mad','BT_w.nsamp');
open(STDOUT,"$out_profile") || croak("$out_profile: $!\n");
undef(%DNCAST) unless defined($CTD{TIME_LAG_DC});
undef(%UPCAST) unless defined($CTD{TIME_LAG_UC});
for (my($bi)=0; $bi<=max($#{$DNCAST{ENSEMBLE}},$#{$UPCAST{ENSEMBLE}},$#{$BT{NSAMP}}); $bi++) {
&antsOut(($bi+0.5)*$output_bin_size, # nominal depth
$DNCAST{MEAN_DEPTH}[$bi],$DNCAST{MEAN_ELAPSED}[$bi],
$DNCAST{N_SAMP}[$bi]>=$min_w_nsamp?$DNCAST{MEDIAN_W}[$bi]:nan,
$DNCAST{MAD_W}[$bi],$DNCAST{N_SAMP}[$bi],
$UPCAST{MEAN_DEPTH}[$bi],$UPCAST{MEAN_ELAPSED}[$bi],
$UPCAST{N_SAMP}[$bi]>=$min_w_nsamp?$UPCAST{MEDIAN_W}[$bi]:nan,
$UPCAST{MAD_W}[$bi],$UPCAST{N_SAMP}[$bi],
$DNCAST{MEAN_ELAPSED}[$bi]/2+$UPCAST{MEAN_ELAPSED}[$bi]/2,
$DNCAST{N_SAMP}[$bi]+$UPCAST{N_SAMP}[$bi]>=$min_w_nsamp ?
($DNCAST{MEDIAN_W}[$bi]*$DNCAST{N_SAMP}[$bi]+$UPCAST{MEDIAN_W}[$bi]*$UPCAST{N_SAMP}[$bi]) / ($DNCAST{N_SAMP}[$bi]+$UPCAST{N_SAMP}[$bi]) :
nan,
$DNCAST{N_SAMP}[$bi]+$UPCAST{N_SAMP}[$bi]>=$min_w_nsamp ?
($DNCAST{MAD_W}[$bi]*$DNCAST{N_SAMP}[$bi]+$UPCAST{MAD_W}[$bi]*$UPCAST{N_SAMP}[$bi]) / ($DNCAST{N_SAMP}[$bi]+$UPCAST{N_SAMP}[$bi]) :
nan,
$DNCAST{N_SAMP}[$bi]+$UPCAST{N_SAMP}[$bi],
$BT{N_SAMP}[$bi]>=$min_w_nsamp?$BT{MEDIAN_W}[$bi]:nan,
$BT{MAD_W}[$bi],$BT{N_SAMP}[$bi]
);
}
&antsOut('EOF'); open(STDOUT,">&2");
}
#--------------------------------------
# write time-series output if requested
#--------------------------------------
if (defined($out_timeseries)) {
progress("Writing time-series data to <$out_timeseries>...\n");
@antsNewLayout = ('ens','elapsed','downcast',
'depth','xducer_sound_speed','pitch','gimbal_pitch','roll','tilt','heading',
'CTD_w','LADCP_reflr_w','LADCP_reflr_w_err','LADCP_reflr_w_nosscorr',
'reflr_ocean_w');
open(STDOUT,"$out_timeseries") || croak("$out_timeseries: $!\n");
for ($ens=$firstGoodEns; $ens<=$realLastGoodEns; $ens++) {
next unless defined($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
my($reflr_oc_w) = defined($LADCP{ENSEMBLE}[$ens]->{REFLR_W})
? $LADCP{ENSEMBLE}[$ens]->{REFLR_W} - $CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]
: undef;
&antsOut($LADCP{ENSEMBLE}[$ens]->{NUMBER},
$CTD{ELAPSED}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
($ens < $LADCP_atbottom) ? 1 : 0,
$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH},
$CTD{SVEL}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$LADCP{ENSEMBLE}[$ens]->{PITCH},
$LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH},
$LADCP{ENSEMBLE}[$ens]->{ROLL},
$LADCP{ENSEMBLE}[$ens]->{TILT},
$LADCP{ENSEMBLE}[$ens]->{HEADING},
$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$LADCP{ENSEMBLE}[$ens]->{REFLR_W},
$LADCP{ENSEMBLE}[$ens]->{REFLR_W_ERR},
$LADCP{ENSEMBLE}[$ens]->{REFLR_W_NOSSCORR},
$reflr_w);
}
&antsOut('EOF'); open(STDOUT,">&2");
}
system("{ sleep 30; ./PostProcess.sh $out_basename $RUN $data_subdir $plot_subdir $log_subdir; }&")
if (-x 'PostProcess.sh');
exit(0);