--- a/HISTORY Sat Jul 24 10:35:41 2021 -0400
+++ b/HISTORY Tue Jun 28 00:23:28 2022 -0400
@@ -1,9 +1,9 @@
======================================================================
H I S T O R Y
doc: Mon Oct 12 16:09:24 2015
- dlm: Sat Jul 24 10:35:19 2021
+ dlm: Sun Aug 8 10:35:31 2021
(c) 2015 A.M. Thurnherr
- uE-Info: 327 25 NIL 0 0 72 3 2 4 NIL ofnI
+ uE-Info: 239 67 NIL 0 0 72 3 2 4 NIL ofnI
======================================================================
----------------------------------------------------------------------
@@ -153,7 +153,7 @@
- [plot_attitude_residuals.pl]
- [plot_residual_profs.pl]
- removed assumption of 1500m/s ADCP soundspeed setting (various files)
- - added correct w12, w34 for earth velocity data
+ - added correct w12, w34 for Earth velocity data
May 19, 2016: - updated to ADCP_tools V1.6 (coord trans interface change)
@@ -236,7 +236,7 @@
Oct 3, 2017: - added -q option to LADCP_w_CTD for reprocessing as 1Hz files
-Oct 12, 2017: - re-wrote code in [LADCP_w_ocean] to deal with earth coordinates
+Oct 12, 2017: - re-wrote code in [LADCP_w_ocean] to deal with Earth coordinates
(MAJOR BUG FIX)
Oct 13, 2017: - bugfix in [edit_data.pl]
--- a/LADCP_VKE Sat Jul 24 10:35:41 2021 -0400
+++ b/LADCP_VKE Tue Jun 28 00:23:28 2022 -0400
@@ -2,9 +2,9 @@
#======================================================================
# L A D C P _ V K E
# doc: Tue Oct 14 11:05:16 2014
-# dlm: Sat Jul 24 09:44:26 2021
+# dlm: Tue May 17 16:46:00 2022
# (c) 2012 A.M. Thurnherr
-# uE-Info: 120 45 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 126 87 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
$antsSummary = 'calculate VKE from LADCP-derived vertical-velocity profiles';
@@ -118,6 +118,12 @@
# that this bug did not matter much in practice
# - changed calibration constant back, because it makes sense, and because it is now
# recorded in the meta-data
+# Oct 19, 2021: - added -h) (min samples)
+# - disabled spectral tests based on data editing implemented for GO-SHIP A20
+# - this allowed return to original $c value for VKE parameterization
+# - doubled default surface layer from 150 to 300m
+# May 17. 2022: - BUG: -h was wrong (used # of files instead of # of spectra)
+# - changed semantics to make output file name based on input file name
# HISTORY END
($ANTSLIB) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
@@ -138,26 +144,28 @@
# Empirical constants & defaults
#----------------------------------------------------------------------
-#my($c) = 0.0215; # Thurnherr et al. (GRL 2015)
-my($c) = 0.026; # increased by 21% for V1.2beta7
+my($c) = 0.0215; # Thurnherr et al. (GRL 2015)
+#my($c) = 0.026; # increased by 21% for V1.2beta7 when spectral filters were introduced; disabled Oct 2021 when filters were disabled vor V2.1
$opt_q = 3; # Equatorial band: little more than a guess based on 2015 P16N
$opt_l = 0; # [W/kg]; cutoff disabled Sep 12, 2019
$opt_a = 0; # assume background dissipation for samples that pass the tests but have eps below -l
$opt_z = 50; # number of w_ocean samples to require (note that the .wprof inputs may have harsher limits)
$opt_o = 0; # remove mean before calculating spectra
-$opt_s = 150; # surface layer to exclude from spectra
+$opt_s = 300; # surface layer to exclude from spectra
$opt_g = 40; # max gap to interpolate over
$opt_c_default = 100; # short-wavelength cutoff
+$opt_h = 2; # require 2 spectra for estimate
#----------------------------------------------------------------------
# Usage
#----------------------------------------------------------------------
$antsSuppressCommonOptions = 1;
-&antsUsage('a:bc:de:f:g:i:k:l:mno:p:q:r:s:tuw:x:yz:',0,
+&antsUsage('a:bc:de:f:g:h:i:k:l:mno:p:q:r:S:s:tuw:x:yz:',0,
"[poly-o)rder <n[$opt_o]> to de-mean data; -1 to disable>] [apply cosine-t)aper]",
'[-d)own/-u)pcast-only] [exclude -b)ottom window]', # LADCP_wspec options
- "[-s)urface <layer depth to exclude[${opt_s}m]>",
+ "[require spectral -s)amples <min[${opt_h}]>",
+ "[-S)urface <layer depth to exclude[${opt_s}m]>",
"[-g)ap <max depth layer to fill with interpolation[${opt_g}m]>]",
'[-w)indow <power-of-two input-records[32]>]',
"[shortwave -c)utoff <kz or lambda[${opt_c_default}m]>]", # LADCP_VKE options
@@ -173,6 +181,20 @@
'[output -p)lot <ps-file[#_VKE.ps]>]',
'[file...]');
+#----------------------------------------------------------------------
+# Determine Output File Name From Input File Names
+#----------------------------------------------------------------------
+
+unless ($ARGV eq '-') {
+ ($basename) = ($ARGV =~ m{([^/]*)\.[^\.]*$});
+ for (my($i)=0; $i<@ARGV; $i++) {
+ my($basename2) = ($ARGV[$i] =~ m{([^/]*)\.[^\.]*$});
+ next if ($basename2 eq $basename);
+ undef($basename);
+ last;
+ }
+}
+
#----------------------------------------------------------------------------
# Calculate VKE spectra with [LADCP_wspec] if input is a set of w_ocean files
#----------------------------------------------------------------------------
@@ -342,16 +364,13 @@
my($id) = &antsRequireParam('profile_id');
-# if ($opt_d) {
-# $opt_p = sprintf('%s/%03d_dc_VKE.ps',$opt_f,$id) unless defined($opt_p);
-# $outfile = sprintf('%s/%03d_dc.VKE',$opt_f,$id);
-# } elsif ($opt_u) {
-# $opt_p = sprintf('%s/%03d_uc_VKE.ps',$opt_f,$id) unless defined($opt_p);
-# $outfile = sprintf('%s/%03d_uc.VKE',$opt_f,$id);
-# } else {
+ if (defined($basename)) {
+ $opt_p = sprintf('%s/%s_VKE.ps',$opt_f,$basename) unless defined($opt_p);
+ $outfile = sprintf('%s/%s.VKE',$opt_f,$basename);
+ } else {
$opt_p = sprintf('%s/%03d_VKE.ps',$opt_f,$id) unless defined($opt_p);
$outfile = sprintf('%s/%03d.VKE',$opt_f,$id);
-# }
+ }
$outfile =~ s@^./@@;
open(STDOUT,">$outfile") || die("$outfile: $!\n");
} elsif (defined($opt_f)) {
@@ -388,7 +407,7 @@
$ants_[$r][$doff] += $specbuf[$bi][$wi][$doff] # update nspec once per input record
if ($f == $fs_fmax); # ... but for all files
}
- $ants_[$r][$f] = ($ns > 0) ? $sum/$ns : nan; # update averaged spectral density
+ $ants_[$r][$f] = ($ants_[$r][$doff] >= $opt_h) ? $sum/$ns : nan; # update averaged spectral density
}
}
@@ -572,6 +591,16 @@
}
#-----------------------------------------------------------------------------------------------------
+ # Except for -z and the low-latitude test, all others were disabled in October 2021 for version 2.1.
+ # This decision was based on two main reasons: 1) The checks were not used for the estimates
+ # in the Lele paper, which is the first major publiction involving the parameterization.
+ # 2) While working on the 2021 QC of the GO-SHIP A20 section, which passes through a region
+ # of extremlely weakbackscatter, I found that regional mean profiles of the filtered
+ # eps.VKE estimates agree very well with (p0/0.022)**2, but the latter are less noisy
+ # because there are signficantly more smaples. Importantly, the filtered eps.VKE were
+ # derived with c=0.026, so as an imporant side effect of disabling the filters I returned
+ # to the c value from the original publication.
+ #
# eps.VKE QC Tests:
# - the following limits were independently derived
# p0fit.rms <= 0.4 primary filter used in Thurnherr et al. (GRL 2015)
@@ -592,9 +621,9 @@
#-----------------------------------------------------------------------------------------------------
if ($latM > $opt_q && # 1) not (too) equatorial
- $ants_[$r][$rmsf] <= 0.4 && # 2) rms spectra misfit <= 0.4 (as in Thurnherr et al., GRL 2015)
- $ants_[$r][$slpf]>=-3 && $ants_[$r][$slpf]<=-1 && # 3) slope consistent with -2 power law
- $ants_[$r][$rf] <= -0.4 && # 4) p and k_z are well correlated
+# $ants_[$r][$rmsf] <= 0.4 && # 2) rms spectra misfit <= 0.4 (as in Thurnherr et al., GRL 2015)
+# $ants_[$r][$slpf]>=-3 && $ants_[$r][$slpf]<=-1 && # 3) slope consistent with -2 power law
+# $ants_[$r][$rf] <= -0.4 && # 4) p and k_z are well correlated
$ants_[$r][$wsf] >= $opt_z) { # 5) minimum # of samples
if (($ants_[$r][$p0f]/$opt_e)**2 >= $opt_l) { # level is above eps.minlim (-l)
$ants_[$r][$wepsf] = ($ants_[$r][$p0f] / $opt_e)**2; # => Thurnherr et al. (GRL 2015)
--- a/LADCP_w_CTD Sat Jul 24 10:35:41 2021 -0400
+++ b/LADCP_w_CTD Tue Jun 28 00:23:28 2022 -0400
@@ -2,9 +2,9 @@
#======================================================================
# L A D C P _ W _ C T D
# doc: Mon Nov 3 17:34:19 2014
-# dlm: Tue Jul 13 11:33:11 2021
+# dlm: Tue May 17 15:39:24 2022
# (c) 2014 A.M. Thurnherr
-# uE-Info: 226 118 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 158 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
$antsSummary = 'pre-process SBE 9plus CTD data for LADCP_w';
@@ -98,6 +98,10 @@
# Jun 30, 2021: - ditto
# Jul 13, 2021: - improved gaps PARAMs
# - added clock/transmission warnings
+# Apr 6, 2022: - added %profile_id (not tested)
+# May 10, 2022: - BUG: non-numeric ids no longer worked
+# - added -d to allow for station and cast numbering
+# May 17, 2022: - made station cast numbering work based on input file name
# HISTORY END
# NOTES:
@@ -124,7 +128,7 @@
$antsParseHeader = 0; # usage
$antsSuppressCommonOptions = 1;
$IS = &antsLoadModel('`','.nminterp','linear');
-&antsUsage("ab:c:fgi:l:morp:qs:v:w:$IS_opts",1,
+&antsUsage("ab:c:d:fgi:l:morp:qs:v:w:$IS_opts",1,
'[-v)erbosity <level[0]>]',
'[use -a)lternate sensor pair]',
'[correct -S)alinity <bias>]',
@@ -132,7 +136,7 @@
'[suppress CTD -m)odulo error correction]',
'[-s)ampling <rate[6Hz]>]',
'[lowpass w_CTD -c)utoff <limit[2s]>] [-w)inch-speed <granularity[10s]>]',
- '[profile -i)d <id>] [station -l)ocation <lat/lon>]',
+ '[profile -i)d <id>] [id -d)igits <#[3]>] [station -l)ocation <lat/lon>]',
'[-p)lot_basenames <[%03d_w_CTD.ps],[%03d_sspd.ps]>]',
'[-q)uiet (no plots)]',
'[-f)ill gaps with linear interpolation]',
@@ -145,8 +149,14 @@
&antsFloatOpt(\$opt_b,0); # salinity bias
&antsCardOpt(\$opt_v,$ENV{VERB}); # support VERB env variable
-$CNVfile = $ARGV[0]; # open CNV file
-open(F,&antsFileArg());
+$CNVfile = $ARGV[0]; # input file
+
+($basename) = ($CNVfile =~ m{([^/]*)\.[^\.]*$}); # determine number of digits to use in profile id
+$opt_d = length($basename) # - default is 3
+ if length($basename)>3 && !defined($opt_d); # - use length of input basename if it is longer than 3
+&antsCardOpt(\$opt_d,3); # - explicit -d overrides
+
+open(F,&antsFileArg()); # open CNV file
&antsAddDeps($CNVfile);
&antsActivateOut(); # activate ANTS file
@@ -406,14 +416,16 @@
$id = defined($opt_i) ? $opt_i : &antsParam('station');
_croak("$CNVfile: no station information in header => -i required\n")
unless defined($id);
-_croak("$CNVfile: non-numeric station information <$id> in header => -i required\n")
- unless numberp($id);
+#_croak("$CNVfile: non-numeric station information <$id> in header => -i required\n")
+# unless numberp($id);
+&antsAddParams('profile_id',$id);
if (-t STDOUT) {
if (numberp($id)) {
- $opt_p = '%03d_w_CTD.ps,%03d_sspd.ps'
+ my($numfmt) = "%0${opt_d}d";
+ $opt_p = "${numfmt}_w_CTD.ps,${numfmt}_sspd.ps"
unless defined($opt_p);
- $outfile = sprintf('%03d.%dHz',$id,$opt_s);
+ $outfile = sprintf("$numfmt.%dHz",$id,$opt_s);
} else {
$opt_p = '%s_w_CTD.ps,%s_sspd.ps'
unless defined($opt_p);
--- a/LADCP_w_ocean Sat Jul 24 10:35:41 2021 -0400
+++ b/LADCP_w_ocean Tue Jun 28 00:23:28 2022 -0400
@@ -2,9 +2,9 @@
#======================================================================
# 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
+# dlm: Mon Oct 18 21:27:57 2021
# (c) 2010 A.M. Thurnherr
-# uE-Info: 323 36 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 2040 92 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# TODO:
@@ -321,6 +321,13 @@
# - BUG: rms BT discrepancy was broken (residuals profile filter)
# - BUG: *_w.var did not respect -k
# Jul 23, 2021: - added %ADCP_type
+# Aug 8, 2021: - cosmetics
+# - added -p to $processing_options
+# Sep 1, 2021: - added Sv to .wprof output (dc_Sv, uc_Sv, Sv.diff)
+# - made .wprof Layout more logical
+# - added {dc,uc}_exposure_time to .wprof output
+# Oct 18, 2021: - BUG: Sv profiles included bins without valid w
+# - moved Sv profile code after additional filters;
# HISTORY END
# CTD REQUIREMENTS
@@ -341,7 +348,7 @@
# 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
+# - 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)
@@ -375,6 +382,8 @@
# - residuals are calculated with respect to down-/upcast medians
# - w12 and w34 2-beam solutions are reported without considering
# -k (min samples)
+# - elapsed times:
+# - {dc,uc}_elapsed are averages estimated BEFORE
# - equivalence, assuming -o 10:
# 1) 004DL.prof dc_w depth
# 2) bindata -Sdowncast:1 -Fw.median,depth -n 20 depth 10 004DL.w
@@ -454,7 +463,6 @@
&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
@@ -485,7 +493,7 @@
}
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 = "-c $opt_c -e $opt_e -g $opt_g -k $opt_k -o $opt_o -p $opt_p -t $opt_t -3 $opt_3";
$processing_options .= ' -l' if defined($opt_l);
if (defined($opt_q)) { # quick-and-dirty
@@ -734,7 +742,7 @@
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")
+ 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++) {
@@ -754,7 +762,7 @@
#------------------------------------------------------------
unless ($LADCP{BEAM_COORDINATES}) {
- progress("Replacing earth- with beam-velocities...\n");
+ 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]} =
@@ -770,7 +778,7 @@
#-------------------------------------------------------------------
-# Calculate earth velocities
+# 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
@@ -780,7 +788,7 @@
#-------------------------------------------------------------------
my($dummy);
-progress("Calculating earth-coordinate velocities...\n");
+progress("Calculating Earth-coordinate velocities...\n");
if ($bad_beam) {
progress("\tdiscarding velocities from beam $bad_beam\n");
&antsAddParams('bad_beam_discarded',$bad_beam);
@@ -836,13 +844,13 @@
error("$LADCP_file: insufficient valid velocities\n") unless ($nvw >= $min_valid_vels);
#----------------------------------------------
-# STEP: Edit earth-coordinate -velocity data
+# 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");
+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)
@@ -942,7 +950,7 @@
# - TILT field is set as a side-effect
#----------------------------------------------------------------------
-progress("Editing additional earth-coordinate velocity data...\n");
+progress("Editing additional Earth-coordinate velocity data...\n");
&antsAddParams('per_bin_valid_frac_lim',$per_bin_valid_frac_lim);
my($first_bad_bin);
@@ -1786,6 +1794,43 @@
} # unless ($opt_q);
+#--------------------------------
+# Add Sv to depth-binned profiles
+#--------------------------------
+
+progress("Binning Sv into depth profiles...");
+
+my(@dc_nSv);
+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 ($DNCAST{N_SAMP}[$bi]>=$opt_k) && numberp($DNCAST{MEDIAN_W}[$bi]);
+ $DNCAST{SV}[$bi] += $LADCP{ENSEMBLE}[$ens]->{SV}[$bin]; $dc_nSv[$bi]++;
+ }
+}
+my(@uc_nSv);
+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 ($UPCAST{N_SAMP}[$bi]>=$opt_k) && numberp($UPCAST{MEDIAN_W}[$bi]);
+ $UPCAST{SV}[$bi] += $LADCP{ENSEMBLE}[$ens]->{SV}[$bin]; $uc_nSv[$bi]++;
+ }
+}
+for (my($bi)=0; $bi<=max($#{$DNCAST{ENSEMBLE}},$#{$UPCAST{ENSEMBLE}}); $bi++) {
+ $DNCAST{SV}[$bi] = $dc_nSv[$bi] ? ($DNCAST{SV}[$bi]/$dc_nSv[$bi]) : nan;
+ $UPCAST{SV}[$bi] = $uc_nSv[$bi] ? ($UPCAST{SV}[$bi]/$uc_nSv[$bi]) : nan;
+}
+
+progress("\n");
+
#----------------------------------------------------------------------
# remove ensembles with large rms residuals
#----------------------------------------------------------------------
@@ -1804,7 +1849,12 @@
}
progress("\tre-binning profile data...\n");
- undef(%DNCAST); undef(%UPCAST);
+ foreach my $k (keys(%DNCAST)) { # keep Sv profiles
+ next if ($k eq 'SV');
+ undef($DNCAST{$k});
+ undef($UPCAST{$k});
+ }
+
$min_depth = 9e99; my($max_depth) = 0;
for ($ens=$firstGoodEns; $ens<$LADCP_atbottom; $ens++) { # downcast
unless (numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH})) {
@@ -1839,6 +1889,8 @@
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{MIN_ELAPSED}[$bi] = min(@{$DNCAST{ELAPSED}[$bi]});
+ $DNCAST{MAX_ELAPSED}[$bi] = max(@{$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]});
@@ -1880,6 +1932,8 @@
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{MIN_ELAPSED}[$bi] = min(@{$UPCAST{ELAPSED}[$bi]});
+ $UPCAST{MAX_ELAPSED}[$bi] = max(@{$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]});
@@ -1926,7 +1980,7 @@
$DNCAST{MEAN_RESIDUAL34}[$bi] = avg(@{$DNCAST{RESIDUAL34}[$bi]});
}
-for (my($bi)=0; $bi<$#{$DNCAST{ENSEMBLE}}; $bi++) { # rms beampair residual in 5-bin-thick layers
+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)]);
@@ -1974,12 +2028,24 @@
warning(2,"large fraction (%d%%) of profile exceeds residuals-profile threshold\n",round(100*$nbrmf))
if ($nbrmf >= 0.1);
-#------------------------
+#-----------------------------------------------
+# Apply Sv Seabed Contamination Filter
+#-----------------------------------------------
+
+progress("Applying seabed echo return filter...\n");
+&antsAddParams('seabed_contamination_filter_limit',$seabed_contamination_Sv_grad_limit);
+
+my($nbrm) = editSeabedContamination($seabed_contamination_Sv_grad_limit);
+progress("\t$nbrm depth bins removed\n");
+warning(2,"seabed contamination threshold exceeded in a large number of profile bins ($nbrm)\n")
+ if ($nbrm > 4);
+
+#----------------------------------
# Profile is now final
-# => calc <w>, w.var
-#------------------------
+# - calculate <w>, w.var metadata
+#----------------------------------
-my($dcw,$ucw,$nw) = (0,0); # <w> (vertically averaged)
+my($dcw,$ucw,$nw) = (0,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];
@@ -1994,7 +2060,7 @@
}
$ucw = $nw ? ($ucw/$nw) : nan;
-my($dcss,$ucss,$nw) = (0,0); # variance(w)
+my($dcss,$ucss,$nw) = (0,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;
@@ -2087,15 +2153,15 @@
}
}
-local($uc_bres12_max_nsamp,$dc_bres12_max_nsamp) = (0,0); # calc rms in block of well sampled bins
+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
+ $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
+ $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);
@@ -2186,9 +2252,9 @@
}
}
-#------------------------
-# output all samples (-w)
-#------------------------
+#----------------------------
+# output all samples (.wsamp)
+#----------------------------
# NB: residual fields are calculated with respect to down-/upcast medians in -o-size bins
@@ -2227,7 +2293,7 @@
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]);
+ next unless numberp($DNCAST{MEDIAN_W}[$bi]); # don't output samples contributing to bad depth bins
&antsOut(
$LADCP{ENSEMBLE}[$ens]->{NUMBER},$bin+1,
$CTD{ELAPSED}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
@@ -2308,9 +2374,9 @@
progress("\n");
}
-#----------------------------
-# Output depth-binned profile
-#----------------------------
+#-------------------------------------
+# Output depth-binned profile (.wprof)
+#-------------------------------------
if (@out_profile) {
progress("Writing vertical-velocity profiles to ");
@@ -2318,13 +2384,15 @@
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)
+ push(@antsNewLayout,'dc_lr_bp_res.rms','dc_Sv', # additional fields (inconsistent w postproc)
+ 'dc_exposure_time',
'dc_depth','dc_w12','dc_w34','dc_v12','dc_v34',
- 'uc_lr_bp_res.rms',
+ 'dc_pitch.mu','dc_roll.mu','dc_tilt.mu',
+ 'uc_lr_bp_res.rms','uc_Sv',
+ 'uc_exposure_time',
'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');
+ 'BT_w','BT_w.mad','BT_w.nsamp','Sv.diff');
foreach my $of (@out_profile) {
progress("<$of> ");
@@ -2350,17 +2418,22 @@
$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{SV}[$bi],
+ $DNCAST{MAX_ELAPSED}[$bi]-$DNCAST{MIN_ELAPSED}[$bi],
$DNCAST{MEAN_DEPTH}[$bi],
$DNCAST{MEDIAN_W12}[$bi],$DNCAST{MEDIAN_W34}[$bi],
$DNCAST{MEDIAN_V12}[$bi],$DNCAST{MEDIAN_V34}[$bi],
+ $DNCAST{MEAN_PITCH}[$bi],$DNCAST{MEAN_ROLL}[$bi],$DNCAST{MEAN_TILT}[$bi],
$UPCAST{LR_RMS_BP_RESIDUAL}[$bi], # remaining uc data
+ $UPCAST{SV}[$bi],
+ $UPCAST{MAX_ELAPSED}[$bi]-$UPCAST{MIN_ELAPSED}[$bi],
$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]
+ $BT{MAD_W}[$bi],$BT{N_SAMP}[$bi],
+ $DNCAST{SV}[$bi]-$UPCAST{SV}[$bi],
);
}
&antsOut('EOF'); open(STDOUT,">&2");
--- a/LADCP_w_postproc Sat Jul 24 10:35:41 2021 -0400
+++ b/LADCP_w_postproc Tue Jun 28 00:23:28 2022 -0400
@@ -2,9 +2,9 @@
#======================================================================
# L A D C P _ W _ P O S T P R O C
# doc: Fri Apr 24 17:15:59 2015
-# dlm: Fri Jul 23 14:03:05 2021
+# dlm: Wed May 18 08:52:28 2022
# (c) 2015 A.M. Thurnherr
-# uE-Info: 740 38 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 193 67 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
$antsSummary = 'edit and re-grid LADCP vertical-velocity samples';
@@ -81,6 +81,13 @@
# Jul 23, 2021: - added summary info to plot
# - added seabed to plot
# - added annotations to plot
+# Aug 7, 2021: - BUG: some wsamp params were wrong
+# Oct 12, 2021: - reduced window thickness for correlations from 500m to 320m
+# Oct 13, 2021: - add filter to include correlation results only if there are
+# less than 20% gaps
+# May 10, 2022: - added -f
+# May 17, 2022: - changed semantics to take output file name from input file names
+# May 18, 2022: - BUG: new semantics did not work (oops)
# HISTORY END
@@ -98,19 +105,30 @@
&antsAddParams('LADCP_w_postproc::version',$VERSION);
$antsSuppressCommonOptions = 1;
-&antsUsage('b:d:i:k:l:o:p:v:w:z',1,
+&antsUsage('b:c:d:f:i:k:l:o:p:v:w:z',1,
'[profile -i)d <id>]',
'[disable -z)eroing of <w> (disable bias correction)]',
'[-o)utput bin <resolution>] [-k) require <min> samples]',
'[-v)alid bins <DL first>,<DL last>[,<UL first>,<UL_last>]',
'[-w) <DL_dc_field>,<DL_uc_field>[,<UL_dc_field>,<UL_uc_field>]',
'[-s)urface-layer <limit[300m]>]',
- '[ouptput -d)ir <name>]',
+ '[-c)orrelation <window size[320m]>]',
+ '[ouptput -d)ir <name> -f)ile <fmt[%03d.wprof]>]',
'[-p)lot <[%03d_wprof.eps]> [-b)t <wprof>]]',
'<DL.wsamp file> [UL.wsamp file] (or only <UL.wsamp file>)');
+($basename) = ($ARGV =~ m{([^/]*)\.[^\.]*$}); # determine output file name
+
+$opt_f = '%03d.wprof' # output file name
+ unless defined($opt_f);
+
$dual_head = (@ARGV==1); # single or dual head
+if ($dual_head) {
+ my($basename2) = ($ARGV[0] =~ m{([^/]*)\.[^\.]*$});
+ undef($basename) unless ($basename2 eq $basename);
+}
+
$id = defined($opt_i) ? $opt_i : &antsParam('profile_id'); # ensure profile id exists
croak("$0: no profile_id in first file => -i required\n")
unless defined($id);
@@ -126,6 +144,7 @@
}
&antsCardOpt(\$opt_s,300); # surface layer depth limit
+&antsCardOpt(\$opc_c,320); # window thickness for correlation estimates
if (defined($opt_v)) { # bin ranges with valid data to use
($fvBin,$lvBin,$UL_fvBin,$UL_lvBin) = split(/,/,$opt_v);
@@ -170,10 +189,16 @@
#----------------------------------------------------------------------
if (-t STDOUT) {
- $opt_p = defined($opt_d) ? "$opt_d/%03d_wprof.ps" : '%03d_wprof.ps'
- unless defined($opt_p);
- $outfile = defined($opt_d) ? sprintf('%s/%03d.wprof',$opt_d,$id)
- : sprintf('%03d.wprof',$id);
+ if (defined($basename)) {
+ $outfile = defined($opt_d) ? "$opt_d/$basename.wprof" : "$basename.wprof";
+ $opt_p = defined($opt_d) ? "$opt_d/${basename}_wprof.ps" : "${basename}_wprof.ps"
+ unless defined($opt_p);
+ } else {
+ $opt_p = defined($opt_d) ? "$opt_d/%03d_wprof.ps" : '%03d_wprof.ps'
+ unless defined($opt_p);
+ $outfile = defined($opt_d) ? sprintf("%s/$opt_f",$opt_d,$id)
+ : sprintf($opt_f,$id);
+ }
open(STDOUT,">$outfile") || die("$outfile: $!\n");
}
@@ -263,6 +288,8 @@
my($ax,$ay) = (0,0);
my($ssq_diff) = 0;
+ return (nan,nan,nan,nan,nan) unless ($li > $fi); # for shallow profiles this fails
+
for (my($bi)=$fi; $bi<=$li; $bi++) {
next unless numberp($DL_dc_median[$bi]) && numberp($UL_dc_median[$bi]);
$n++;
@@ -270,7 +297,8 @@
$ay += $UL_dc_median[$bi];
$ssq_diff += ($DL_dc_median[$bi] - $UL_dc_median[$bi])**2;
}
- return (nan,nan,nan,nan,nan) unless ($n > 2);
+ return (nan,nan,nan,nan,nan) unless ($n > 0.8*($li-$fi+1));
+ die("$n,$li,$fi [$n > 0.8*($li-$fi+1)]\n") unless ($n > 0);
$ax /= $n;
$ay /= $n;
@@ -362,10 +390,10 @@
if (defined($opt_p)) { # begin summary plot
$xmin = -0.1; $x2min = -700;
$xmax = 0.35; $x2max = 500;
- $ymin = antsParam('min_depth');
+ $ymin = antsParam('depth.min');
$ymin = round($ymin-25,50);
$ymax = antsParam('water_depth');
- $ymax = antsRequireParam('max_depth') unless numberp($ymax);
+ $ymax = antsRequireParam('depth.max') unless numberp($ymax);
$ymax = round($ymax+25,50);
$plotsize = 13;
$R = "-R$xmin/$xmax/$ymin/$ymax";
@@ -535,27 +563,14 @@
'dc_wdiff.rms',$dc_rms_wdiff,'uc_wdiff.rms',$uc_rms_wdiff);
my($last_depth,$last_bi,$dc_sumsq_res,$dc_n,$uc_sumsq_res,$uc_n); # window correlation
+ my($window_size) = 320;
for (my($bi)=0; $bi<=$#dcw1; $bi++) {
($dc_R[$bi],$dc_var[$bi],$dummy,$dummy,$dc_rms_wdiff[$bi]) =
- &dc_corr(max(0,$bi-int(250/$opt_o)),min($#dcw1,$bi+int(250/$opt_o)));
+ &dc_corr(max(0,$bi-int($window_size/2/$opt_o)),min($#dcw1,$bi+int($window_size/2/$opt_o)));
($uc_R[$bi],$uc_var[$bi],$dummy,$dummy,$uc_rms_wdiff[$bi]) =
- &uc_corr(max(0,$bi-int(250/$opt_o)),min($#dcw1,$bi+int(250/$opt_o)));
+ &uc_corr(max(0,$bi-int($window_size/2/$opt_o)),min($#dcw1,$bi+int($window_size/2/$opt_o)));
}
-# my($depth) = ($bi+0.5) * $opt_o;
-# unless (defined($last_bi)) {
-# $last_bi = $bi;
-# $last_depth = $depth;
-# next;
-# }
-# if ($depth >= $last_depth+500 || $bi == $#dcw1) {
-# ($dc_R[$bi],$dc_rms[$bi],$dc_rms_wdiff[$bi]) = &dc_corr($last_bi,$bi);
-# ($uc_R[$bi],$uc_rms[$bi],$uc_rms_wdiff[$bi]) = &uc_corr($last_bi,$bi);
-# undef($last_bi);
-# next;
-# }
-# }
-
if (defined($opt_p)) { # plot 2nd-instrument profiles
GMT_psxy('-W1,coral,.');
for (my($bi)=0; $bi<=$#dcw1; $bi++) {
@@ -634,11 +649,18 @@
GMT_psxy('-W1.5,coral'); # median profiles
for (my($bi)=0; $bi<=$#dcw; $bi++) {
- printf(GMT "%f %f\n",(($dcns[$bi]>=$opt_k)?$dcwm[$bi]:nan),($bi+0.5)*$opt_o);
+# printf(GMT "%f %f\n",(($dcns[$bi]>=$opt_k)?$dcwm[$bi]:nan),($bi+0.5)*$opt_o);
+ printf(GMT "%f %f\n",(numberp($DL_dc_median[$bi]) &&
+ numberp($UL_dc_median[$bi]) &&
+ ($dcns[$bi]>=$opt_k)?$dcwm[$bi]:nan)
+ ,($bi+0.5)*$opt_o);
}
GMT_psxy('-W1.5,SeaGreen');
for (my($bi)=0; $bi<=$#ucw; $bi++) {
- printf(GMT "%f %f\n",(($ucns[$bi]>=$opt_k)?$ucwm[$bi]:nan),($bi+0.5)*$opt_o);
+ printf(GMT "%f %f\n",(numberp($DL_uc_median[$bi]) &&
+ numberp($UL_uc_median[$bi]) &&
+ ($ucns[$bi]>=$opt_k)?$ucwm[$bi]:nan)
+ ,($bi+0.5)*$opt_o);
}
GMT_psxy('-Sc0.1 -Gcoral'); # m.a.d. profiles
--- a/acoustic_backscatter.pl Sat Jul 24 10:35:41 2021 -0400
+++ b/acoustic_backscatter.pl Tue Jun 28 00:23:28 2022 -0400
@@ -1,9 +1,9 @@
#======================================================================
# A C O U S T I C _ B A C K S C A T T E R . P L
# doc: Wed Oct 20 13:02:27 2010
-# dlm: Thu Jul 1 09:37:40 2021
+# dlm: Mon Oct 18 19:28:18 2021
# (c) 2010 A.M. Thurnherr
-# uE-Info: 32 46 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 33 56 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -30,6 +30,7 @@
# May 18, 2016: - improved logging
# May 24, 2016: - calc_binDepths() -> binDepths()
# Jul 1, 2021: - made %PARAMs more standard
+# Oct 18, 2021: - made edit high-residuals not edit Sv
# HISTORY END
--- a/defaults.pl Sat Jul 24 10:35:41 2021 -0400
+++ b/defaults.pl Tue Jun 28 00:23:28 2022 -0400
@@ -1,9 +1,9 @@
#======================================================================
# D E F A U L T S . P L
# doc: Tue Oct 11 17:11:21 2011
-# dlm: Fri Jul 9 13:30:48 2021
+# dlm: Tue Oct 26 12:12:27 2021
# (c) 2011 A.M. Thurnherr
-# uE-Info: 360 33 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 183 72 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -90,6 +90,7 @@
# May 16, 2020: - updated doc
# Jun 30, 2021: - ditto
# Jul 9, 2021: - added $layer_residuals_rms_max
+# Sep 1, 2021: - added $seabed_contamination_Sv_grad_limit
# HISTORY END
#======================================================================
@@ -178,7 +179,8 @@
&antsFloatOpt(\$opt_r,'0.06,0.06');
-# By default, ensembles with uncertain time-lagging are discarded.
+# By default, ensembles with uncertain time-lagging are discarded,
+# unless the CTD time series data have been corrected for dropped scans.
# This allows profiles with dropped CTD scans to be processed without
# manual intervention. For profiles collected in very calm conditions
# (e.g. near the ice off Antarctica) time lagging is highly uncertain
@@ -328,6 +330,16 @@
$vessel_draft = 6; # in meters
+# Based on SR1b/2004 data, there are sidelobe editing does not remove
+# all seabed contamination in all profiles. This could be due to rough
+# and/or sloping seabed or instrument tilt measurement errors. The
+# contamination can easily be detected from vertical Sv gradients.
+# Based on the SR1b/2004 data, a limiting dSv/dz value of 0.1db/m
+# is suitable for WH300 instruments.
+
+$seabed_contamination_Sv_grad_limit = 0.1;
+
+
# The following function, which is called after the LADCP data have been
# read, must return the maximum horizontal reference-layer
# speed that is allowed. The following values are based on 2018 GO-SHIP
--- a/edit_data.pl Sat Jul 24 10:35:41 2021 -0400
+++ b/edit_data.pl Tue Jun 28 00:23:28 2022 -0400
@@ -1,9 +1,9 @@
#======================================================================
# E D I T _ D A T A . P L
# doc: Sat May 22 21:35:55 2010
-# dlm: Fri Jul 9 13:41:36 2021
+# dlm: Mon Oct 18 21:39:56 2021
# (c) 2010 A.M. Thurnherr
-# uE-Info: 597 63 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 55 60 NIL 0 0 72 72 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -16,7 +16,7 @@
# Dec 25, 2010: - adapted to changes in [LADCP_w]
# Aug 3, 2011: - added editTruncRange()
# Oct 10, 2011: - added editFalsePositives()
-# - BUG: when earth velocities were edited, all were
+# - BUG: when Earth velocities were edited, all were
# counted, not just those between first and lastBin
# Oct 11, 2011: - moved defaults to [defaults.pl]
# Oct 12, 2011: - added &editSurfLayer()
@@ -46,6 +46,13 @@
# Nov 17, 2018: - BUG: spurious letter "z" had crept in at some stage
# Mar 23, 2021: - updated PPI doc
# Jul 9, 2021: - added editHighResidualLayers()
+# Sep 1, 2021: - added Sv editing to editHighResidualLayers()
+# - modified sidelobe editing to include instrument tilt
+# Oct 15, 2021: - BUG: new sidelobe editing was stupid, because sidelobe
+# contamination works in the time domain and is, therefore,
+# independent of tilt
+# Oct 18, 2021: - BUG: seabed contamination was missing abs() and did not
+# work correctly with missing Sv data
# END OF HISTORY
# NOTES:
@@ -268,17 +275,20 @@
return $nrm;
}
-#======================================================================
+#===========================================================================================
# ($nvrm,$nerm) = editSideLobes($fromEns,$toEns,$water_depth)
#
# NOTES:
-# 1) When this code is executed the sound speed is known. No attempt is made to correct for
-# along-beam soundspeed variation, but the soundspeed at the transducer is accounted for.
-# 2) for surface sidelobes, water_depth == undef; surface sidelobes include the
-# vessel draft
+# - When this code is executed the sound speed is known. No attempt is made to correct for
+# along-beam soundspeed variation, but the soundspeed at the transducer is accounted for.
+# - for surface sidelobes, water_depth == undef; surface sidelobes include the
+# vessel draft
# - all velocities are counted, even those outside valid bin range,
# because the %age is not reported
-#======================================================================
+# - while this filter removes the sidelobe contamination in most profiles
+# there are still profiles with Sv.diff anomalies near the seabed
+# (based on SR1b/2004 data); for these editSeabedContamination has been implemented
+#==========================================================================================
sub editSideLobes($$$)
{
@@ -290,17 +300,20 @@
my($range) = defined($wd) ? $wd - $LADCP{ENSEMBLE}[$e]->{CTD_DEPTH}
: $LADCP{ENSEMBLE}[$e]->{CTD_DEPTH} - $vessel_draft;
$range = 0 if ($range < 0);
+
+# from UH code
my($sscorr) = $CTD{SVEL}[$LADCP{ENSEMBLE}[$e]->{CTD_SCAN}] / $LADCP{ENSEMBLE}[$e]->{SPEED_OF_SOUND};
- my($goodBins) = ($range - $sscorr*$LADCP{DISTANCE_TO_BIN1_CENTER}) * cos(rad($LADCP{BEAM_ANGLE}))
+ my($firstBadBin) = ($range - $sscorr*$LADCP{DISTANCE_TO_BIN1_CENTER}) * cos(rad($LADCP{BEAM_ANGLE}))
/ ($sscorr*$LADCP{BIN_LENGTH})
- 1.5;
my($dirty) = 0;
- for (my($bin)=int($goodBins); $bin<$LADCP{N_BINS}; $bin++) { # NB: 2 good bins implies that bin 2 is bad
+ for (my($bin)=int($firstBadBin); $bin<$LADCP{N_BINS}; $bin++) {
next unless ($bin>=0 && defined($LADCP{ENSEMBLE}[$e]->{W}[$bin]));
$dirty = 1;
$nvrm++;
undef($LADCP{ENSEMBLE}[$e]->{W}[$bin]);
+ debugmsg("sidelobe at range=$range firstBadBin=$firstBadBin ens=$e bin=$bin at CTD depth = $LADCP{ENSEMBLE}[$e]->{CTD_DEPTH}\n");
}
$nerm += $dirty;
@@ -573,7 +586,9 @@
#======================================================================
# $nbrm = editHighResidualLayers($max_lr_res)
-# - filter applied after gridding
+# - filter applied after depth binning
+# - while filter only removes values from profiles, the corresponding
+# samples are not output to .wsamp, because MEDIAN_W is not defined
# - filter based on observation that profiles of beam-pair residuals
# are good indicators of bad data, but very noisy
# - current version uses estimates in 5-bin-thick layers (200m by
@@ -590,16 +605,54 @@
for (my($bi)=0; $bi<=$#{$DNCAST{LR_RMS_BP_RESIDUAL}}; $bi++) {
next unless ($DNCAST{LR_RMS_BP_RESIDUAL}[$bi] > $limit);
$DNCAST{MEDIAN_W}[$bi] = $DNCAST{MEDIAN_W12}[$bi] = $DNCAST{MEDIAN_W34}[$bi] = nan;
+# $DNCAST{SV}[$bi] = nan;
$nbrm++;
}
for (my($bi)=0; $bi<=$#{$UPCAST{LR_RMS_BP_RESIDUAL}}; $bi++) {
next unless ($UPCAST{LR_RMS_BP_RESIDUAL}[$bi] > $limit);
$UPCAST{MEDIAN_W}[$bi] = $UPCAST{MEDIAN_W12}[$bi] = $UPCAST{MEDIAN_W34}[$bi] = nan;
+# $UPCAST{SV}[$bi] = nan;
$nbrm++;
}
return $nbrm;
}
#======================================================================
+# $nbrm = editSeabedContamination($max_lr_res)
+# - filter applied after depth binning
+# - while filter only removes values from profiles, the corresponding
+# samples are not output to .wsamp, because MEDIAN_W is not defined
+# - filter based on SR1b/2004 observation that some profiles of Sv
+# show clear anomalies near the seabed
+# - anomalies are easily detected in dSv/dz plots, with
+# $seabed_contamination_Sv_grad_limit = 0.1 being a suitable limiting
+# value for WH300 ADCPs
+#======================================================================
+
+sub editSeabedContamination($)
+{
+ my($limit) = @_;
+ my($nbrm) = 0;
+
+ for (my($bi)=$#{$DNCAST{SV}}; $bi>0; $bi--) {
+ next unless numberp($DNCAST{SV}[$bi]) && numberp($DNCAST{SV}[$bi-1]);
+ last if (abs($DNCAST{SV}[$bi]-$DNCAST{SV}[$bi-1])/$opt_o < $limit);
+ $DNCAST{MEDIAN_W}[$bi] = $DNCAST{MEDIAN_W12}[$bi] = $DNCAST{MEDIAN_W34}[$bi] = nan;
+ $DNCAST{SV}[$bi] = nan;
+ $nbrm++;
+ }
+
+ for (my($bi)=$#{$UPCAST{SV}}; $bi>0; $bi--) {
+ next unless numberp($UPCAST{SV}[$bi]) && numberp($UPCAST{SV}[$bi-1]);
+ last if (abs($UPCAST{SV}[$bi]-$UPCAST{SV}[$bi-1])/$opt_o < $limit);
+ $UPCAST{MEDIAN_W}[$bi] = $UPCAST{MEDIAN_W12}[$bi] = $UPCAST{MEDIAN_W34}[$bi] = nan;
+ $UPCAST{SV}[$bi] = nan;
+ $nbrm++;
+ }
+
+ return $nbrm;
+}
+
+#======================================================================
1;
--- a/time_lag.pl Sat Jul 24 10:35:41 2021 -0400
+++ b/time_lag.pl Tue Jun 28 00:23:28 2022 -0400
@@ -1,9 +1,9 @@
#======================================================================
# T I M E _ L A G . P L
# doc: Fri Dec 17 21:59:07 2010
-# dlm: Thu Jul 1 09:39:53 2021
+# dlm: Sun Aug 8 11:06:48 2021
# (c) 2010 A.M. Thurnherr
-# uE-Info: 81 46 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 82 60 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -79,6 +79,8 @@
# Oct 4, 2018: - added timelagging debug code
# Oct 16, 2018: - removed debug code
# Jul 1, 2021: - made %PARAMs more standard
+# Aug 8, 2021: - BUG: empty upcast made time-lagging bomb
+# HISTORY END
# DIFFICULT STATIONS:
# NBP0901#131 this requires the search-radius doubling heuristic
@@ -205,6 +207,12 @@
$n_valid_windows++;
$nBest{$so}++; $madBest{$so} += $mad;
}
+
+ unless ($n_valid_windows) {
+ $failed = 1;
+ goto CONTINUE;
+ }
+
my($maxN) = 0;
foreach my $i (keys(%nBest)) {
$maxN = $nBest{$i} if ($nBest{$i} > $maxN);
@@ -217,14 +225,6 @@
$nBest{$lag} = 0; $madBest{$lag} = 9e99;
}
-## my($med_mad) = median(values(%madBest)); # remove lags with large mads
-## my($mad_mad) = mad2($med_mad,values(%madBest));
-## foreach my $lag (keys(%nBest)) {
-## next if ($madBest{$lag} <= $med_mad+$mad_mad);
-## $n_valid_windows -= $nBest{$lag};
-## $nBest{$lag} = 0;
-## }
-
my($min_mad) = min(values(%madBest)); # remove lags with large mads
foreach my $lag (keys(%nBest)) {
next if ($madBest{$lag} <= 3*$min_mad);
@@ -302,6 +302,8 @@
# Here, either $failed is set, or we have a valid lag.
#----------------------------------------------------
+CONTINUE:
+
# if ($failed) {
# for (my($wi)=0; $wi<$n_windows; $wi++) {
# print(STDERR "$wi $so[$wi] $mad[$wi]\n");
--- a/version.pl Sat Jul 24 10:35:41 2021 -0400
+++ b/version.pl Tue Jun 28 00:23:28 2022 -0400
@@ -1,9 +1,9 @@
#======================================================================
# V E R S I O N . P L
# doc: Tue Oct 13 10:40:57 2015
-# dlm: Sat Jul 24 09:55:02 2021
+# dlm: Wed Sep 1 13:32:23 2021
# (c) 2015 A.M. Thurnherr
-# uE-Info: 32 68 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 34 68 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -30,12 +30,16 @@
# Jun 29, 2021: - updated ANTSlib to V7.2
# Jul 1, 2021: - updated ANTSlib to V7.3
# Jul 24, 2021: - updated to V2.0 (major improvements) for release
+# Aug 6, 2021: - ARGH! never actually updated version :(
+# Sep 1, 2021: - updated to V2.1 (Sv profiles & better sidelobes)
#$VERSION = '1.1'; # Jan 4, 2016
#$VERSION = '1.2'; # May 12, 2016
#$VERSION = '1.3'; # Mar 15, 2017
#$VERSION = '1.4'; # Nov 28, 2017
-$VERSION = '1.5'; # Sep 12, 2018
+#$VERSION = '1.5'; # Sep 12, 2018
+#$VERSION = '2.0'; # Aug 6, 2021
+$VERSION = '2.1'; # Sep 1, 2021
$antsMinLibVersion = 7.3;
$ADCP_tools_minVersion = 2.4;