V2.0
- lots of bug fixes
- major new features:
- dropped CTD scans handled correctly (no more apparent clock drifts)
- support for data files collected with Nortek Signature instruments
- much improved data editing
- significant changes:
- no minimum limit for eps_VKE
- updated for GMT6
- much better data quality information in summary plots
#!/usr/bin/perl
#======================================================================
# L A D C P _ W _ O C E A N
# doc: Fri Dec 17 18:11:13 2010
# dlm: Fri Jul 23 08:50:23 2021
# (c) 2010 A.M. Thurnherr
# uE-Info: 323 36 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# TODO:
# ! use instrument tilt in sidelobe editing
# - use instrument tilt in e.g. calculation of Sv
# - use instrument tilt in PPI editing
#
# - plots:
# - avoid over-plotting axis labels
# - allow for different nsamp magnitudes
# - add "dc" "uc" labels
# - add seabed to eps.VKE profile
# - add seabed to LADCP_w_postproc .wprof output
#
# - worry about water-depth differences (disabled warning)
# - make upcast-flag valid for yoyo casts
# - make diagnostic output 3-beam field work for Earth coordinates
$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
# - HG COMMIT
# Feb 13, 2013: - BUG: CTD_neg_press_offset did not work for CTD time series with -ve starting depth
# Mar 23, 2013: - cosmetics
# - HG COMMIT
# Apr 22, 2013: - removed opt_? variable aliases
# May 8, 2013: - replaced default run label (default to profiles) to make more readable directory structure
# May 14, 2013: - opt_m => w_max_lim
# - BUG: time-series output had two messed-up fields
# May 15, 2013: - added w12 and w34 (2-beam solutions) to profile output
# May 16, 2013: - added CTD_w_tt to time-series & all-sample outputs
# - -a => -d (CTD depth offset)
# - implemented pressure-sensor acceleration correction (-a)
# - added re-gridding of full profile after ping-coherent error removal
# Jun 5, 2013: - renamed $discard_velocities_from_beam to $bad_beam
# - BUG: $bad_beam did not discard BT_VELOCITY data
# Jun 6, 2013: - BUG: error message had -a instead of -d
# Sep 5, 2013: - BUG: w12/w34 do not work for earth-coordinate data, of course
# Apr 17, 2014: - BUG: edit_tilt was never called when all recorded bins are valid
# Apr 21, 2014: - updated comments
# May 19, 2014: - began adding support for PPI filtering
# May 20, 2014: - changed volume_scattering_coeff to Sv in output
# - added editPPI()
# Jul 6, 2014: - BUG: nan water depth had been interpreted as known
# - BUG: sVelProf[] was not allowed to have any gaps
# Jul 9, 2014: - BUG: Jul 6 bug fixes had been applied to older
# version
# - BUG: code meant to ensure gap-free svel profiles did not work correctly
# Jul 12, 2014: - finally made output files executable
# Apr 5, 2015: - added check for required software
# - BUG: removed dc/uc mean w fields from .prof again
# Apr 7, 2015: - made LADCP_w callable from installation directory
# - BUG: -v default was wrong in usage message
# - replaced 'ens' in output files by 'ensemble'
# Apr 16, 2015: - turned output specifies into lists (re-design of
# plotting sub-system)
# - removed 30s sleep from PostProcess.sh call
# - disabled active output when ANTS are not available
# - removed /bin/ksh requirement
# - BUG: error messages were not reported in the log file
# - made seabed detection code more flexible
# - made reference-layer w_ocean warning more severe
# - BUG: info() did not write anything into logfile when format string was used
# - BUG: LADCP_lastBin was set 1 too low
# Apr 20, 2015: - improved comments
# - improved diagnostic messages (warning on missing CTD temp)
# - added support for empirical backscatter correction
# Apr 21, 2015: - improved screen log output
# Apr 22, 2015: - BUG: $realLastGoodEns could be undefined, breaking plotting routines
# Apr 24, 2015: - removed @ARGS
# - added %profile_id
# May 13, 2015: - loosened "insufficient valid velocities" test for short casts
# May 15, 2015: - added $min_valid_vels
# - BUG: LADCP_atbottom could be less than firstGoodEns
# May 18, 2015: - added %LADCP_pulse_length, %dnXX
# May 20, 2015: - added $PROF as newer alternative to $STN
# - replaced "require ProcessingParams" by "do ProcessingParams"
# Jun 15, 2015: - clarified warning message
# Jul 26, 2015: - added %output_grid_dz %output_grid_minsamp for use by [LADCP_w_regrid]
# - began work on support for [libGMT.pl]
# - -v usage message had wrong default
# - added $outGrid_firstBin & $outGrid_lastBin
# Jul 28, 2015: - added GMT plot support for all @out_ lists
# Jul 29, 2015: - continue adaptation of code for new plotting system
# Jul 30, 2015: - cosmetics
# Sep 3, 2015: - changed out_w to out_wsamp
# Sep 26, 2015: - replaced numberp($water_depth) by defined($water_depth) for consistency
# - implemented secondary sidelobe editing
# - allow for -h <filename>
# - allow $ID as alias for $PROF and $STN
# Sep 27, 2015: - BUG: -h filename was broken when water_depth is not known
# Oct 12, 2015: - upgraded missing water-depth warnings to L2
# - added -V)ersion
# - require ANTSlibs V6.2 for release
# Oct 13, 2015: - adapted to [version.pl]
# Nov 25, 2015: - made warning disappear on setting $ANTS_TOOLS_AVAILABLE
# Nov 27, 2015: - changed RDI_BB_READ.pl to RDI_PD0_IO.pl
# Jan 4, 2016: - decreased default vertical resolution to 20m
# Jan 5, 2016: - adapted to [ADCP_tools_lib.pl]
# Jan 22, 2016: - updated for improved mean_residuals plot
# - added per-bin residual quality check
# Jan 25, 2016: - added antsAddParams() with version number
# Jan 26, 2016: - added %processing_params, many others
# - expunged -d
# - implemented outGrid_firstBin eq '*' (also lastBin)
# Jan 27, 2016: - BUG: outGrid_lastBin eq '*' did not work
# - removed large ref-lr w warning
# Feb 14, 2016: - fiddled with ping-coherent residuals code, fixing
# 0-2 potential (minor?) bugs related to outGrid_{first,last}Bin;
# output of new code checked against old code => identical
# Feb 16, 2016: - BUG: per-bin-residual QC could cause division by zero
# Feb 19, 2016: - added -l (disable time-lag filtering)
# Mar 7, 2016: - added error message when -h is neither number nor file
# - BUG: -ve depth error message referred to obsolete -d
# - BUG: dn field name did not use zero filling for year number
# Mar 8, 2016: - removed L0 water-depth-difference warning
# - added test for 1500m/s sound speed
# Mar 9, 2016: - added hab field to .wprof output
# Mar 13, 2016: - cleaned up warnings created before LADCP_file is defined
# - added sanity checks and warnings
# Mar 17, 2016: - added {dc,uc}_rms_{tilt,w_pkg}
# - replaced a couple of **2 by &SQR()
# - replaced %ADCP_orientation values by DL & UL
# - {dc,uc}_rms_w_pkg => {dc,uc}_rms_accel_pkg for V1.2beta6
# Mar 18, 2016: - added -l to %processing_options
# Mar 24, 2016: - generalized plotting system (renamed variables)
# - made plotting system search current dir before LADCP_w dir
# - BUG: ProcessingParams syntax errors were not flagged
# - BUG: progress message in rarely used filter
# Mar 25, 2016: - added -r)esidual rms filter
# - added -q -r to processing_options
# Mar 26, 2016: - BUG: water_depth < CTD_maxdepth was allowed
# - round(CTD_depth)
# Mar 29, 2016: - split [default_paths.pl] from [defaults.pl] to allow defaults
# to be read before usage
# - renamed _subdir to _dir variables
# - renamed post-process hook script to LADCP_w.PostProcess
# - added -r default (0.04m/s)
# Mar 31, 2016: - changed version %PARAM
# Apr 6, 2016: - BUG: GMT test looked for psxy, rather than gmt binary
# May 12, 2106: - added v12, v34 to wsamp output
# May 17, 2016: - added residual12, residual34 to wsamp output
# - changed pitch field in wsamp output to gimbal pitch
# - added w_t to .wsamp output for completeness
# May 18, 2016: - replaced rms_tilt by mean_tilt
# - add (gimbal-)pitch, roll and residual{12,34} to wprof output
# - expunged super-confusing $realLastGoodEns
# - removed assumption of default 1500m/s sspd setting during acquisition
# - added w12, w34 for earth-coordinates
# May 19, 2016: - updated to ADCP_tools V1.6
# May 24, 2016: - calc_binDepths() -> binDepths()
# May 26, 2016: - added -d)isable bin interpolation
# - adapted to w_CTD, w_winch
# - re-enabled disabled L0 warning
# - suppressing L0 warnings on screen
# May 27, 2016: - cosmetics
# - added reflr_w to .wsamp output
# May 29, 2016: - BUG: do ProcessingParams did not handle all errors correctly
# Jun 1, 2016: - made empirical accoustic backscatter correction optional
# - split [default_paths.pl] to same plus [default_output.pl]
# Jun 2, 2016: - added applyTiltCorrection()
# - BUG: PPI-editing was no longer enabled automatically
# for 150kHz instruments
# Jun 3, 2016: - BUG: earth-coordinate data did not have gimbal-pitch
# calculated
# Jun 6, 2016: - removed applyTiltCorrection()
# Jun 11, 2016: - BUG: w12, w34 were wrong for Earth-coord data
# Jul 7, 2016: - BUG: finally fixed with velEarthToBPw() from [RDI_Coords.pl]
# Jul 31, 2016: - BUG: -d did not work because it was handled in [defaults.pl]
# Oct 16, 2016: - cosmetics
# Dec 22, 2016: - moved $opt_p to [defaults.pl]
# Dec 23, 2016: - BUG: -u did not set required variables to proceed
# Mar 6, 2017: - BUG: division by zero when water-depth ~ max(CTD_depth)
# Oct 12, 2017: - BUG: beampair w did not work for earth-coord vels; major re-write
# of earthcoord code to unify processing
# Nov 26, 2017: - BUG: $bad_beam did not work correctly with bin interpolation
# - BUG: ping-coherent residual removal did not respect missing values
# Nov 28, 2017: - added $initial_time_lag
# - expanded semantics of -q to disable time-lagging and residual filters
# Dec 9, 2017: - added $antsSuppressCommonOptions = 1;
# Dec 17, 2017: - added dependencies
# Apr 24, 2018: - added support for $water_depth_db_cmd
# May 1, 2018: - added threshold for reference-layer horizontal speed
# - added ambiguity velocity check
# May 2, 2018: - BUG: ref-lr threshold did not work
# - BUG: BT code was called for UL when -h was used
# - replaced $PPI_seabed_editing_required by &PPI_seabed_editing_required
# - BUG: surface PPI editing code could not be enabled; added &PPI_surface_editing_required
# Nov 2, 2018: - BUG: for 3-beam solutions, residual{12,34} with affected beam was wrong
# Apr 12, 2019: - added logfile error message for non-existing file
# Apr 18, 2019: - added addtional RDI coordinate transformation %PARAMs
# Apr 21, 2019: - BUG: removed duplicate warning about lacking water-depth info
# Sep 10, 2019: - BUG: -i was erroneously listed as the initial CTD time offset
# Sep 12, 2019: - modified .prof output layout to be consistent with [LADCP_w_postproc]
# May 15, 2020: - BUG: re-binning on -r did not work for v12 and v34
# - BUG: before re-binning, upcast v12 and v34 were actuall from DC data
# - BUG: all beam-pair velocity stats were contaminated by missing values
# - made bin residuals per beam pair
# May 16, 2020: - BUG: write_residual_profs() clobbered layout
# Jun 22, 2020: - sscorr_w renamed to remove conflict with [RDI_Coords.pl]
# Mar 17, 2021: - turned multi-ping error into warning (for Nortek files)
# Mar 23, 3021: - adapted to ADCP_tools V2.4 (Nortek PD0 files)
# Jun 30, 2021: - BUG: rms instead of mean tilt was reported in fig
# - added pitch/roll/tilt to output profile
# Jul 1, 2021: - made %PARAMs more standard
# - added %dc_w.mu, %uc_w.mu
# Jul 6, 2021: - BUG: %uc_w.mu was off by factor 2
# Jul 9, 2021: - added LR_RMS_BP_RESIDUAL, dc_lr_bp_res.rms, uc_lr_bp_res.rms
# - added residuals profile filter
# Jul 11, 2021: - added %dc_w.var, %uc_w.var
# - reduced warning associated with residuals profile filter
# - disable time-lag filtering for transmission-error corrected CTD files
# Jul 13, 2021: - set dc_w.mu, dc_w.var to nan on insufficient data
# - BUG: residuals profile filter bombed when no data were removed
# - BUG: rms BT discrepancy was broken (residuals profile filter)
# - BUG: *_w.var did not respect -k
# Jul 23, 2021: - added %ADCP_type
# HISTORY END
# CTD REQUIREMENTS
# - elapsed elapsed seconds; see note below
# - depth
# - ss sound speed
# - w[_CTD] ddepth/dt
# - w_winch OPTIONAL; winch-speed estimate
# - temp OPTIONAL; used for backscatter calculation
# - %lat/%lon OPTIONAL
# - %CTD_transmission_errors OPTIONAL; if defined, -l is implicitly enabled
# OUTPUT PITCH
# - gimbal pitch is used
# OUTPUT EQUIVALENCES:
# 1) list -Q %uc_w.var UL/042.wprof
# list dc_w='$dc_w-%dc_w.mu' UL/042.wprof | avg -Qm dc_w
# 2-BEAM SOLUTIONS
# - both for beam- and earth-coordinate data, two separate two-beam
# solutions (w12 & w34) are calculated:
# - w12 corresponds to ROLL axis (plotted with dashed lines)
# - w34 corresponds to PITCH axis (plotted with dotted lines)
# - the 2-beam solutions are largely independent of each other
# => key for detecting data problems, such as package wakes
# 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
# BINNED-PROFILE OUTPUT NOTES
# - residuals are calculated with respect to down-/upcast medians
# - w12 and w34 2-beam solutions are reported without considering
# -k (min samples)
# - 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 $w_max_lim catches those errors
# quite nicely
# - even when the errors are not filtered with $w_max_lim, they do not
# affect the w profiles, as long as the median bin values are used
($ANTS) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
($ADCP_TOOLS) = (`which mkProfile` =~ m{^(.*)/[^/]*$});
($WCALC) = ($0 =~ m{^(.*)/[^/]*$});
$WCALC = '.' if ($WCALC eq '');
$ANTS_TOOLS_AVAILABLE = (`which list 2>/dev/null` ne '');
die("$0: Generic Mapping Tools (GMT) required but not found (bad \$PATH?)\n")
unless (`which gmt` ne '');
die("$0: ANTSlib required but not found (bad \$PATH?)\n")
unless ($ANTS ne '');
die("$0: ADCP Tools required but not found (bad \$PATH?)\n")
unless ($ADCP_TOOLS ne '');
require "$WCALC/version.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 "$ANTS/ants.pl";
require "$ANTS/libstats.pl";
require "$ADCP_TOOLS/ADCP_tools_lib.pl";
&antsAddParams('LADCP_w_ocean::version',$VERSION);
use IO::Handle;
my($TINY) = 1e-16;
#------
# Usage
#------
require "$WCALC/defaults.pl"; # load default/option parameters
$antsParseHeader = 0;
$antsSuppressCommonOptions = 1;
&antsUsage('3:4a:b:c:de:g:h:i:k:lm:n:o:p:qr:s:t:uv:Vw:x:',0,
"[print software -V)ersion] [-v)erbosity <level[$opt_v]>]",
"[require -4)-beam solutions] [-d)isable bin interpolation] [apply beamvel-m)ask <file> if it exists]",
"[valid LADCP -b)ins <bin,bin[$opt_b]>",
"[-c)orrelation <min[$opt_c counts]>] [-t)ilt <max[$opt_t deg]> [-e)rr-vel <max[$opt_e m/s]>]",
"[max -r)esidual <rms.max[,delta.max][$opt_r m/s]>]",
"[-h water <depth|filename>]",
"[max LADCP time-series -g)ap <length[$opt_g s]>]",
"[-i)nitial LADCP time offset <guestimate> [-u)se as final]]",
"[calculate -n) <lags,lags[$opt_n]>] [lag -w)indow <sz,sz[$opt_w s]>] [lag-p)iece <CTD_elapsed_min|+[,...]>]",
"[require top-3) lags to account for <frac[$opt_3]> of all]",
"[disable time-l)ag filtering]",
"[pressure-sensor -a)cceleration-derivative correction <residual/CTD_w_tt>]",
"[-o)utput bin <resolution[$opt_o m]>] [-k) require <min[$opt_k]> samples]",
"[e-x)ecute <perl-expr>]",
"[-q)uick-and-dirty (no single-ping denoising, residual and time-lagging filters)]",
"<profile-id> [run-label]");
if ($opt_V) {
printf(STDERR "+------------------------------+\n");
printf(STDERR "| LADCP_w Software V%s |\n",$VERSION);
printf(STDERR "| (c) 2015-2021 A.M. Thurnherr |\n");
printf(STDERR "+------------------------------+\n");
exit(0);
}
&antsUsageError() if ($opt_u && !defined($opt_i));
&antsUsageError() unless (@ARGV==1 || @ARGV==2);
&antsCardOpt(\$opt_s,0); # skip # initial ensembles
#$opt_p = '+' unless defined($opt_p); # separate uc/dc time lagging
&antsFloatOpt(\$opt_a,1); # pressure acceleration correction
$ID = $PROF = $STN = &antsCardArg(); # station id ($STN for compatibility)
if (@ARGV) { # run label
$RUN = $ARGV[0];
shift;
} else {
$RUN = 'profiles';
}
#----------------------------------------------------------------
# Handle Processing Parameters
# - paths need to be read first to define ProcessingParams file
# and output directories
# - processing params file is read next to allow setting
# plotting_level
# - default_output.pl is read finally to set any missing output
# variables to default values
#----------------------------------------------------------------
require "$WCALC/default_paths.pl"; # define default input/output paths
my($retval) = do $processing_param_file; # load processing parameters
if (!defined($retval) && $@ ne '') {
croak("$processing_param_file: $@\n");
} elsif (!defined($retval) && $! != 0) {
croak("$processing_param_file: $!\n");
}
require "$WCALC/default_output.pl"; # set default output plots and files
$processing_options = "-k $opt_k -o $opt_o -c $opt_c -t $opt_t -e $opt_e -g $opt_g -3 $opt_3";
$processing_options .= ' -l' if defined($opt_l);
if (defined($opt_q)) { # quick-and-dirty
$processing_options .= ' -q';
$opt_l = 1; # disable time-lagging filter
}
if (defined($opt_i)) { # set initial time lag
$processing_options .= " -i $opt_i";
$initial_time_lag = &antsFloatOpt($opt_i);
}
if (defined($opt_x)) { # eval cmd-line expression to override anything
$processing_options .= " -x $opt_x";
eval($opt_x);
}
if ($opt_4) { # disallow 3-beam solutions
$processing_options .= " -4";
$RDI_Coords::minValidVels = 4;
}
if ($opt_d) { # disable bin mapping
$processing_options .= ' -d';
$RDI_Coords::binMapping = 'none';
}
if (defined($opt_r)) { # residuals filters
$processing_options .= " -r $opt_r";
($residuals_rms_max,$residuals_delta_max) = split(',',$opt_r);
croak("$0: cannot decode -r $opt_r\n") unless ($residuals_rms_max > 0);
}
($LADCP_firstBin,$LADCP_lastBin) = split(',',$opt_b); # select valid bins
croak("$0: cannot decode -b $opt_b\n")
unless (numberp($LADCP_firstBin) &&
($LADCP_lastBin eq '*' || numberp($LADCP_lastBin)));
$processing_options .= " -b $opt_b";
$outGrid_firstBin = $LADCP_firstBin if ($outGrid_firstBin eq '*'); # default bins to use in gridded output
$outGrid_lastBin = $LADCP_lastBin if ($outGrid_lastBin eq '*'); # NB: can still be '*'!!!
@number_of_timelag_windows = split(',',$opt_n);
croak("$0: cannot decode -n $opt_n\n")
unless numberp($number_of_timelag_windows[0]) && numberp($number_of_timelag_windows[1]);
$processing_options .= " -n $opt_n";
@length_of_timelag_windows = split(',',$opt_w);
croak("$0: cannot decode -w $opt_w\n")
unless numberp($length_of_timelag_windows[0]) && numberp($length_of_timelag_windows[1]);
$processing_options .= " -w $opt_w";
croak("$0: \$out_basename undefined\n") # plotting routines use this to label the plots
unless defined($out_basename);
#&antsAddParams('RDI_Coords::binMapping',$RDI_Coords::binMapping); # must be set below since binmapping depends on coords
&antsAddParams('processing_options',$processing_options);
&antsAddParams('out_basename',$out_basename);
&antsAddParams('profile_id',$PROF,'run_label',$RUN);
#----------------------------------------------------------------------
# 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 ($opt_v > 1);
}
sub info(@)
{
print(LOGF "\t"),printf(LOGF @_) if defined($out_log);
printf(STDERR @_) if ($opt_v > 1);
}
sub warning(@)
{
my($lvl,@msg) = @_;
if (defined($out_log)) {
print(LOGF "\nWARNING (L$lvl): ");
printf(LOGF @msg);
print(LOGF "\n");
}
return if ($opt_v == 0) || ($lvl == 0);
if ($opt_v == 1) {
print(STDERR "$LADCP_file: ") if defined($LADCP_file);
print(STDERR "WARNING (L$lvl): ");
printf(STDERR @msg);
} else {
print(STDERR "\n-------------\nWARNING (L$lvl): ");
printf(STDERR @msg);
print(STDERR "-------------\n\n")
}
}
sub error($)
{
print(LOGF "ABORT: @_") if defined($out_log);
croak("ABORT: @_");
}
sub debugmsg(@)
{ printf(STDERR @_) if ($opt_v > 2); }
#---------------------------
# Read LADCP data
#
# sanity checks:
# - single-ping ensembles
#
# warnings:
# - high bandwidth
# - low transmit power
#---------------------------
progress("Reading LADCP data from <$LADCP_file>...\n");
error("cannot open <$LADCP_file> for reading\n")
unless (-r $LADCP_file);
readData($LADCP_file,\%LADCP);
&antsAddDeps($LADCP_file);
if ($LADCP{BEAM_COORDINATES}) {
progress("\t%d ensembles (beam coordinates)\n",scalar(@{$LADCP{ENSEMBLE}}));
} else {
progress("\t%d ensembles (Earth coordinates)\n",scalar(@{$LADCP{ENSEMBLE}}));
}
if ($valid_ensemble_range[0] > 0) { # remove leading invalid records
my($ens) = 0;
while ($ens < @{$LADCP{ENSEMBLE}} && $LADCP{ENSEMBLE}[$ens]->{NUMBER}<$valid_ensemble_range[0]) { $ens++ }
splice(@{$LADCP{ENSEMBLE}},0,$ens);
progress("\t%d invalid leading ensembles removed\n",$ens)
}
if ($valid_ensemble_range[1] > 0) { # remove trailing invalid records
my($ens) = 0;
while ($ens < @{$LADCP{ENSEMBLE}} && $LADCP{ENSEMBLE}[$ens]->{NUMBER}<$valid_ensemble_range[1]) { $ens++ }
splice(@{$LADCP{ENSEMBLE}},$ens);
progress("\t%d invalid trailing ensembles removed\n",@{$LADCP{ENSEMBLE}}-$ens)
}
warning(2,"$LADCP_file: multi-ping ensembles\n")
unless ($LADCP{PINGS_PER_ENSEMBLE} == 1);
warning(2,"$LADCP_file: wide-bandwidth setting\n")
if ($LADCP{WIDE_BANDWIDTH});
warning(2,"$LADCP_file: low transmit-power setting\n")
unless ($LADCP{TRANSMIT_POWER_MAX});
error("$LADCP_file: not enough LADCP bins ($LADCP{N_BINS}) for choice of -r\n")
unless ($LADCP{N_BINS} >= $refLr_lastBin);
error("$0: first reference-layer bin outside valid range\n")
unless ($refLr_firstBin>=1 && $refLr_firstBin<=$LADCP{N_BINS});
error("$0: last reference-layer bin outside valid range\n")
unless ($refLr_lastBin>=1 && $refLr_lastBin<=$LADCP{N_BINS});
error("$0: first reference-layer bin > last reference-layer bin\n")
unless ($refLr_firstBin <= $refLr_lastBin);
&antsAddParams('refLr_firstBin',$refLr_firstBin,'refLr_lastBin',$refLr_lastBin);
$LADCP_lastBin = $LADCP{N_BINS} if ($LADCP_lastBin eq '*');
$outGrid_lastBin = $LADCP{N_BINS} if ($outGrid_lastBin eq '*');
error("$0: first valid LADCP bin outside valid range\n")
unless ($LADCP_firstBin>=1 && $LADCP_firstBin<=$LADCP{N_BINS});
error("$0: last valid LADCP bin outside valid range\n")
unless ($LADCP_lastBin>=1 && $LADCP_lastBin<=$LADCP{N_BINS});
error("$0: first valid LADCP bin > last valid LADCP bin\n")
unless ($LADCP_firstBin <= $LADCP_lastBin);
&antsAddParams('LADCP_firstBin',$LADCP_firstBin,'LADCP_lastBin',$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_pulse_length',$LADCP{TRANSMITTED_PULSE_LENGTH},
'ADCP_type',$LADCP{INSTRUMENT_TYPE},
'ADCP_frequency',$LADCP{BEAM_FREQUENCY},
'ADCP_blanking_distance',$LADCP{BLANKING_DISTANCE});
#------------------------------------------------------------
# Edit velocity data
# beam coords:
# 1) beam-vel mask on -m
# mask file has three columns: from_ens to_ens ignore_beam
# 2) correlation threshold
# Earth coords:
# 1) correlation threshold
#------------------------------------------------------------
if ($LADCP{BEAM_COORDINATES}) {
progress("Editing beam-velocity data...\n");
if (defined($opt_m) && -r $opt_m) {
progress("\tmasking beam velocities according to $opt_m:");
my($nee) = 0;
open(BVM,$opt_m) || die("$opt_m: $!\n");
while (<BVM>) {
s/#.*//;
s/^\s*$//;
next if ($_ eq '');
my($fe,$te,$db) = split;
die("$opt_m: cannot decode $_\n")
unless (numberp($fe) && numberp($te) && $te>=$fe && $db>=1 && $db<=4);
die("$0: assertion failed")
unless ($LADCP{ENSEMBLE}[$fe-1]->{NUMBER} == $fe &&
$LADCP{ENSEMBLE}[$te-1]->{NUMBER} == $te);
for (my($ens)=$fe-1; $ens<=$te-1; $ens++) {
$nee++;
for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
undef($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$db-1]);
}
}
}
progress(" $nee ensembles edited\n");
close(BVM);
}
$nvv = $cte = 0;
if ($pitch_bias || $roll_bias || $heading_bias) {
&antsAddParams('pitch_bias',$pitch_bias,
'roll_bias',$roll_bias,
'heading_bias',$heading_bias);
for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
correctAttitude($ens,$pitch_bias,$roll_bias,$heading_bias);
}
}
for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
$nvv += countValidBeamVels($ens);
$cte += editCorr($ens,$opt_c);
}
error("$LADCP_file: no valid data\n") unless ($nvv > 0);
progress("\tcorrelation threshold (-c %d counts): %d velocites removed (%d%% of total)\n",$opt_c,$cte,round(100*$cte/$nvv));
} else { # if BEAM_COORDINATES
progress("Editing velocity data...\n");
error("$LADCP_file: cannot apply beamvel-mask $opt_m to earth-coordinate data\n")
if defined($opt_m);
$nvv = $cte = 0;
for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
$nvv += countValidBeamVels($ens);
$cte += editCorr_Earthcoords($ens,$opt_c);
}
error("$LADCP_file: no valid data\n") unless ($nvv > 0);
progress("\tcorrelation threshold (-c %d counts): %d velocites removed (%d%% of total)\n",$opt_c,$cte,round(100*$cte/$nvv));
}
#------------------------------------------------------------
# Create beam coordinate velocities for Earth-velocity data
# - velocities are replaced "in place"
# - BT velocities are treated separately in [find_seabed.pl]
# - this transformation will remove all 3-beam solutions
# - disable bin mapping because Earth coords are typically bin-remapped
#------------------------------------------------------------
unless ($LADCP{BEAM_COORDINATES}) {
progress("Replacing earth- with beam-velocities...\n");
for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
@{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]} =
&velEarthToBeam(\%LADCP,$ens,@{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]});
}
}
$RDI_Coords::binMapping = 'none';
}
&antsAddParams('RDI_Coords::binMapping', $RDI_Coords::binMapping, # finally, bin mapping is known
'RDI_Coords::minValidVels', $RDI_Coords::minValidVels,
'RDI_Coords::beamTransformation',$RDI_Coords::beamTransformation);
#-------------------------------------------------------------------
# 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)
#-------------------------------------------------------------------
my($dummy);
progress("Calculating earth-coordinate velocities...\n");
if ($bad_beam) {
progress("\tdiscarding velocities from beam $bad_beam\n");
&antsAddParams('bad_beam_discarded',$bad_beam);
}
$nvw = 0;
for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
$LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH} =
gimbal_pitch($LADCP{ENSEMBLE}[$ens]->{PITCH},$LADCP{ENSEMBLE}[$ens]->{ROLL});
if ($bad_beam) {
for (my($bin)=0; $bin<=$LADCP{N_BINS}-1; $bin++) {
undef($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$bad_beam-1]);
undef($LADCP{ENSEMBLE}[$ens]->{BT_VELOCITY}[$bin][$bad_beam-1]);
}
}
for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
($LADCP{ENSEMBLE}[$ens]->{INTERP_U}[$bin],
$LADCP{ENSEMBLE}[$ens]->{INTERP_V}[$bin],
$LADCP{ENSEMBLE}[$ens]->{INTERP_W}[$bin],
$LADCP{ENSEMBLE}[$ens]->{INTERP_ERRVEL}[$bin]) = earthVels(\%LADCP,$ens,$bin);
($LADCP{ENSEMBLE}[$ens]->{V12}[$bin],$LADCP{ENSEMBLE}[$ens]->{W12}[$bin],
$LADCP{ENSEMBLE}[$ens]->{V34}[$bin],$LADCP{ENSEMBLE}[$ens]->{W34}[$bin]) = BPEarthVels(\%LADCP,$ens,$bin);
}
for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
$LADCP{ENSEMBLE}[$ens]->{U}[$bin] = $LADCP{ENSEMBLE}[$ens]->{INTERP_U}[$bin];
$LADCP{ENSEMBLE}[$ens]->{V}[$bin] = $LADCP{ENSEMBLE}[$ens]->{INTERP_V}[$bin];
$LADCP{ENSEMBLE}[$ens]->{W}[$bin] = $LADCP{ENSEMBLE}[$ens]->{INTERP_W}[$bin];
$LADCP{ENSEMBLE}[$ens]->{ERRVEL}[$bin] = $LADCP{ENSEMBLE}[$ens]->{INTERP_ERRVEL}[$bin];
undef($LADCP{ENSEMBLE}[$ens]->{INTERP_U}[$bin]);
undef($LADCP{ENSEMBLE}[$ens]->{INTERP_V}[$bin]);
undef($LADCP{ENSEMBLE}[$ens]->{INTERP_W}[$bin]);
undef($LADCP{ENSEMBLE}[$ens]->{INTERP_ERRVEL}[$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");
progress("\t3-beam solutions : $RDI_Coords::threeBeam_1 " .
"$RDI_Coords::threeBeam_2 " .
"$RDI_Coords::threeBeam_3 " .
"$RDI_Coords::threeBeam_4\n")
unless ($opt_4);
error("$LADCP_file: insufficient valid velocities\n") unless ($nvw >= $min_valid_vels);
#----------------------------------------------
# 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");
&antsAddParams('per_ens_outliers_mad_limit',$per_ens_outliers_mad_limit)
&antsAddParams('farthest_valid_bins_truncated',$truncate_farthest_valid_bins)
if ($truncate_farthest_valid_bins);
$evrm = $trrm = $worm = 0;
for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
$evrm += editErrVel($ens,$opt_e);
$worm += editWOutliers($ens,$per_ens_outliers_mad_limit);
$trrm += editTruncRange($ens,$truncate_farthest_valid_bins) if ($truncate_farthest_valid_bins);
}
progress("\terror-velocity threshold (-e %.2f m/s): %d velocites removed (%d%% of total in bins $LADCP_firstBin-$LADCP_lastBin)\n",
$opt_e,$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,$opt_g);
error("$LADCP_file: insufficient valid data\n")
unless defined($firstGoodEns) && ($lastGoodEns>$firstGoodEns) && ($LADCP_atbottom>=$firstGoodEns);
my($cast_duration) = $LADCP{ENSEMBLE}[$lastGoodEns]->{ELAPSED} -
$LADCP{ENSEMBLE}[$firstGoodEns]->{ELAPSED};
error("$0: implausibly short cast ($cast_duration seconds)\n")
unless ($cast_duration > 600);
my($year) = $LADCP{ENSEMBLE}[$firstGoodEns]->{YEAR} % 100;
&antsAddParams(sprintf('dn%02d',$year),$LADCP{ENSEMBLE}[$firstGoodEns]->{DAYNO});
$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 (@out_LADCPtis) {
progress("Writing LADCP-timeseries to ");
@antsNewLayout = ('ensemble','elapsed','reflr_w','reflr_w.stddev','reflr_w.nsamp','depth');
foreach my $of (@out_LADCP) {
progress("<$of> ");
my($sub,$arg) = ($of =~ /^([^\(]+)\(([^\)]+)\)$/); # sub(arg), e.g. plot_wprof(DL/003_wprof.ps)
if (defined($arg)) { # NB: exactly one arg
if (-f "./${sub}.pl") { # # NB: don't quote arg
require "./${sub}.pl";
} else {
require "$WCALC/${sub}.pl";
}
&{$sub}($arg);
next;
}
$of = ">$of" unless ($of =~ /^$|^\s*\|/);
open(STDOUT,$of) || error("$of: $!\n");
undef($antsActiveHeader) unless ($ANTS_TOOLS_AVAILABLE);
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');
}
progress("\n");
}
error("deepest depth is at end of cast (no upcast data)\n")
if ($lastGoodEns-$LADCP_atbottom < 100);
#----------------------------------------------------------------------
# More editing
# - attitude threshold
# - data in far bins (beyond reliable range)
# - at this stage ${first,last}GoodEns are known
# - TILT field is set as a side-effect
#----------------------------------------------------------------------
progress("Editing additional earth-coordinate velocity data...\n");
&antsAddParams('per_bin_valid_frac_lim',$per_bin_valid_frac_lim);
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;
}
my($dc_sumTilts,$dc_n,$uc_sumTilts,$uc_n) = (0,0,0,0);
$fprm = $pte = 0;
for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
my($ne) = editTilt($ens,$opt_t); # this sets the TILT field
if ($ne == 0) {
if ($ens <= $LADCP_atbottom) {
$dc_sumTilts += $LADCP{ENSEMBLE}[$ens]->{TILT};
$dc_n++;
} else {
$uc_sumTilts += $LADCP{ENSEMBLE}[$ens]->{TILT};
$uc_n++;
}
} else {
$pte += $ne;
}
$fprm += editFarBins($ens,$first_bad_bin) if defined($first_bad_bin);
}
my($dc_mean_tilt) = ($dc_n > 0) ? $dc_sumTilts/$dc_n : nan;
my($uc_mean_tilt) = ($uc_n > 0) ? $uc_sumTilts/$uc_n : nan;
&antsAddParams('dc_tilt.mu',$dc_mean_tilt,'uc_tilt.mu',$uc_mean_tilt);
progress("\tattitude threshold (max_tilt = %d deg): %d velocites removed (%d%% of total)\n",
$opt_t,$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) || error("$CTD_file: $!\n");
error("$CTD_file: no data\n") unless (&antsIn());
undef($antsOldHeaders);
$opt_l = 1 # disable time-lag filtering
if defined($P{CTD_transmission_errors}); # for transmission-corrected CTD files
&antsAddDeps($CTD_file);
&antsAddParams('lat',$P{lat}) if defined($P{lat});
&antsAddParams('lon',$P{lon}) if defined($P{lon});
($CTD_elapsed,$CTD_depth,$CTD_svel) = &fnr('elapsed','depth','ss');
$CTD_w = &fnrNoErr('w_CTD');
$CTD_w = &fnr('w') unless defined($CTD_w);
$CTD_winchvel = &fnrNoErr('w_winch');
$CTD_temp = &fnrNoErr('temp');
warning(0,"no CTD temperature --- using ADCP temperature instead => Sv degraded!\n",$s)
unless defined($CTD_temp);
$CTD_maxdepth = -1;
do { # read data
error("$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);
push(@{$CTD{W_WINCH}},$ants_[0][$CTD_winchvel]) if defined($CTD_winchvel);
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});
$CTD{W_t}[0] = $CTD{W_t} [@{$CTD{ELAPSED}}] = nan; # calculate w-derivatives
$CTD{W_tt}[0] = $CTD{W_tt}[@{$CTD{ELAPSED}}] = nan;
for (my($s)=1; $s<@{$CTD{ELAPSED}}-1; $s++) { # centered differences for both
$CTD{W_t} [$s] = ($CTD{W}[$s+1] - $CTD{W}[$s-1]) / (2*$CTD{DT});
$CTD{W_tt}[$s] = ($CTD{W}[$s+1] + $CTD{W}[$s-1] - 2*$CTD{W}[$s]) / &SQR($CTD{DT});
}
$CTD_maxdepth = round($CTD_maxdepth);
error("$0: CTD start depth must be numeric\n")
unless numberp($CTD{DEPTH}[0]);
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
# For PPI filtering, a sound speed profile without gaps is required.
#--------------------------------------------------------------------
progress("Constructing sound-speed correction profile\n");
my($scans_per_sec) = int(1/$CTD{DT}+0.5);
my($min_depth) = 9e99;
for (my($s)=0; $s<=$CTD_atbottom; $s+=$scans_per_sec) {
next unless ($CTD{DEPTH}[$s] >= 0 && numberp($CTD{SVEL}[$s]));
my($d) = int($CTD{DEPTH}[$s]);
$min_depth = $d if ($d < $min_depth);
$sVelProf[$d] = $CTD{SVEL}[$s];
}
while ($min_depth > 0) { # fill surface gap
$sVelProf[$min_depth-1] = $sVelProf[$min_depth];
$min_depth--;
}
for (my($d)=$min_depth+1; $d<=$#sVelProf; $d++) { # fill interior gaps
$sVelProf[$d] = $sVelProf[$d-1]
unless defined($sVelProf[$d]);
}
#-------------------
# Determine time lag
#-------------------
if (defined($initial_time_lag)) {
progress("Setting initial time lag...\n");
$CTD{TIME_LAG} = $initial_time_lag;
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");
$CTD_time_lag[0] = $CTD{TIME_LAG};
$CTD_tl_toEns[0] = $lastGoodEns;
} else {
#------------------------
# stage 1: 1Hz; full cast
#------------------------
$CTD{TIME_LAG} =
calc_lag($number_of_timelag_windows[0],$length_of_timelag_windows[0],
int(1/$CTD{DT}+0.5),$firstGoodEns,$lastGoodEns);
error("$0: Cannot proceed without valid lag!\n") unless defined($CTD{TIME_LAG});
progress("\telapsed(CTD) ~ elapsed(LADCP) + %.2fs\n",$CTD{TIME_LAG});
#---------------------------------
# stage 2: piece-wise time lagging
#---------------------------------
my(@splits);
push(@splits,$firstGoodEns);
foreach my $s (split(/,/,$opt_p)) {
if ($s eq '+') {
push(@splits,$LADCP_atbottom);
} elsif (numberp($s)) { # find ens for given LADCP elapsed time in min
my($ela) = $s*60;
my($ens) = $firstGoodEns + round($ela / $LADCP{MEAN_DT}); # first, guess
while ($LADCP{ENSEMBLE}[$ens]->{ELAPSED} > $ela) { $ens--; } # then, correct
while ($LADCP{ENSEMBLE}[$ens]->{ELAPSED} < $ela) { $ens++; }
push(@splits,$ens);
} else {
warning(2,"ignoring non-numeric -p argument %s",$s);
}
}
push(@splits,$lastGoodEns);
my($valid_lag); # determine piece-wise lags
while (@splits > 1) {
push(@CTD_tl_fromEns,$splits[0]);
push(@CTD_tl_toEns,$splits[1]);
debugmsg("lag($splits[0],$splits[1]) = ");
my($lag) = calc_lag($number_of_timelag_windows[1],$length_of_timelag_windows[1],
1,$splits[0],$splits[1]);
debugmsg("$lag\n");
if (defined($lag)) {
progress("\tcast-piece: elapsed(CTD) = elapsed(LADCP) + %.2fs\n",$lag);
push(@CTD_time_lag,$lag);
$valid_lag = $#CTD_time_lag;
} else {
progress("\tcast-piece: time-lagging failed\n");
push(@CTD_time_lag,nan);
}
shift(@splits);
}
error("$0: Cannot proceed without at least one lag!\n") # fill failed lag with surrounding data
unless defined($valid_lag);
while ($valid_lag < $#CTD_time_lag) { # forward
$CTD_time_lag[$valid_lag+1] = $CTD_time_lag[$valid_lag];
$valid_lag++;
}
while ($valid_lag > 0) { # backward
$CTD_time_lag[$valid_lag-1] = $CTD_time_lag[$valid_lag]
unless (numberp($CTD_time_lag[$valid_lag-1]));
$valid_lag--;
}
}
&antsAddParams('CTD_time_lags',"@CTD_time_lag");
undef($CTD{TIME_LAG}); # to avoid confusion, clear old 1Hz-based lag
#------------------------------------------------
# Merge (associate) CTD with LADCP data
# - after this step, reflr w is sound-speed corrected!!!
#------------------------------------------------
progress("Merging CTD with LADCP data...\n");
&antsAddParams('w_lim.max',$w_max_lim);
my($cli) = 0; # current-lag index
my($lag) = $CTD_time_lag[$cli]; # current lag
my($dc_sumsq,$dc_n,$uc_sumsq,$uc_n) = (0,0,0,0); # for a_pkg calcs
for (my($skipped)=0,my($ens)=$firstGoodEns; $ens<=$lastGoodEns; $ens++) {
if ($ens > $CTD_tl_toEns[$cli]) { # use correct lag piece
$cli++;
die("assertion failed!\n\ttest: \$cli(=$cli) <= \$#CTD_time_lag(=$#CTD_time_lag) at ens=$ens\n")
unless ($cli <= $#CTD_time_lag);
$lag = $CTD_time_lag[$cli];
}
my($scan) =
int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED} + $lag - $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("\t$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("\t%d final LADCP ensembles skipped because CTD data end with LADCP in water\n",
$lastGoodEns-$ens+1));
$lastGoodEns = $ens-1;
last;
}
$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} *= # correct for sound-speed variations at source
$CTD{SVEL}[$scan] / $LADCP{ENSEMBLE}[$ens]->{SPEED_OF_SOUND};
error(sprintf("\n$0: negative depth (%.1fm) in CTD file at elapsed(CTD) = %.1fs\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; # used for other CTD fields
if (numberp($CTD{W_t}[$scan])) {
if ($ens <= $LADCP_atbottom) {
$dc_sumsq += &SQR($CTD{W_t}[$scan]); $dc_n++;
} else {
$uc_sumsq += &SQR($CTD{W_t}[$scan]); $uc_n++;
}
}
my($reflr_ocean_w) = $LADCP{ENSEMBLE}[$ens]->{REFLR_W} - $CTD{W}[$scan];
if (abs($reflr_ocean_w) <= $w_max_lim) {
$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
}
}
my($dc_rms_apkg) = ($dc_n > 0) ? sqrt($dc_sumsq / $dc_n) : nan; # rms package acceleration
my($uc_rms_apkg) = ($uc_n > 0) ? sqrt($uc_sumsq / $uc_n) : nan;
&antsAddParams('dc_accel_pkg.rms',$dc_rms_apkg,'uc_accel_pkg.rms',$uc_rms_apkg);
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));
if (sqrt($sumWsqI/$nWsqI) > 0.05) {
warning(0,"%.2f cm/s (large) reference-layer w_ocean away from boundaries\n",100*sqrt($sumWsqI/$nWsqI));
} elsif (sqrt($sumWsqI/$nWsqI) > 0.15) {
warning(2,"%.2f cm/s (implausibly large) reference-layer w_ocean away from boundaries\n",100*sqrt($sumWsqI/$nWsqI));
}
} 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 {
error("$0: no valid vertical velocities\n");
}
#----------------------------------------------------------------------
# Calculate Volume Scattering Coefficients
#----------------------------------------------------------------------
progress("Calculating volume-scattering coefficients (Sv)...\n");
calc_backscatter_profs($firstGoodEns,$lastGoodEns);
if (defined($Sv_ref_bin) && @nSv) {
progress("\tapplying empirical Sv correction\n");
&antsAddParams('Sv_correction','empirical');
correct_backscatter($firstGoodEns,$lastGoodEns);
}
#----------------------------------------------------------------------------
# Find Seabed & Edit data
# 1) contaminated by sidelobe reflection from seabed and sea surface
# 2) PPI
#----------------------------------------------------------------------------
error("$0: conflicting water-depth information provided by user\n") # only happens when $water_depth is set explicitly
if defined($opt_h) && defined($water_depth);
if (defined($opt_h)) { # handle user-provided water-depth info
if (numberp($opt_h)) {
$water_depth = $opt_h;
} elsif (-f $opt_h) {
open(WDF,$opt_h) || error("$opt_h: $_\n");
$water_depth = &antsFileScanParam(WDF,'water_depth');
close(WDF);
undef($water_depth) unless numberp($water_depth);
} else {
error("$0: -h $opt_h defines neither number nor existing file\n");
}
}
if (!defined($water_depth) && # find seabed in data
$LADCP{ENSEMBLE}[$LADCP_atbottom]->{XDUCER_FACING_DOWN}) {
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) && defined($water_depth_BT)) {
my($dd) = abs($water_depth_BT - $water_depth);
# warning(0,sprintf("Large instrument vs. backscatter-derived water-depth difference (%.1fm)\n",$dd))
# if ($dd > 10);
}
if (!$SS_use_BT && !defined($water_depth) && defined($water_depth_BT)) { # warn if use of BT was not
warning(1,"using water_depth from ADCP BT data\n"); # explicitly requested
$SS_use_BT = 1;
}
if ($SS_use_BT && numberp($water_depth_BT)) { # water depth from BT data
&antsAddParams('water_depth_from','BT_data');
$water_depth = $water_depth_BT;
$sig_water_depth = $sig_water_depth_BT;
} elsif (defined($water_depth)) { # water depth from WT data
&antsAddParams('water_depth_from','echo_amplitudes');
}
}
if (defined($water_depth)) { # 1 or 2 water depths available
if (defined($water_depth_BT)) {
progress("\t%.1f(%.1f) m water depth (%.1f(%.1f)m from ADCP BT data)\n",
$water_depth,$sig_water_depth,$water_depth_BT,$sig_water_depth_BT);
} else {
progress("\t%.1f(%.1f) m water depth (no seabed found in BT data)\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});
if ($water_depth < $CTD_maxdepth) {
warning(2,"water depth ($water_depth m) < max CTD depth ($CTD_maxdepth m) ignored\n");
undef($water_depth);
}
}
if (!defined($water_depth) && defined($water_depth_db_cmd)) { # set water depth from data base
error("$0: lat/lon required for running $water_depth_db_cmd\n")
unless numbersp($P{lat},$P{lon});
chomp($water_depth = `$water_depth_db_cmd $P{lon} $P{lat}`);
error("$0: command '$water_depth_db_cmd $P{lon} $P{lat}' did not return valid water depth\n")
unless numberp($water_depth);
&antsAddParams('water_depth_from',"$water_depth_db_cmd $P{lon} $P{lat}");
}
if (defined($water_depth)) { # set %PARAMs
&antsAddParams('water_depth',$water_depth,'water_depth.sig',$sig_water_depth);
} else {
warning(2,"unknown water depth --- cannot edit sidelobes or PPI near the seabed\n");
&antsAddParams('water_depth','unknown','water_depth.sig','nan');
}
if ($LADCP{ENSEMBLE}[$LADCP_atbottom]->{XDUCER_FACING_DOWN}) { # DOWNLOOKER
&antsAddParams('ADCP_orientation','DL');
if (defined($water_depth)) {
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");
if ($sidelobe_editing_DL_surface) {
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");
&antsAddParams('sidelobe_editing','surface+seabed','vessel_draft',$vessel_draft);
} else {
&antsAddParams('sidelobe_editing','seabed');
}
# Setting $DEBUG_DL_SURFACE_PPI_EDITING to 1 enables sea-surface PPI editing for
# downlookers (but disabled seabed editing). This was used on 2021 A20 for testing.
if ($DEBUG_DL_SURFACE_PPI_EDITING) {
&antsAddParams('PPI_editing','surface');
&antsAddParams('PPI_extend_upper_limit',$PPI_extend_upper_limit)
if numberp($PPI_extend_upper_limit);
progress("Editing data to remove PPI from sea surface...\n");
progress("\tConstructing depth-average soundspeed profile...\n");
$DASSprof[0] = my($sum) = 0;
for (my($d)=1; $d<=$#sVelProf; $d++) {
die("assertion failed") unless numberp($sVelProf[$d]);
$sum += $sVelProf[$d];
$DASSprof[$d] = $sum/$d;
}
($nvrm,$nerm) = editPPI($firstGoodEns,$lastGoodEns,undef);
progress("\t$nvrm velocities from $nerm ensembles removed\n");
} else {
if (&PPI_seabed_editing_required()) {
&antsAddParams('PPI_editing','seabed');
&antsAddParams('PPI_extend_upper_limit',$PPI_extend_upper_limit)
if numberp($PPI_extend_upper_limit);
progress("Editing data to remove PPI from seabed...\n");
progress("\tConstructing depth-average soundspeed profile...\n");
die("assertion failed") unless defined($water_depth);
my($dz) = $water_depth - $#sVelProf; # $#sVelProf = max_depth(profile) in meters
my($sum) = $dz * $sVelProf[$#sVelProf];
if ($dz == 0) { # water-depth <= max CTD depth
warning(1,"inferred water depth very close to max(CTD depth)\n");
$DASSprof[$#sVelProf] = $sVelProf[$#sVelProf];
} else {
die("assertion failed") unless defined($dz > 0);
$DASSprof[$#sVelProf] = $sum/$dz;
}
for (my($d)=$#sVelProf-1; $d>=0; $d--) {
die("assertion failed (d=$d, #sVelProf=$#sVelProf)") unless numberp($sVelProf[$d]);
$sum += $sVelProf[$d];
$dz++;
$DASSprof[$d] = $sum/$dz;
}
($nvrm,$nerm) = editPPI($firstGoodEns,$lastGoodEns,$water_depth);
progress("\t$nvrm velocities from $nerm ensembles removed\n");
}
}
}
} else { # UPLOOKER
&antsAddParams('ADCP_orientation','UL');
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");
if ($sidelobe_editing_UL_seabed) {
if (defined($water_depth)) {
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 {
# WARNING ALREADY PRODUCED ABOVE
# warning(2,"unknown water depth --- cannot edit UL data for sidelobe interference from seabed\n");
}
&antsAddParams('sidelobe_editing','surface+seabed','vessel_draft',$vessel_draft);
} else {
&antsAddParams('sidelobe_editing','surface','vessel_draft',$vessel_draft);
}
if (&PPI_surface_editing_required()) {
&antsAddParams('PPI_editing','surface');
&antsAddParams('PPI_extend_upper_limit',$PPI_extend_upper_limit)
if numberp($PPI_extend_upper_limit);
progress("Editing data to remove PPI from sea surface...\n");
progress("\tConstructing depth-average soundspeed profile...\n");
$DASSprof[0] = my($sum) = 0;
for (my($d)=1; $d<=$#sVelProf; $d++) {
die("assertion failed") unless numberp($sVelProf[$d]);
$sum += $sVelProf[$d];
$DASSprof[$d] = $sum/$d;
}
($nvrm,$nerm) = editPPI($firstGoodEns,$lastGoodEns,undef);
progress("\t$nvrm velocities from $nerm ensembles removed\n");
}
}
#----------------------------------------------------------------------
# Check For Ambiguity Velocity Problems
#----------------------------------------------------------------------
progress("Checking for ambiguity velocity violations...\n");
my($ambiguity_velocity) = ambiguity_velocity($LADCP{BEAM_FREQUENCY},
$LADCP{BEAM_ANGLE},
$LADCP{SPEED_OF_SOUND},
$LADCP{TRANSMIT_LAG_DISTANCE});
&antsAddParams('ambiguity_velocity',$ambiguity_velocity);
my($nbad) = 0;
for (my($ens)=$firstGoodEns; $ens<=$lastGoodEns; $ens++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
next unless ($CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}] > $ambiguity_velocity);
$nbad++;
}
my($badf) = $nbad / ($lastGoodEns - $firstGoodEns + 1); # fraction of bad values
if ($bad > 0.01) { # allow 1% violations before warning is triggered
warning(2,"%d ensembles (%d%% of total) have CTD_w > ambiguity velocity of %.1 m/s\n",
$nbad,round(100*$badf),$ambiguity_velocity);
} elsif ($nbad > 0) {
info("\t%d ensembles (%d%% of total) have CTD_w > ambiguity velocity of %.1 m/s",
$nbad,round(100*$badf),$ambiguity_velocity);
} else {
info("\tnone found\n");
}
#----------------------------------------------------------------------
# Data Editing after LADCP and CTD data have been merged
# 1) surface layer editing
# 2) reference-layer horizontal velocity threshold
# 3) Execute user-supplied $edit_data_hook
#----------------------------------------------------------------------
progress("Removing data from instrument at surface...\n");
&antsAddParams('surface_layer_depth',$surface_layer_depth);
$nerm = editSurfLayer($firstGoodEns,$lastGoodEns,$surface_layer_depth);
progress("\t$nerm ensembles removed\n");
progress("Removing data collected at large horizontal package speeds...\n");
$max_hspeed = &max_hspeed(); # defined in [defaults.h]
&antsAddParams('max_hspeed',$max_hspeed);
$nerm = editLargeHSpeeds($firstGoodEns,$lastGoodEns,$max_hspeed);
my($nermf) = $nerm / ($lastGoodEns - $firstGoodEns + 1);
info("\treference-layer horizontal speed threshold (max_hspeed = %g m/s): %d ensembles removed (%d%% of total)\n",
$max_hspeed,$nerm,round(100*$nermf));
warning(2,"large fraction (%d%%) of samples exceed reference-layer horizontal speed threshold\n",round(100*$nermf))
if ($nermf > 0.05);
if (defined($post_merge_hook)) {
progress("Executing user-supplied \$post_merge_hook...\n");
&{$post_merge_hook}($firstGoodEns,$lastGoodEns);
}
#----------------------------------------------------------------------
# apply pressure-sensor acceleration correction (-a)
# - valid, but useless except for troubleshooting data sets (e.g. DIMES UK2)
# - useless, because ping-coherent error removal below does same job but better
# - it does not make sense to to use this with ping-coherent error removal enabled
#----------------------------------------------------------------------
if ($opt_a != 1) {
progress("Applying pressure-sensor acceleration correction...\n");
for ($ens=$firstGoodEns; $ens<=$lastGoodEns; $ens++) {
next unless numberp($CTD{W_tt}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]);
my($acorr) = $opt_a * $CTD{W_tt}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}];
if (numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH})) {
for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
$LADCP{ENSEMBLE}[$ens]->{W}[$bin] -= $acorr;
}
}
if (numberp($LADCP{ENSEMBLE}[$ens]->{REFLR_W})) {
$LADCP{ENSEMBLE}[$ens]->{REFLR_W} -= $acorr;
}
}
}
#---------------------------------------------------------------------------
# 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("Creating binned profiles at ${opt_o}m resolution...\n");
&antsAddParams('outgrid_dz',$opt_o,'outgrid_minsamp',$opt_k); # used by LADCP_w_postproc
&antsAddParams('outgrid_firstbin',$outGrid_firstBin,'outgrid_lastbin',$outGrid_lastBin);
my($min_depth) = 9e99;
my($max_depth) = 0;
local($realLastGoodEns);
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;
}
$realLastGoodEns = $ens;
my(@bindepth) = binDepths($ens);
for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
next if ($bin<$outGrid_firstBin-1 || $bin>$outGrid_lastBin-1);
$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_LADCP_w($LADCP{ENSEMBLE}[$ens]->{W}[$bin],
$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$LADCP{ENSEMBLE}[$ens]->{SPEED_OF_SOUND},
$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH},
$bindepth[$bin]);
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin] =
sscorr_LADCP_w($LADCP{ENSEMBLE}[$ens]->{W12}[$bin],
$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$LADCP{ENSEMBLE}[$ens]->{SPEED_OF_SOUND},
$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH},
$bindepth[$bin])
if numberp($LADCP{ENSEMBLE}[$ens]->{W12}[$bin]);
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin] =
sscorr_LADCP_w($LADCP{ENSEMBLE}[$ens]->{W34}[$bin],
$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$LADCP{ENSEMBLE}[$ens]->{SPEED_OF_SOUND},
$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}];
$LADCP{ENSEMBLE}[$ens]->{V12}[$bin] *= $CTD{SVEL}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]/$LADCP{ENSEMBLE}[$ens]->{SPEED_OF_SOUND}
if defined($LADCP{ENSEMBLE}[$ens]->{V12}[$bin]);
$LADCP{ENSEMBLE}[$ens]->{V34}[$bin] *= $CTD{SVEL}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]/$LADCP{ENSEMBLE}[$ens]->{SPEED_OF_SOUND}
if defined($LADCP{ENSEMBLE}[$ens]->{V34}[$bin]);
my($bi) = $bindepth[$bin]/$opt_o;
push(@{$DNCAST{ENSEMBLE}[$bi]},$ens);
push(@{$DNCAST{ELAPSED}[$bi]},$CTD{ELAPSED}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]);
push(@{$DNCAST{DEPTH}[$bi]},$bindepth[$bin]);
push(@{$DNCAST{W}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin]);
push(@{$DNCAST{W12}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin])
if defined($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin]);
push(@{$DNCAST{W34}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin])
if defined($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin]);
push(@{$DNCAST{V12}[$bi]},$LADCP{ENSEMBLE}[$ens]->{V12}[$bin])
if defined($LADCP{ENSEMBLE}[$ens]->{V12}[$bin]);
push(@{$DNCAST{V34}[$bi]},$LADCP{ENSEMBLE}[$ens]->{V34}[$bin])
if defined($LADCP{ENSEMBLE}[$ens]->{V34}[$bin]);
push(@{$DNCAST{PITCH}[$bi]},$LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH});
push(@{$DNCAST{ROLL}[$bi]},$LADCP{ENSEMBLE}[$ens]->{ROLL});
push(@{$DNCAST{TILT}[$bi]},$LADCP{ENSEMBLE}[$ens]->{TILT});
}
}
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{MEDIAN_W12}[$bi] = median(@{$DNCAST{W12}[$bi]});
$DNCAST{MEDIAN_W34}[$bi] = median(@{$DNCAST{W34}[$bi]});
$DNCAST{MEDIAN_V12}[$bi] = median(@{$DNCAST{V12}[$bi]});
$DNCAST{MEDIAN_V34}[$bi] = median(@{$DNCAST{V34}[$bi]});
$DNCAST{MEAN_PITCH}[$bi] = avg(@{$DNCAST{PITCH}[$bi]});
$DNCAST{MEAN_ROLL}[$bi] = avg(@{$DNCAST{ROLL}[$bi]});
$DNCAST{MEAN_TILT}[$bi] = avg(@{$DNCAST{TILT}[$bi]});
$DNCAST{MAD_W}[$bi] = mad2($DNCAST{MEDIAN_W}[$bi],@{$DNCAST{W}[$bi]});
$DNCAST{N_SAMP}[$bi] = @{$DNCAST{W}[$bi]};
}
progress("\tupcast...\n"); # upcast
for ($ens=$LADCP_atbottom; $ens<=$lastGoodEns; $ens++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
$realLastGoodEns = $ens;
my(@bindepth) = binDepths($ens);
for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
next if ($bin<$outGrid_firstBin-1 || $bin>$outGrid_lastBin-1);
$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_LADCP_w($LADCP{ENSEMBLE}[$ens]->{W}[$bin],
$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$LADCP{ENSEMBLE}[$ens]->{SPEED_OF_SOUND},
$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH},
$bindepth[$bin]);
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin] =
sscorr_LADCP_w($LADCP{ENSEMBLE}[$ens]->{W12}[$bin],
$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$LADCP{ENSEMBLE}[$ens]->{SPEED_OF_SOUND},
$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH},
$bindepth[$bin])
if numberp($LADCP{ENSEMBLE}[$ens]->{W12}[$bin]);
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin] =
sscorr_LADCP_w($LADCP{ENSEMBLE}[$ens]->{W34}[$bin],
$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$LADCP{ENSEMBLE}[$ens]->{SPEED_OF_SOUND},
$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}];
$LADCP{ENSEMBLE}[$ens]->{V12}[$bin] *= $CTD{SVEL}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]/$LADCP{ENSEMBLE}[$ens]->{SPEED_OF_SOUND}
if defined($LADCP{ENSEMBLE}[$ens]->{V12}[$bin]);
$LADCP{ENSEMBLE}[$ens]->{V34}[$bin] *= $CTD{SVEL}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]/$LADCP{ENSEMBLE}[$ens]->{SPEED_OF_SOUND}
if defined($LADCP{ENSEMBLE}[$ens]->{V34}[$bin]);
my($bi) = $bindepth[$bin]/$opt_o;
push(@{$UPCAST{ENSEMBLE}[$bi]},$ens);
push(@{$UPCAST{ELAPSED}[$bi]},$CTD{ELAPSED}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]);
push(@{$UPCAST{DEPTH}[$bi]},$bindepth[$bin]);
push(@{$UPCAST{W}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin]);
push(@{$UPCAST{W12}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin])
if defined($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin]);
push(@{$UPCAST{W34}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin])
if defined($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin]);
push(@{$UPCAST{V12}[$bi]},$LADCP{ENSEMBLE}[$ens]->{V12}[$bin])
if defined($LADCP{ENSEMBLE}[$ens]->{V12}[$bin]);
push(@{$UPCAST{V34}[$bi]},$LADCP{ENSEMBLE}[$ens]->{V34}[$bin])
if defined($LADCP{ENSEMBLE}[$ens]->{V34}[$bin]);
push(@{$UPCAST{PITCH}[$bi]},$LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH});
push(@{$UPCAST{ROLL}[$bi]},$LADCP{ENSEMBLE}[$ens]->{ROLL});
push(@{$UPCAST{TILT}[$bi]},$LADCP{ENSEMBLE}[$ens]->{TILT});
}
}
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{MEDIAN_W12}[$bi] = median(@{$UPCAST{W12}[$bi]});
$UPCAST{MEDIAN_W34}[$bi] = median(@{$UPCAST{W34}[$bi]});
$UPCAST{MEDIAN_V12}[$bi] = median(@{$UPCAST{V12}[$bi]});
$UPCAST{MEDIAN_V34}[$bi] = median(@{$UPCAST{V34}[$bi]});
$UPCAST{MEAN_PITCH}[$bi] = avg(@{$UPCAST{PITCH}[$bi]});
$UPCAST{MEAN_ROLL}[$bi] = avg(@{$UPCAST{ROLL}[$bi]});
$UPCAST{MEAN_TILT}[$bi] = avg(@{$UPCAST{TILT}[$bi]});
$UPCAST{MAD_W}[$bi] = mad2($UPCAST{MEDIAN_W}[$bi],@{$UPCAST{W}[$bi]});
$UPCAST{N_SAMP}[$bi] = @{$UPCAST{W}[$bi]};
}
#------------------------------------------------------------------------------------------------------
# remove ping-coherent residuals
# - potential error sources:
# 1) acceleration-dependence of Paroscientific pressure measurements; O(10cm/s) [IWISE 28]
# 2) residual CTD/LADCP mismatch errors; O(1cm/s) [Thurnherr, CWTMC 2011]
# 3) ADCP short-term variability; O(1cm/s) for vertical?
# - multiple sub-steps:
# - remove ensemble-median residual from all ocean vertical velocities
# - re-bin profiles
# - calculate new residuals
# - improves solution but can introduce errors in ensembles with strong outliers => median instead of
# mean, and outlier filter above
#------------------------------------------------------------------------------------------------------
unless ($opt_q) {
progress("Subtracting ping-coherent residuals...\n");
for ($ens=$firstGoodEns; $ens<=$lastGoodEns; $ens++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
my(@residuals) = (); # calc median of residuals
my(@bindepth) = binDepths($ens);
for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
next if ($bin<$outGrid_firstBin-1 || $bin>$outGrid_lastBin-1);
my($bi) = $bindepth[$bin]/$opt_o;
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); # NB: can be nan!
next unless numberp($LADCP{ENSEMBLE}[$ens]->{MEDIAN_RESIDUAL_W});
for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) { # subtract from ocean velocities
next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin] -=
$LADCP{ENSEMBLE}[$ens]->{MEDIAN_RESIDUAL_W}
if numberp($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin]);
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin] -=
$LADCP{ENSEMBLE}[$ens]->{MEDIAN_RESIDUAL_W}
if numberp($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin]);
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin] -=
$LADCP{ENSEMBLE}[$ens]->{MEDIAN_RESIDUAL_W}
if numberp($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin]);
}
$LADCP{ENSEMBLE}[$ens]->{REFLR_W} -= # NB: this can be nan here
$LADCP{ENSEMBLE}[$ens]->{MEDIAN_RESIDUAL_W}
if numberp($LADCP{ENSEMBLE}[$ens]->{REFLR_W});
}
#------------------------------------------------------------
progress("\tre-binning profile data...\n");
for (my($bi)=0; $bi<=$#{$DNCAST{ENSEMBLE}}; $bi++) { # bin data
for (my($i)=0; $i<@{$DNCAST{W}[$bi]}; $i++) { # code works if MEDIAN_RESIDUAL_W is nan (possible?)
$DNCAST{W} [$bi][$i] -= $LADCP{ENSEMBLE}[$DNCAST{ENSEMBLE}[$bi][$i]]->{MEDIAN_RESIDUAL_W}
if numberp($DNCAST{W}[$bi][$i]);
$DNCAST{W12}[$bi][$i] -= $LADCP{ENSEMBLE}[$DNCAST{ENSEMBLE}[$bi][$i]]->{MEDIAN_RESIDUAL_W}
if numberp($DNCAST{W12}[$bi][$i]);
$DNCAST{W34}[$bi][$i] -= $LADCP{ENSEMBLE}[$DNCAST{ENSEMBLE}[$bi][$i]]->{MEDIAN_RESIDUAL_W}
if numberp($DNCAST{W34}[$bi][$i]);
}
$DNCAST{MEDIAN_W} [$bi] = median(@{$DNCAST{W}[$bi]}) if numberp($DNCAST{MEDIAN_W} [$bi]);
$DNCAST{MEDIAN_W12}[$bi] = median(@{$DNCAST{W12}[$bi]}) if numberp($DNCAST{MEDIAN_W12}[$bi]);
$DNCAST{MEDIAN_W34}[$bi] = median(@{$DNCAST{W34}[$bi]}) if numberp($DNCAST{MEDIAN_W34}[$bi]);
$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}
if numberp($UPCAST{W}[$bi][$i]);
$UPCAST{W12}[$bi][$i] -= $LADCP{ENSEMBLE}[$UPCAST{ENSEMBLE}[$bi][$i]]->{MEDIAN_RESIDUAL_W}
if numberp($UPCAST{W12}[$bi][$i]);
$UPCAST{W34}[$bi][$i] -= $LADCP{ENSEMBLE}[$UPCAST{ENSEMBLE}[$bi][$i]]->{MEDIAN_RESIDUAL_W}
if numberp($UPCAST{W34}[$bi][$i]);
}
$UPCAST{MEDIAN_W} [$bi] = median(@{$UPCAST{W}[$bi]}) if numberp($UPCAST{MEDIAN_W} [$bi]);
$UPCAST{MEDIAN_W12}[$bi] = median(@{$UPCAST{W12}[$bi]}) if numberp($UPCAST{MEDIAN_W12}[$bi]);
$UPCAST{MEDIAN_W34}[$bi] = median(@{$UPCAST{W34}[$bi]}) if numberp($UPCAST{MEDIAN_W34}[$bi]);
$UPCAST{MAD_W} [$bi] = mad2($UPCAST{MEDIAN_W}[$bi],@{$UPCAST{W}[$bi]});
}
} # unless ($opt_q);
#----------------------------------------------------------------------
# remove ensembles with large rms residuals
#----------------------------------------------------------------------
if (defined($opt_r) && !$opt_q) {
progress("Applying residuals filters...\n");
progress("\trms residuals > $residuals_rms_max: ");
my($nerm) = editResiduals_rmsMax($firstGoodEns,$lastGoodEns,$residuals_rms_max);
progress("$nerm ensembles removed (%d%% of total)\n",round(100*$nerm/($lastGoodEns-$firstGoodEns+1)));
if (defined($residuals_delta_max)) { # filter based on difference between
progress("\tbeampair residuals difference > $residuals_delta_max [m/s]: ");
my($nvrm) = editResiduals_deltaMax($firstGoodEns,$lastGoodEns,$residuals_delta_max);
progress("$nvrm velocities removed (%d%% of total in bins $LADCP_firstBin-$LADCP_lastBin)\n",round(100*$nvrm/$nvw));
}
progress("\tre-binning profile data...\n");
undef(%DNCAST); undef(%UPCAST);
$min_depth = 9e99; my($max_depth) = 0;
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;
}
$realLastGoodEns = $ens;
my(@bindepth) = binDepths($ens);
for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
next if ($bin<$outGrid_firstBin-1 || $bin>$outGrid_lastBin-1);
$min_depth = $bindepth[$bin] if ($bindepth[$bin] < $min_depth);
$max_depth = $bindepth[$bin] if ($bindepth[$bin] > $max_depth);
my($bi) = $bindepth[$bin]/$opt_o;
push(@{$DNCAST{ENSEMBLE}[$bi]},$ens);
push(@{$DNCAST{ELAPSED}[$bi]},$CTD{ELAPSED}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]);
push(@{$DNCAST{DEPTH}[$bi]},$bindepth[$bin]);
push(@{$DNCAST{W}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin]);
push(@{$DNCAST{W12}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin])
if defined($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin]);
push(@{$DNCAST{W34}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin])
if defined($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin]);
push(@{$DNCAST{V12}[$bi]},$LADCP{ENSEMBLE}[$ens]->{V12}[$bin])
if defined($LADCP{ENSEMBLE}[$ens]->{V12}[$bin]);
push(@{$DNCAST{V34}[$bi]},$LADCP{ENSEMBLE}[$ens]->{V34}[$bin])
if defined($LADCP{ENSEMBLE}[$ens]->{V34}[$bin]);
push(@{$DNCAST{PITCH}[$bi]},$LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH});
push(@{$DNCAST{ROLL}[$bi]},$LADCP{ENSEMBLE}[$ens]->{ROLL});
push(@{$DNCAST{TILT}[$bi]},$LADCP{ENSEMBLE}[$ens]->{TILT});
}
}
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{MEDIAN_W12}[$bi] = median(@{$DNCAST{W12}[$bi]});
$DNCAST{MEDIAN_W34}[$bi] = median(@{$DNCAST{W34}[$bi]});
$DNCAST{MEDIAN_V12}[$bi] = median(@{$DNCAST{V12}[$bi]});
$DNCAST{MEDIAN_V34}[$bi] = median(@{$DNCAST{V34}[$bi]});
$DNCAST{MEAN_PITCH}[$bi] = avg(@{$DNCAST{PITCH}[$bi]});
$DNCAST{MEAN_ROLL}[$bi] = avg(@{$DNCAST{ROLL}[$bi]});
$DNCAST{MEAN_TILT}[$bi] = avg(@{$DNCAST{TILT}[$bi]});
$DNCAST{MAD_W}[$bi] = mad2($DNCAST{MEDIAN_W}[$bi],@{$DNCAST{W}[$bi]});
$DNCAST{N_SAMP}[$bi] = @{$DNCAST{W}[$bi]};
}
for ($ens=$LADCP_atbottom; $ens<=$lastGoodEns; $ens++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
$realLastGoodEns = $ens;
my(@bindepth) = binDepths($ens);
for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
next if ($bin<$outGrid_firstBin-1 || $bin>$outGrid_lastBin-1);
$min_depth = $bindepth[$bin] if ($bindepth[$bin] < $min_depth);
$max_depth = $bindepth[$bin] if ($bindepth[$bin] > $max_depth);
my($bi) = $bindepth[$bin]/$opt_o;
push(@{$UPCAST{ENSEMBLE}[$bi]},$ens);
push(@{$UPCAST{ELAPSED}[$bi]},$CTD{ELAPSED}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]);
push(@{$UPCAST{DEPTH}[$bi]},$bindepth[$bin]);
push(@{$UPCAST{W}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin]);
push(@{$UPCAST{W12}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin])
if defined($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin]);
push(@{$UPCAST{W34}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin])
if defined($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin]);
push(@{$UPCAST{V12}[$bi]},$LADCP{ENSEMBLE}[$ens]->{V12}[$bin])
if defined($LADCP{ENSEMBLE}[$ens]->{V12}[$bin]);
push(@{$UPCAST{V34}[$bi]},$LADCP{ENSEMBLE}[$ens]->{V34}[$bin])
if defined($LADCP{ENSEMBLE}[$ens]->{V34}[$bin]);
push(@{$UPCAST{PITCH}[$bi]},$LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH});
push(@{$UPCAST{ROLL}<[$bi]},$LADCP{ENSEMBLE}[$ens]->{ROLL});
push(@{$UPCAST{TILT}[$bi]},$LADCP{ENSEMBLE}[$ens]->{TILT});
}
}
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{MEDIAN_W12}[$bi] = median(@{$UPCAST{W12}[$bi]});
$UPCAST{MEDIAN_W34}[$bi] = median(@{$UPCAST{W34}[$bi]});
$UPCAST{MEDIAN_V12}[$bi] = median(@{$UPCAST{V12}[$bi]});
$UPCAST{MEDIAN_V34}[$bi] = median(@{$UPCAST{V34}[$bi]});
$UPCAST{MEAN_PITCH}[$bi] = avg(@{$UPCAST{PITCH}[$bi]});
$UPCAST{MEAN_ROLL}[$bi] = avg(@{$UPCAST{ROLL}[$bi]});
$UPCAST{MEAN_TILT}[$bi] = avg(@{$UPCAST{TILT}[$bi]});
$UPCAST{MAD_W}[$bi] = mad2($UPCAST{MEDIAN_W}[$bi],@{$UPCAST{W}[$bi]});
$UPCAST{N_SAMP}[$bi] = @{$UPCAST{W}[$bi]};
}
} # if defined($opt_r)
$lastGoodEns = $realLastGoodEns; # now the ensemble range is final
&antsAddParams('depth.min',$min_depth,'depth.max',$max_depth, # plot range limits
'ens.min',$LADCP{ENSEMBLE}[$firstGoodEns]->{NUMBER},
'ens.max',$LADCP{ENSEMBLE}[$lastGoodEns]->{NUMBER},
'elapsed.min',$LADCP{ENSEMBLE}[$firstGoodEns]->{CTD_ELAPSED},
'elapsed.max',$LADCP{ENSEMBLE}[$lastGoodEns]->{CTD_ELAPSED});
#-----------------------------------------------
# Add residuals to depth-binned profiles (wprof)
#-----------------------------------------------
progress("Binning vertical-velocity residuals into depth profiles...");
for ($ens=$firstGoodEns; $ens<$LADCP_atbottom; $ens++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
my(@bindepth) = binDepths($ens);
for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
next unless ($bin+1>=$outGrid_firstBin && $bin+1<=$outGrid_lastBin);
next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
my($bi) = $bindepth[$bin]/$opt_o;
next unless numberp($DNCAST{MEDIAN_W}[$bi]);
push(@{$DNCAST{RESIDUAL12}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin]-$DNCAST{MEDIAN_W}[$bi])
if numberp($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin]);
push(@{$DNCAST{RESIDUAL34}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin]-$DNCAST{MEDIAN_W}[$bi])
if numberp($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin]);
}
}
for (my($bi)=0; $bi<=$#{$DNCAST{ENSEMBLE}}; $bi++) {
$DNCAST{MEAN_RESIDUAL12}[$bi] = avg(@{$DNCAST{RESIDUAL12}[$bi]});
$DNCAST{MEAN_RESIDUAL34}[$bi] = avg(@{$DNCAST{RESIDUAL34}[$bi]});
}
for (my($bi)=0; $bi<$#{$DNCAST{ENSEMBLE}}; $bi++) { # rms beampair residual in 5-bin-thick layers
$DNCAST{LR_RMS_BP_RESIDUAL}[$bi] =
rms(@{$DNCAST{MEAN_RESIDUAL12}}[max(0,$bi-2)..min($#{$DNCAST{ENSEMBLE}},$bi+2)],
@{$DNCAST{MEAN_RESIDUAL34}}[max(0,$bi-2)..min($#{$DNCAST{ENSEMBLE}},$bi+2)]);
}
for ($ens=$LADCP_atbottom; $ens<=$lastGoodEns; $ens++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
my(@bindepth) = binDepths($ens);
for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
next unless ($bin+1>=$outGrid_firstBin && $bin+1<=$outGrid_lastBin);
next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
my($bi) = $bindepth[$bin]/$opt_o;
next unless numberp($UPCAST{MEDIAN_W}[$bi]);
push(@{$UPCAST{RESIDUAL12}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin]-$UPCAST{MEDIAN_W}[$bi])
if numberp($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin]);
push(@{$UPCAST{RESIDUAL34}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin]-$UPCAST{MEDIAN_W}[$bi])
if numberp($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin]);
}
}
for (my($bi)=0; $bi<=$#{$UPCAST{ENSEMBLE}}; $bi++) {
$UPCAST{MEAN_RESIDUAL12}[$bi] = avg(@{$UPCAST{RESIDUAL12}[$bi]});
$UPCAST{MEAN_RESIDUAL34}[$bi] = avg(@{$UPCAST{RESIDUAL34}[$bi]});
}
for (my($bi)=0; $bi<$#{$UPCAST{ENSEMBLE}}; $bi++) { # rms beampair residual in 5-bin-thick layers
$UPCAST{LR_RMS_BP_RESIDUAL}[$bi] =
rms(@{$UPCAST{MEAN_RESIDUAL12}}[max(0,$bi-2)..min($#{$UPCAST{ENSEMBLE}},$bi+2)],
@{$UPCAST{MEAN_RESIDUAL34}}[max(0,$bi-2)..min($#{$UPCAST{ENSEMBLE}},$bi+2)]);
}
progress("\n");
#-----------------------------------------------
# Apply high-layer-residuals filter
#-----------------------------------------------
progress("Applying residuals profile filter...\n");
&antsAddParams('layer_residuals_rms.max',$layer_residuals_rms_max);
progress("\t5-bin rms residuals > $layer_residuals_rms_max: ");
my($nbrm) = editHighResidualLayers($layer_residuals_rms_max);
my($nbrmf) = 0;
$nbrmf = $nbrm / (@{$DNCAST{LR_RMS_BP_RESIDUAL}}+@{$UPCAST{LR_RMS_BP_RESIDUAL}})
if ($nbrm > 0);
progress("$nbrm bins removed (%d%% of profile)\n",round(100*$nbrmf));
warning(2,"large fraction (%d%%) of profile exceeds residuals-profile threshold\n",round(100*$nbrmf))
if ($nbrmf >= 0.1);
#------------------------
# Profile is now final
# => calc <w>, w.var
#------------------------
my($dcw,$ucw,$nw) = (0,0); # <w> (vertically averaged)
for (my($bi)=0; $bi<@{$DNCAST{ENSEMBLE}}; $bi++) {
next unless ($DNCAST{N_SAMP}[$bi]>=$opt_k) && numberp($DNCAST{MEDIAN_W}[$bi]);
$dcw += $DNCAST{MEDIAN_W}[$bi];
$nw++;
}
$dcw = $nw ? ($dcw/$nw) : nan;
$nw = 0;
for (my($bi)=0; $bi<@{$UPCAST{ENSEMBLE}}; $bi++) {
next unless ($UPCAST{N_SAMP}[$bi]>=$opt_k) && numberp($UPCAST{MEDIAN_W}[$bi]);
$ucw += $UPCAST{MEDIAN_W}[$bi];
$nw++;
}
$ucw = $nw ? ($ucw/$nw) : nan;
my($dcss,$ucss,$nw) = (0,0); # variance(w)
for (my($bi)=0; $bi<@{$DNCAST{ENSEMBLE}}; $bi++) {
next unless ($DNCAST{N_SAMP}[$bi]>=$opt_k) && numberp($DNCAST{MEDIAN_W}[$bi]);
$dcss += ($DNCAST{MEDIAN_W}[$bi] - $dcw)**2;
$nw++;
}
$dcss = ($nw > 1) ? ($dcss/$nw) : nan;
$nw = 0;
for (my($bi)=0; $bi<@{$UPCAST{ENSEMBLE}}; $bi++) {
next unless ($UPCAST{N_SAMP}[$bi]>=$opt_k) && numberp($UPCAST{MEDIAN_W}[$bi]);
$ucss += ($UPCAST{MEDIAN_W}[$bi] - $ucw)**2;
$nw++;
}
$ucss = ($nw > 1) ? ($ucss/$nw) : nan;
&antsAddParams('dc_w.mu',$dcw,'dc_w.var',$dcss,'uc_w.mu',$ucw,'uc_w.var',$ucss);
#-------------------------------------------------
# Calculate and Output ADCP-Bin-Averaged Residuals
#-------------------------------------------------
progress("Binning vertical-velocity residuals wrt. ADCP bin number...");
local(@dc_bres12,@uc_bres12,@dc_bres34,@uc_bres34);
for ($ens=$firstGoodEns; $ens<$LADCP_atbottom; $ens++) { # DOWNCAST
next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
my(@bindepth) = binDepths($ens);
for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
next unless numberp($DNCAST{MEDIAN_W}[$bindepth[$bin]/$opt_o]);
push(@{$dc_bres12[$bin]},
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin]
- $DNCAST{MEDIAN_W12}[$bindepth[$bin]/$opt_o])
if numberp($DNCAST{MEDIAN_W12}[$bindepth[$bin]/$opt_o]);
push(@{$dc_bres34[$bin]},
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin]
- $DNCAST{MEDIAN_W34}[$bindepth[$bin]/$opt_o])
if numberp($DNCAST{MEDIAN_W34}[$bindepth[$bin]/$opt_o]);
}
}
for ($ens=$LADCP_atbottom; $ens<=$lastGoodEns; $ens++) { # UPCAST
next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
my(@bindepth) = binDepths($ens);
for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
next unless numberp($UPCAST{MEDIAN_W}[$bindepth[$bin]/$opt_o]);
push(@{$uc_bres12[$bin]},
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin]
- $UPCAST{MEDIAN_W12}[$bindepth[$bin]/$opt_o])
if numberp($UPCAST{MEDIAN_W12}[$bindepth[$bin]/$opt_o]);
push(@{$uc_bres34[$bin]},
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin]
- $UPCAST{MEDIAN_W34}[$bindepth[$bin]/$opt_o])
if numberp($UPCAST{MEDIAN_W34}[$bindepth[$bin]/$opt_o]);
}
}
local(@dc_avg_bres12,@uc_avg_bres12,@dc_avg_bres12_nsamp,@uc_avg_bres12_nsamp); # means/stddevs
local(@dc_avg_bres34,@uc_avg_bres34,@dc_avg_bres34_nsamp,@uc_avg_bres34_nsamp);
for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
$dc_bres12_nsamp[$bin] = @{$dc_bres12[$bin]};
if ($dc_bres12_nsamp[$bin] > 0) {
$dc_avg_bres12[$bin] = avg(@{$dc_bres12[$bin]});
$dc_sig_bres12[$bin] = stddev2($dc_avg_bres12[$bin],@{$dc_bres12[$bin]});
} else {
$dc_avg_bres12[$bin] = nan;
$dc_sig_bres12[$bin] = nan;
}
$uc_bres12_nsamp[$bin] = @{$uc_bres12[$bin]};
if ($uc_bres12_nsamp[$bin] > 0) {
$uc_avg_bres12[$bin] = avg(@{$uc_bres12[$bin]});
$uc_sig_bres12[$bin] = stddev2($uc_avg_bres12[$bin],@{$uc_bres12[$bin]});
} else {
$uc_avg_bres12[$bin] = nan;
$uc_sig_bres12[$bin] = nan;
}
$dc_bres34_nsamp[$bin] = @{$dc_bres34[$bin]};
if ($dc_bres34_nsamp[$bin] > 0) {
$dc_avg_bres34[$bin] = avg(@{$dc_bres34[$bin]});
$dc_sig_bres34[$bin] = stddev2($dc_avg_bres34[$bin],@{$dc_bres34[$bin]});
} else {
$dc_avg_bres34[$bin] = nan;
$dc_sig_bres34[$bin] = nan;
}
$uc_bres34_nsamp[$bin] = @{$uc_bres34[$bin]};
if ($uc_bres34_nsamp[$bin] > 0) {
$uc_avg_bres34[$bin] = avg(@{$uc_bres34[$bin]});
$uc_sig_bres34[$bin] = stddev2($uc_avg_bres34[$bin],@{$uc_bres34[$bin]});
} else {
$uc_avg_bres34[$bin] = nan;
$uc_sig_bres34[$bin] = nan;
}
}
local($uc_bres12_max_nsamp,$dc_bres12_max_nsamp) = (0,0); # calc rms in block of well sampled bins
local($uc_bres34_max_nsamp,$dc_bres34_max_nsamp) = (0,0);
for (my($bin)=$LADCP_firstBin; $bin<=$LADCP_lastBin-1; $bin++) { # SKIP 1ST BIN!!!
next if ($bin+1<$outGrid_firstBin || $bin+1>$outGrid_lastBin); # skip bins not included in gridded output
$dc_bres12_max_nsamp = $dc_bres12_nsamp[$bin] # nsamp in best sampled bin
if ($dc_bres12_nsamp[$bin] > $dc_bres12_max_nsamp);
$uc_bres12_max_nsamp = $uc_bres12_nsamp[$bin]
if ($uc_bres12_nsamp[$bin] > $uc_bres12_max_nsamp);
$dc_bres34_max_nsamp = $dc_bres34_nsamp[$bin] # nsamp in best sampled bin
if ($dc_bres34_nsamp[$bin] > $dc_bres34_max_nsamp);
$uc_bres34_max_nsamp = $uc_bres34_nsamp[$bin]
if ($uc_bres34_nsamp[$bin] > $uc_bres34_max_nsamp);
}
my($dc_sumsq12,$uc_sumsq12,$dc_n12,$uc_n12) = (0,0,0,0); # calc rms residuals
my($dc_sumsq34,$uc_sumsq34,$dc_n34,$uc_n34) = (0,0,0,0);
for (my($bin)=$LADCP_firstBin; $bin<=$LADCP_lastBin-1; $bin++) { # SKIP 1ST BIN
next if ($bin+1<$outGrid_firstBin || $bin+1>$outGrid_lastBin); # skip bins not included in gridded output
if ($dc_bres12_nsamp[$bin] >= $dc_bres12_max_nsamp/3) { # skip bins with < 1/3 max(nsamp)
$dc_sumsq12 += $dc_avg_bres12[$bin]**2;
$dc_n12++;
}
if ($uc_bres12_nsamp[$bin] >= $uc_bres12_max_nsamp/3) {
$uc_sumsq12 += $uc_avg_bres12[$bin]**2;
$uc_n12++;
}
if ($dc_bres34_nsamp[$bin] >= $dc_bres34_max_nsamp/3) { # skip bins with < 1/3 max(nsamp)
$dc_sumsq34 += $dc_avg_bres34[$bin]**2;
$dc_n34++;
}
if ($uc_bres34_nsamp[$bin] >= $uc_bres34_max_nsamp/3) {
$uc_sumsq34 += $uc_avg_bres34[$bin]**2;
$uc_n34++;
}
}
$dc_bres12_rms = ($dc_n12 > 0) ? sqrt($dc_sumsq12/$dc_n12) : nan;
$uc_bres12_rms = ($uc_n12 > 0) ? sqrt($uc_sumsq12/$uc_n12) : nan;
$dc_bres34_rms = ($dc_n34 > 0) ? sqrt($dc_sumsq34/$dc_n34) : nan;
$uc_bres34_rms = ($uc_n34 > 0) ? sqrt($uc_sumsq34/$uc_n34) : nan;
&antsAddParams('dc_bin_residuals12.rms',$dc_bres12_rms,
'dc_bin_residuals34.rms',$dc_bres34_rms,
'uc_bin_residuals12.rms',$uc_bres12_rms,
'uc_bin_residuals34.rms',$uc_bres34_rms);
progress("\n\twriting binned residuals to ");
my($saveParams) = $antsCurParams;
@antsNewLayout = ('bin','dc_residual12','dc_residual12.sig','dc_residual12.nsamp',
'uc_residual12','uc_residual12.sig','uc_residual12.nsamp',
'dc_residual34','dc_residual34.sig','dc_residual34.nsamp',
'uc_residual34','uc_residual34.sig','uc_residual34.nsamp');
&antsAddParams('BR_bin.max',max(scalar(@dc_bres12),scalar(@uc_bres12),
scalar(@dc_bres34),scalar(@uc_bres34)));
foreach my $of (@out_BR) {
progress("<$of> ");
my($sub,$arg) = ($of =~ /^([^\(]+)\(([^\)]+)\)$/); # plot_sub(out_file)
if (defined($arg)) {
require "$WCALC/${sub}.pl";
&{$sub}($arg);
next;
}
$of = ">$of" unless ($of =~ /^$|^\s*\|/);
open(STDOUT,$of) || error("$of: $!\n");
undef($antsActiveHeader) unless ($ANTS_TOOLS_AVAILABLE);
for (my($bin)=0; $bin<max(scalar(@dc_bres12),scalar(@uc_bres12)); $bin++) {
&antsOut($bin+1,$dc_avg_bres12[$bin],$dc_sig_bres12[$bin],$dc_bres12_nsamp[$bin],
$uc_avg_bres12[$bin],$uc_sig_bres12[$bin],$uc_bres12_nsamp[$bin],
$dc_avg_bres34[$bin],$dc_sig_bres34[$bin],$dc_bres34_nsamp[$bin],
$uc_avg_bres34[$bin],$uc_sig_bres34[$bin],$uc_bres34_nsamp[$bin]);
}
&antsOut('EOF'); open(STDOUT,">&2");
}
$antsCurParams = $saveParams;
progress("\n");
#--------------------------------------------------
# Calculate BT-referenced vertical-velocity profile
#--------------------------------------------------
if ($LADCP{ENSEMBLE}[$LADCP_atbottom]->{XDUCER_FACING_DOWN} && defined($water_depth)) {
progress("Calculating BT-referenced vertical velocities\n");
calc_BTprof($firstGoodEns,$lastGoodEns,$water_depth,$sig_water_depth);
my($sumSq) = my($n) = 0;
for (my($bi)=0; $bi<=$#{$BT{MEDIAN_W}}; $bi++) {
next unless numbersp($BT{MEDIAN_W}[$bi],$DNCAST{MEDIAN_W}[$bi],$UPCAST{MEDIAN_W}[$bi]);
next unless ($BT{N_SAMP}[$bi]>=$opt_k && $DNCAST{N_SAMP}[$bi]>=$opt_k && $UPCAST{N_SAMP}[$bi]>=$opt_k);
$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_w_discrepancy.rms',$rms);
progress("\t$rms m/s rms vertical-velocity discrepancy\n");
}
}
#------------------------
# output all samples (-w)
#------------------------
# NB: residual fields are calculated with respect to down-/upcast medians in -o-size bins
if (@out_wsamp) {
progress("Writing vertical-velocity data to ");
foreach my $of (@out_wsamp) {
progress("<$of> ");
my($sub,$arg) = ($of =~ /^([^\(]+)\(([^\)]+)\)$/); # plot_sub(out_file)
if (defined($arg)) {
require "$WCALC/${sub}.pl";
&{$sub}($arg);
next;
}
@antsNewLayout = ('ensemble','bin','elapsed','depth','CTD_depth','downcast',
'w','w12','w34','v12','v34','residual','residual12','residual34',
'CTD_w','CTD_w_t','CTD_w_tt','LADCP_w','LADCP_reflr_w','winch_w',
'errvel','correlation','echo_amplitude','Sv',
'pitch','roll','tilt','heading','3_beam','svel');
$of = ">$of" unless ($of =~ /^$|^\s*\|/);
open(STDOUT,$of) || error("$of: $!\n");
undef($antsActiveHeader) unless ($ANTS_TOOLS_AVAILABLE);
sub res($$)
{
my($meas,$mean) = @_;
return numberp($meas) ? ($meas - $mean) : undef;
}
for ($ens=$firstGoodEns; $ens<$LADCP_atbottom; $ens++) { # downcast
next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
my(@bindepth) = binDepths($ens);
for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
my($bi) = $bindepth[$bin]/$opt_o;
next unless numberp($DNCAST{MEDIAN_W}[$bi]);
&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{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin],
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin],
$LADCP{ENSEMBLE}[$ens]->{V12}[$bin],$LADCP{ENSEMBLE}[$ens]->{V34}[$bin],
res($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin],$DNCAST{MEDIAN_W}[$bi]),
res($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin],$DNCAST{MEDIAN_W}[$bi]),
res($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin],$DNCAST{MEDIAN_W}[$bi]),
$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$CTD{W_t}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$CTD{W_tt}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_W}[$bin],
$LADCP{ENSEMBLE}[$ens]->{REFLR_W},
$CTD{W_WINCH}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$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]->{GIMBAL_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
for ($ens=$LADCP_atbottom; $ens<=$lastGoodEns; $ens++) { # upcast
next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
my(@bindepth) = binDepths($ens);
for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
my($bi) = $bindepth[$bin]/$opt_o;
next unless numberp($UPCAST{MEDIAN_W}[$bi]);
&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{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin],
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin],
$LADCP{ENSEMBLE}[$ens]->{V12}[$bin],$LADCP{ENSEMBLE}[$ens]->{V34}[$bin],
res($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin],$UPCAST{MEDIAN_W}[$bi]),
res($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin],$UPCAST{MEDIAN_W}[$bi]),
res($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin],$UPCAST{MEDIAN_W}[$bi]),
$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$CTD{W_t}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$CTD{W_tt}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_W}[$bin],
$LADCP{ENSEMBLE}[$ens]->{REFLR_W},
$CTD{W_WINCH}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$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]->{GIMBAL_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
&antsOut('EOF'); open(STDOUT,">&2");
}
progress("\n");
}
#----------------------------
# Output depth-binned profile
#----------------------------
if (@out_profile) {
progress("Writing vertical-velocity profiles to ");
undef(@antsNewLayout);
push(@antsNewLayout,'depth','hab', # consistent with [LADCP_w_postproc]
'dc_elapsed','dc_w','dc_w.mad','dc_w.nsamp',
'uc_elapsed','uc_w','uc_w.mad','uc_w.nsamp');
push(@antsNewLayout,'dc_lr_bp_res.rms', # additional fields (inconsistent w postproc)
'dc_depth','dc_w12','dc_w34','dc_v12','dc_v34',
'uc_lr_bp_res.rms',
'uc_depth','uc_w12','uc_w34','uc_v12','uc_v34',
'dc_pitch.mu','dc_roll.mu','dc_tilt.mu',
'uc_pitch.mu','uc_roll.mu','uc_tilt.mu',
'BT_w','BT_w.mad','BT_w.nsamp');
foreach my $of (@out_profile) {
progress("<$of> ");
my($sub,$arg) = ($of =~ /^([^\(]+)\(([^\)]+)\)$/); # plot_sub(out_file)
if (defined($arg)) {
require "$WCALC/${sub}.pl";
&{$sub}($arg);
next;
}
$of = ">$of" unless ($of =~ /^$|^\s*\|/); # pipe or file output
open(STDOUT,$of) || error("$of: $!\n");
undef($antsActiveHeader) unless ($ANTS_TOOLS_AVAILABLE);
for (my($bi)=0; $bi<=max($#{$DNCAST{ENSEMBLE}},$#{$UPCAST{ENSEMBLE}},$#{$BT{NSAMP}}); $bi++) {
&antsOut(($bi+0.5)*$opt_o, # nominal depth
defined($water_depth)?$water_depth-($bi+0.5)*$opt_o:nan, # nominal hab
$DNCAST{MEAN_ELAPSED}[$bi], # dc data consistent w postproc
$DNCAST{N_SAMP}[$bi]>=$opt_k?$DNCAST{MEDIAN_W}[$bi]:nan,
$DNCAST{MAD_W}[$bi],$DNCAST{N_SAMP}[$bi],
$UPCAST{MEAN_ELAPSED}[$bi], # uc data consistent w postproc
$UPCAST{N_SAMP}[$bi]>=$opt_k?$UPCAST{MEDIAN_W}[$bi]:nan,
$UPCAST{MAD_W}[$bi],$UPCAST{N_SAMP}[$bi],
$DNCAST{LR_RMS_BP_RESIDUAL}[$bi], # remaining dc data
$DNCAST{MEAN_DEPTH}[$bi],
$DNCAST{MEDIAN_W12}[$bi],$DNCAST{MEDIAN_W34}[$bi],
$DNCAST{MEDIAN_V12}[$bi],$DNCAST{MEDIAN_V34}[$bi],
$UPCAST{LR_RMS_BP_RESIDUAL}[$bi], # remaining uc data
$UPCAST{MEAN_DEPTH}[$bi],
$UPCAST{MEDIAN_W12}[$bi],$UPCAST{MEDIAN_W34}[$bi],
$UPCAST{MEDIAN_V12}[$bi],$UPCAST{MEDIAN_V34}[$bi],
$DNCAST{MEAN_PITCH}[$bi],$DNCAST{MEAN_ROLL}[$bi],$DNCAST{MEAN_TILT}[$bi], # instrument tilt
$UPCAST{MEAN_PITCH}[$bi],$UPCAST{MEAN_ROLL}[$bi],$UPCAST{MEAN_TILT}[$bi],
$BT{N_SAMP}[$bi]>=$opt_k?$BT{MEDIAN_W}[$bi]:nan, # BT data
$BT{MAD_W}[$bi],$BT{N_SAMP}[$bi]
);
}
&antsOut('EOF'); open(STDOUT,">&2");
}
progress("\n");
}
#--------------------------------------
# write time-series output if requested
#--------------------------------------
if (@out_timeseries) {
progress("Writing time-series data to ");
@antsNewLayout = ('ensemble','elapsed','downcast',
'depth','xducer_sound_speed','pitch','gimbal_pitch','roll','tilt','heading',
'CTD_w','CTD_w_tt','LADCP_reflr_w','LADCP_reflr_w.sig',
'reflr_ocean_w');
foreach my $of (@out_timeseries) {
progress("<$of> ");
my($sub,$arg) = ($of =~ /^([^\(]+)\(([^\)]+)\)$/); # plot_sub(out_file)
if (defined($arg)) {
require "$WCALC/${sub}.pl";
&{$sub}($arg);
next;
}
$of = ">$of" unless ($of =~ /^$|^\s*\|/);
open(STDOUT,$of) || error("$of: $!\n");
undef($antsActiveHeader) unless ($ANTS_TOOLS_AVAILABLE);
for ($ens=$firstGoodEns; $ens<=$lastGoodEns; $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}],
$CTD{W_tt}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$LADCP{ENSEMBLE}[$ens]->{REFLR_W},
$LADCP{ENSEMBLE}[$ens]->{REFLR_W_STDDEV},
$reflr_oc_w);
}
&antsOut('EOF'); open(STDOUT,">&2");
}
progress("\n");
}
system("{ ./LADCP_w.PostProcess $out_basename $RUN $data_dir $plot_dir $log_dir; }&")
if (-x 'LADCP_w.PostProcess');
exit(0);