--- a/HISTORY Thu Mar 17 07:50:24 2016 -0400
+++ b/HISTORY Tue Mar 29 15:03:23 2016 -0400
@@ -1,9 +1,9 @@
======================================================================
H I S T O R Y
doc: Mon Oct 12 16:09:24 2015
- dlm: Thu Mar 17 07:50:06 2016
+ dlm: Tue Mar 29 15:02:04 2016
(c) 2015 A.M. Thurnherr
- uE-Info: 67 19 NIL 0 0 72 3 2 4 NIL ofnI
+ uE-Info: 84 66 NIL 0 0 72 3 2 4 NIL ofnI
======================================================================
V1.0:
@@ -65,4 +65,21 @@
- various plot improvements
- updated howto
- published
-
+ - V1.2beta6 [version.pl] [.hg/hgrc]
+ - changed surface-wave stat in [LADCP_w_ocean]
+ Mar 18-29, 2016: V1.2beta6
+ - [LADCP_VKE]:
+ - added -k to supply external eps
+ - re-designed QC checks (slope estimates)
+ - several other minor improvements
+ - [LADCP_wspec]:
+ - added median # of samples to output
+ - [LADCP_w_CTD]:
+ - added support for $ENV{VERB}
+ - [LADCP_w_ocean]
+ - added -r)ms residual filter with 4cm/s default cutoff
+ - fixed a couple of minor bugs (Sv correction had been disabled)
+ - replaced CTD w with CTD acceleration in wprof plots
+ Mar 29, 2016: V1.2beta6
+ - update antsMinLib to 6.6, perl-tools to 1.5 [version.pl]
+ - updated howto
--- a/LADCP_VKE Thu Mar 17 07:50:24 2016 -0400
+++ b/LADCP_VKE Tue Mar 29 15:03:23 2016 -0400
@@ -2,13 +2,16 @@
#======================================================================
# L A D C P _ V K E
# doc: Tue Oct 14 11:05:16 2014
-# dlm: Thu Mar 17 06:57:09 2016
+# dlm: Tue Mar 29 14:42:24 2016
# (c) 2012 A.M. Thurnherr
-# uE-Info: 429 49 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 70 65 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
$antsSummary = 'calculate VKE from LADCP-derived vertical-velocity profiles';
+# TOOD:
+# ! verify that p0fit.slope.sig is correct (-x scale factor)
+
# HISTORY:
# Oct 14, 2014: - created from [LADCPfs]
# Oct 15, 2014: - added parameterization output
@@ -50,6 +53,21 @@
# - removed ./ from figure label
# - update ANTSlib to 6.3
# Mar 16, 2016: - adapted to gmt5
+# Mar 17, 2016: - added version to plot
+# - added -k)e dissipation
+# - reduced eps scale
+# Mar 26, 2016: - added slope
+# Mar 27, 2016: - added slope stddev
+# - added opt_x
+# - re-designed QC checks
+# - added -z)ap
+# Mar 28, 2016: - made QC tests consistent with I08S observations based on p0fit.rms == 0.4
+# - added support for w.nsamp.avg
+# - removed -z from LDCP_wspec
+# - reduced -z default to 1, because 2 noise should affect only shotest scales
+# Mar 29, 2016: - changed default of -x to 1 & removed from usage
+# - BUG: -l/-a did not work
+# - removed p0fit.nsamp from output (is constant)
($ANTSLIB) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
($WCALC) = ($0 =~ m{^(.*)/[^/]*$});
@@ -66,32 +84,39 @@
use IPC::Open2;
#----------------------------------------------------------------------
+# Empirical constants
+#----------------------------------------------------------------------
+
+my($c) = 0.0215; # Thurnherr et al. (GRL 2015)
+
+#----------------------------------------------------------------------
# Usage
#----------------------------------------------------------------------
-my($c) = 0.0215; # Thurnherr et al. (GRL 2015)
-
-&antsUsage('a:bc:de:f:g:i:l:mno:p:r:s:tuw:',0,
+&antsUsage('a:bc:de:f:g:i:k:l:mno:p:r:s:tuw:x:z:',0,
'[poly-o)rder <n[0]> 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[150m]>',
'[-g)ap <max depth layer to fill with interpolation[40m]>]',
'[-w)indow <power-of-two input-records[32]>]',
'[shortwave -c)utoff <kz or lambda[100m]>]', # LADCP_VKE options
+ '[-z) ignore velocities derived from fewer than <N[1]> samples]',
'[o-m)it spectral correction] [spectral-tilt-correction -r)ange <max[0m]>]',
'[-l) override ADCP bin <length>] [-a) override pulse <length>]',
"[-e)ps-parameterization <scale[$c]>",
'[output -f)iles in <dir>]',
'[output -i)ndividual spectra <basename>]',
- '[output -p)lot <ps-file>',
+ '[output -p)lot <ps-file>] [-k)e dissipation <file:field>]',
'[file]');
+&antsCardOpt(\$opt_z,1); # number of w samples to require
+
#----------------------------------------------------------------------
# Calculate VKE spectra with [LADCP_wspec] if input is a w_ocean file
#----------------------------------------------------------------------
-if (defined(fnrNoErr('dc_w'))) { # pre-process with LADCP_wspec when handed vertical-velocity input
- my($opts);
+if (defined(fnrNoErr('dc_w'))) { # pre-process with LADCP_wspec when handed vertical-velocity input
+ my($opts); # pass options
$opts .= ' -d' if ($opt_d);
$opts .= ' -u' if ($opt_u);
$opts .= ' -b' if ($opt_b);
@@ -100,20 +125,25 @@
$opts .= " -g $opt_g" if defined($opt_g);
$opts .= " -w $opt_w" if defined($opt_w);
$opts .= " -o $opt_o" if defined($opt_o);
- open2(\*FROMCLD,\*TOCLD,"LADCP_wspec $opts") || # spawn sub-process
+ open2(\*FROMCLD,\*TOCLD,"LADCP_wspec $opts") || # spawn sub-process
croak("LADCP_wspec $opts: $!\n");
open(STDIN,"<&FROMCLD") || croak("dup(FROMCLD): $!\n");
close(FROMCLD);
- print(TOCLD $antsOldHeaders); undef($antsOldHeaders); # feed already gobbled header & first record to child
- print(TOCLD $antsPeekBuffer); undef($antsPeekBuffer);
- undef(%P); # shouldn't matter, because we'll get the same %PARAMs back
- undef(@antsLayout); # shouldn't matter, because it will get overwritten
- while (<>) { print(TOCLD $_); } # feed remaining data
+ print(TOCLD $antsOldHeaders); undef($antsOldHeaders); # feed already gobbled header
+ if (defined($opt_l) || defined($opt_a)) { # override bin size and/or pulse length
+ &antsAddParams('ADCP_bin_length',$opt_l) if defined($opt_l);
+ &antsAddParams('ADCP_pulse_length',$opt_a) if defined($opt_a);
+ print(TOCLD $antsCurParams);
+ }
+ print(TOCLD $antsPeekBuffer); undef($antsPeekBuffer); # feed first record
+ undef(%P); # shouldn't matter, because we'll get the same %PARAMs back
+ undef(@antsLayout); # shouldn't matter, because it will get overwritten
+ while (<>) { print(TOCLD $_); } # feed remaining data
close(TOCLD);
} elsif (defined(fnrNoErr('pwrdens.1'))) {
- croak("$0: -d, -u, -b, -w, -s meaningless when $0 used with spectral input\n")
+ croak("$0: -a, -l, -d, -u, -b, -w, -s meaningless when $0 used with spectral input\n")
if ($opt_d || $opt_u || $opt_b || defined($opt_w) ||
- defined($opt_s) || defined($opt_g));
+ defined($opt_s) || defined($opt_g)) || defined($opt_a) || defined($opt_l);
} else {
if ($ARGV[0]) {
croak("$ARGV[0]: no such file or directory\n");
@@ -126,27 +156,25 @@
# Handle LADCP_VKE usage & read data
#----------------------------------------------------------------------
-&antsAddParams('ADCP_bin_length',$opt_l) # override bin length
- if defined($opt_l);
-&antsAddParams('ADCP_pulse_length',$opt_a) # override pulse length
- if defined($opt_a);
+&antsAddParams('LADCP_VKE::wsamp.min',$opt_z); # require min # of w samples
-&antsFloatOpt(\$opt_e,$c); # default parameterization
+&antsFloatOpt(\$opt_e,$c); # default parameterization
+&antsFloatOpt(\$opt_x,1); # spectral fit stddev scale factor
-if (defined($opt_c)) { # shortwave cutoff supplied
+if (defined($opt_c)) { # shortwave cutoff supplied
$lmin = ($opt_c < 1) ? 2*$PI/$opt_c : $opt_c;
&antsAddParams('LADCP_VKE::shortwave_cutoff',2*$PI/$lmin); # ensure eps.VKE is calculated below
-} elsif (defined(antsParam('shortwave_cutoff'))) { # cutoff already applied
+} elsif (defined(antsParam('shortwave_cutoff'))) { # cutoff already applied
$lmin = 2*$PI/antsParam('shortwave_cutoff');
-} else { # use 100m default cutoff
+} else { # use 100m default cutoff
$lmin = 100;
&antsAddParams('LADCP_VKE::shortwave_cutoff',2*$PI/$lmin); # ensure eps.VKE is calculated below
}
-$lmax = 9e99; # no longwave cutoff implemented yet
+$lmax = 9e99; # no longwave cutoff implemented yet
-&antsInstallBufFull(0); # load entire file
+&antsInstallBufFull(0); # load entire file
&antsIn();
-my($Hbuf) = $antsOldHeaders; # save for later
+my($Hbuf) = $antsOldHeaders; # save for later
&antsRequireParam('profile_id');
&antsRequireParam('lambda.1');
@@ -158,7 +186,7 @@
&antsAddParams('ADCP_pulse_length',antsParam('ADCP_bin_length'))
unless defined(antsParam('ADCP_pulse_length'));
-$imin = 0; # find frequency bin limits
+$imin = 0; # find frequency bin limits
for ($nfreq=1; defined(antsParam("lambda.$nfreq")); $nfreq++) {
$imin = $nfreq if ($imin==0 && antsParam("lambda.$nfreq")<=$lmax);
$imax = $nfreq if (antsParam("lambda.$nfreq") >= $lmin);
@@ -173,6 +201,7 @@
$widxf = fnr('widx');
$mindf = fnr('depth.min');
$maxdf = fnr('depth.max');
+$wsf = fnr('w.nsamp.avg');
#----------------------------------------------------------------------
# Redirect STDOUT & create plot if STDOUT is a tty
@@ -216,7 +245,7 @@
}
-sub fit_universal_w_spec($) # vertical velocity => p0
+sub fit_universal_w_spec($) # vertical velocity => p0
{
my($r) = @_;
my($nsamp) = $fs_fmax - $fs_fmin + 1;
@@ -225,10 +254,11 @@
# fit slope-2 line in log-log space (main estimator)
#---------------------------------------------------
- if ($nsamp >= 2) { # require min 2 wavenumber samples
+ if ($nsamp >= 2) { # require min 2 wavenumber samples
- my($sumd,$sumx,$sumy) = (0,0,0);
- for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) { # loop over wavenumbers
+ my($DOF) = ($opt_d || $opt_u) ? 1 : 2; # degrees of freedom
+ my($sumd,$sumx,$sumy) = (0,0,0); # fit kz^-2 power law
+ for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) {
my($i) = $f - $pg_fmin;
$sumd += log10($ants_[$r][$f]) + 2*log10(antsRequireParam("k.$i"));
$sumx += log10(antsParam("k.$i"));
@@ -237,24 +267,68 @@
my($p0) = $sumd/$nsamp;
$ants_[$r][$p0f] = 10**$p0;
- my($avgx) = $sumx/$nsamp; # avg for r calc
+ my($avgx) = $sumx/$nsamp; # avg for r calc
my($avgy) = $sumy/$nsamp;
-
- my($sumsq,$sxx,$syy,$sxy) = (0,0,0,0); # calc rms misfit & correlation coeff
+ my($sumsqerr,$sxx,$syy,$sxy,$sumsqxt,$sumx,$sumy,$sumwt) = # r, rms error, pwrlaw slope
+ (0,0,0,0,0,0,0,0);
for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) {
- my($i) = $f - $pg_fmin;
- my($xt) = log10(antsParam("k.$i")) - $avgx;
- my($yt) = log10($ants_[$r][$f]) - $avgy;
- $sumsq += ($p0 - 2*log10(antsParam("k.$i")) - log10($ants_[$r][$f]))**2;
- $sxx += $xt * $xt; $syy += $yt * $yt; $sxy += $xt * $yt;
+ my($i) = $f - $pg_fmin;
+ my($x) = log10(&antsParam("k.$i"));
+ my($y) = log10($ants_[$r][$f]);
+ my($ysig) = $opt_x * $y / sqrt($DOF);
+ my($xt) = $x - $avgx; $sxx += &SQR($xt); # correlation coeff (r)
+ my($yt) = $y - $avgy; $syy += &SQR($yt); $sxy += $xt * $yt;
+ my($wt) = 1 / &SQR($ysig); $sumwt += $wt; # slope (linear fit in log-log space)
+ $sumx += $x * $wt;
+ $sumy += $y * $wt;
+ $sumsqerr += ($p0 - 2*$x - $y)**2; # rms error
}
- $ants_[$r][$rmsf] = sqrt($sumsq/$nsamp);
- $ants_[$r][$rf] = $sxy/(sqrt($sxx * $syy) + $SMALL_AMOUNT);
+ my($midx) = $sumx / $sumwt;
+# print(STDERR "$sumx:$sumy:$sumwt\n");
+# print(STDERR "$midx\n");
+ my($sumsqdx,$sumslp) = (0,0);
+ for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) {
+ my($i) = $f - $pg_fmin;
+ my($x) = log10(&antsParam("k.$i"));
+ my($y) = log10($ants_[$r][$f]);
+ my($ysig) = $opt_x * $y / sqrt($DOF);
+ my($dx) = ($x - $midx) / $ysig;
+ $sumsqdx += &SQR($dx);
+ $sumslp += $dx * $y / $ysig;
+ }
+ $ants_[$r][$rmsf] = sqrt($sumsqerr/$nsamp);
+ $ants_[$r][$rf] = $sxy/(sqrt($sxx * $syy) + $SMALL_AMOUNT);
+ $ants_[$r][$slpf] = $sumslp / $sumsqdx;
+ $ants_[$r][$sslpf] = sqrt(1 / $sumsqdx);
} else {
&antsInfo("WARNING: no fit --- need min 2 samples");
}
- $ants_[$r][$sampf] = $nsamp;
+}
+
+#----------------------------------------------------------------------
+# Load Dissipation Data
+#----------------------------------------------------------------------
+
+my(@eps_ms,@depth_ms); # output variables
+if (defined($opt_k)) {
+ my($file,$field) = split(':',$opt_k);
+ open(ADDF,"$file") || croak("$file: $!\n"); # open file
+ my(@afl) = &antsFileLayout(ADDF); # read layout
+ my($akf,$aef);
+ for (my($f)=0; $f<=$#afl; $f++) { # find depth & eps fields
+ $akf = $f if ($afl[$f] eq 'depth');
+ $aef = $f if ($afl[$f] eq $field);
+ }
+ croak("$file: fields 'depth' and '$field' required\n")
+ unless defined($akf) && defined($aef);
+ while (1) { # load entire profile
+ my(@ar);
+ last unless @ar = &antsFileIn(ADDF);
+ next unless numberp($ar[$akf]) && numberp($ar[$aef]);
+ push(@depth_ms,$ar[$akf]); push(@eps_ms,$ar[$aef]);
+ }
+ close(ADDF);
}
#----------------------------------------------------------------------
@@ -263,11 +337,14 @@
$fspwrf = &antsNewField('pwr.fs'); # derived fields
-$p0f = &antsNewField('p0'); # VKE parameterization
-$rmsf = &antsNewField('p0fit.rms');
-$sampf = &antsNewField('p0fit.nsamp');
-$rf = &antsNewField('p0fit.r');
-$wepsf = &antsNewField('eps.VKE');
+$p0f = &antsNewField('p0'); # VKE density
+$rmsf = &antsNewField('p0fit.rms'); # rms misfit
+$rf = &antsNewField('p0fit.r'); # correlation coefficient
+$slpf = &antsNewField('p0fit.slope'); # power-law slope
+$sslpf = &antsNewField('p0fit.slope.sig'); # power-law slope stddev
+$wepsf = &antsNewField('eps.VKE'); # epsilon from VKE
+$msepsf = &antsNewField('eps.ms') # externally supplied microstructure eps
+ if defined($opt_k);
my(@outLayout) = @antsNewLayout; # save for later
for ($f=0; $f<@outLayout; $f++) { # determine last spectral field in input
@@ -321,18 +398,54 @@
integrate_fs_power($r); # calculate total finescale power
fit_universal_w_spec($r); # fit kz^-2 spectrum
- if ($latM > 5 && # Tests: 0) not equatorial
- $ants_[$r][$rmsf] <= 0.4 && # 1) rms from kz fit < 0.4 (as in paper)
- $ants_[$r][$rf] < 0 && # 2) raw spectra are red spectrum
- $ants_[$r][$p0f] > 1e-7) { # 3) exclude low-p0 tail apparent at eps<2e-10 W/kg
- $ants_[$r][$wepsf] = ($ants_[$r][$p0f] / $opt_e)**2;
+ if (numberp($ants_[$r][$p0f])) { # update min/max depth
+ $min_depth = $ants_[$r][$mindf] if ($ants_[$r][$mindf] < $min_depth);
+ $max_depth = $ants_[$r][$maxdf] if ($ants_[$r][$maxdf] > $max_depth);
+ }
+
+ #-----------------------------------------------------------------------------------------------------
+ # QC Tests:
+ # - the following limits were independently derived
+ # p0fit.rms <= 0.4 primary filter used in Thurnherr et al. (GRL 2015)
+ # -3 <= p0fit.slope <= -1 based largely on 2016 I08S data with sufficient/insufficient range
+ # p0fit.r <= -0.5 based largely on 2016 I08S data with sufficient/insufficient range
+ # w.nsamp.avg >= 50 based on observations in many data sets with weak backscatter,
+ # including DoMORE, GO-SHIP P16S, GO-SHIP I08S
+ # - then, I plotted slope & r vs. rms and found that rms = 0.4 corresponds to, on average
+ # -3 <= p0fit.slope <= -1
+ # p0fit.r <= -0.4
+ # - in a plot of rms vs nsamp the limiting value of 0.4 is hit at 50 samples
+ #
+ # => FULL SET OF MUTUALLY CONSISTENT CRITERIA
+ #
+ # Additional Empirical Filters:
+ # - latitude > 5deg guess based on Thurnherr et al., 2015 (and some Gregg et al., 2003)
+ # - p0 > 10^-7 m^s/s^2/(rad/m) based on Fig.2 in Thurnherr et al., 2015
+ #-----------------------------------------------------------------------------------------------------
+
+ if ($latM > 5 && # 1) not 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][$wsf] >= $opt_z && # 5) on average at least 50 samples
+ $ants_[$r][$p0f] > 1e-7) { # 6) exclude low-p0 tail apparent at eps<2e-10 W/kg
+ $ants_[$r][$wepsf] = ($ants_[$r][$p0f] / $opt_e)**2;
} else {
$ants_[$r][$wepsf] = nan;
}
- if (numberp($ants_[$r][$p0f])) { # update min/max depth
- $min_depth = $ants_[$r][$mindf] if ($ants_[$r][$mindf] < $min_depth);
- $max_depth = $ants_[$r][$maxdf] if ($ants_[$r][$maxdf] > $max_depth);
+ #-------------------------------------------------
+ # average external microstructure eps, if supplied
+ #-------------------------------------------------
+
+ if (defined($opt_k)) {
+ my($sum,$n) = (0,0);
+ for (my($i)=0; $i<@eps_ms; $i++) { # linearly search all eps records
+ next unless $depth_ms[$i] >= $ants_[$r][$mindf] &&
+ $depth_ms[$i] <= $ants_[$r][$maxdf];
+ $sum += $eps_ms[$i]; $n++;
+ }
+ $ants_[$r][$msepsf] = $n ? $sum / $n : nan;
}
#---------------
@@ -436,23 +549,28 @@
'FONT_LABEL 10','MAP_LABEL_OFFSET 0.2c');
$min_depth = 400 if ($min_depth > 1e4);
$max_depth = 1500 if ($max_depth < 0);
- GMT_setR(sprintf("-R5e-13/1e-7/%d/%d",round($min_depth-250,500),round($max_depth+250,500)));
+ GMT_setR(sprintf("-R1e-11/1e-7/%d/%d",round($min_depth-250,500),round($max_depth+250,500)));
GMT_setJ(sprintf('-JX%fl/-%f',$plotsize/2.2,$plotsize/2.2));
GMT_psbasemap('-X0.3 -Y0.3 -Bg1a1f3p:"@~e@~@-VKE@- [W kg@+-1@+]":/g1000a500f100:"Depth [m]":wEsN');
GMT_psxy();
for (my($r)=0; $r<@ants_; $r++) {
- next unless numberp($ants_[$r][$wepsf]);
- my($R) = 0; my($G) = int(200*(1-$r/@ants_));
+ my($R) = 0; my($G) = int(200*(1-$r/@ants_)); # calculate color ramp
my($B) = ($r < @ants_/2) ? 150 : int(100+100*(1-$r/@ants_));
- print(GMT "> -W2,$R/$G/$B\n");
- print(GMT "$ants_[$r][$wepsf] $ants_[$r][$mindf]\n");
- print(GMT "$ants_[$r][$wepsf] $ants_[$r][$maxdf]\n");
+ if (numberp($ants_[$r][$msepsf])) { # microstructure eps
+ print(GMT "> -W1.5,DarkOrange2\n");
+ print(GMT "$ants_[$r][$msepsf] $ants_[$r][$mindf]\n");
+ print(GMT "$ants_[$r][$msepsf] $ants_[$r][$maxdf]\n");
+ }
+ if (numberp($ants_[$r][$wepsf])) {
+ print(GMT "> -W2,$R/$G/$B\n"); # plot eps.w in blue
+ print(GMT "$ants_[$r][$wepsf] $ants_[$r][$mindf]\n");
+ print(GMT "$ants_[$r][$wepsf] $ants_[$r][$maxdf]\n");
+ }
}
-
GMT_end(); # finish plot
}
-@antsNewLayout = @outLayout;
+@antsNewLayout = @outLayout; # restore layout
$antsOldHeaders = $Hbuf;
$antsHeadersPrinted = 0;
--- a/LADCP_w_CTD Thu Mar 17 07:50:24 2016 -0400
+++ b/LADCP_w_CTD Tue Mar 29 15:03:23 2016 -0400
@@ -2,9 +2,9 @@
#======================================================================
# L A D C P _ W _ C T D
# doc: Mon Nov 3 17:34:19 2014
-# dlm: Wed Mar 16 17:06:21 2016
+# dlm: Sat Mar 19 18:19:43 2016
# (c) 2014 A.M. Thurnherr
-# uE-Info: 598 40 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 60 27 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
$antsSummary = 'pre-process SBE 9plus CTD data for LADCP_w';
@@ -56,6 +56,8 @@
# - BUG: Layout error when input had no lat/lon
# - added lon to simple ASCII format
# Mar 16, 2016: - adapted to gmt5
+# Mar 19, 2016: - added support for $VERB environment var
+# - added $libSBE_quiet flag
($ANTS) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
($WCALC) = ($0 =~ m{^(.*)/[^/]*$});
@@ -82,7 +84,7 @@
&antsFloatOpt(\$opt_l,2); # defaults
&antsCardOpt(\$opt_s,6);
-&antsCardOpt(\$opt_v,0);
+&antsCardOpt(\$opt_v,$ENV{VERB}); # support VERB env variable
$CNVfile = $ARGV[0]; # open CNV file
open(F,&antsFileArg());
@@ -99,6 +101,7 @@
unless defined($rec);
if ($rec =~ /^\*/) { # SBE CNV file
+ $libSBE_quiet = 1; # suppress diagnostic messages
($nfields,$nscans,$sampint,$badval,$ftype,$lat,$lon) = # decode SBE header
SBE_parseHeader(F,0,0); # SBE field names, no time check
@@ -220,6 +223,7 @@
#----------------------------------------
my($trimmed) = 0; # trim leading nan pressures
shift(@ants_),$trimmed++
+# ,printf(STDERR "-> p=$ants_[0][$pressF] c=$ants_[0][$condF]\n")
until !@ants_ ||
numberp($ants_[0][$pressF]) &&
numberp($ants_[0][$condF]) &&
Binary file LADCP_w_howto.pdf has changed
--- a/LADCP_w_ocean Thu Mar 17 07:50:24 2016 -0400
+++ b/LADCP_w_ocean Tue Mar 29 15:03:23 2016 -0400
@@ -2,231 +2,246 @@
#======================================================================
# L A D C P _ W _ O C E A N
# doc: Fri Dec 17 18:11:13 2010
-# dlm: Thu Mar 17 05:55:25 2016
+# dlm: Tue Mar 29 07:28:54 2016
# (c) 2010 A.M. Thurnherr
-# uE-Info: 229 64 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 371 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# TODO:
-# ! plots:
-# - avoid over-plotting axis labels
-# - allow for different nsamp magnitudes
-# - add "dc" "uc" labels
-# - add software version
-# - eps.VKE limits
-# - add seabed to eps.VKE profile
-# - add seabed to LADCP_w_postproc .wprof output
-# - use instrument tilt in sidelobe editing
-# - worry about water-depth differences (disabled warning)
-# - make upcast-flag valid for yoyo casts
-# - make diagnostic output 3-beam field work for Earth coordinates
-# ? remove water-depth from BT code, which is not really used and a bit of an outlier
-# because mean and stddev are used instead of median/mad
+# ! plots:
+# - avoid over-plotting axis labels
+# - allow for different nsamp magnitudes
+# - add "dc" "uc" labels
+# - add software version
+# - add seabed to eps.VKE profile
+# - add seabed to LADCP_w_postproc .wprof output
+# ! use instrument tilt in sidelobe editing
+# ! change w_tt output to w_t output; change -a to use w_t or remove -a?
+# - worry about water-depth differences (disabled warning)
+# - make upcast-flag valid for yoyo casts
+# - make diagnostic output 3-beam field work for Earth coordinates
+# ? remove water-depth from BT code, which is not really used and a bit of an outlier
+# because mean and stddev are used instead of median/mad
$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
+# 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
+# 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)
# HISTORY END
# CTD REQUIREMENTS
@@ -314,56 +329,64 @@
# Usage
#------
+require "$WCALC/defaults.pl"; # load default/option parameters
+
$antsParseHeader = 0;
-&antsUsage('3:4a:b:c:e:g:h:i:k:lm:n:o:p:qs:t:uv:Vw:x:',0,
- '[print software -V)ersion] [-v)erbosity <level[1]>]',
- '[-q)uick (no single-ping denoising)]',
- '[require -4)-beam solutions] [apply beamvel-m)ask <file> if it exists]',
- '[valid LADCP -b)ins <bin[2],bin[*]>',
- '[-c)orrelation <min[70]>] [-t)ilt <max[10deg]> [-e)rr-vel <max[0.1m/s]>]',
- '[-h water <depth|filename>]',
- '[max LADCP time-series -g)ap <length[60s]>]',
- '[-i)nitial CTD time offset <guestimate> [-u)se as final]]',
- '[calculate -n) <lags,lags[10,100]>] [lag -w)indow <sz,sz[240s,20s]>] [lag-p)iece <CTD_elapsed_min|+[,...]>]',
- '[require top-3) lags to account for <frac[0.6]> of all]',
- '[disable time-l)ag filtering]',
- '[pressure-sensor -a)cceleration correction <residual/CTD_w_tt>]',
- '[-o)utput bin <resolution[20m]>] [-k) require <min[20]> samples]',
- '[e-x)ecute <perl-expr>]',
- '<profile-id> [run-label]');
+&antsUsage('3:4a:b:c:e:g:h:i:k:lm:n:o:p:qr:s:t:uv:Vw:x:',0,
+ "[print software -V)ersion] [-v)erbosity <level[$opt_v]>]",
+ "[-q)uick (no single-ping denoising)]",
+ "[require -4)-beam solutions] [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]>]",
+ "[-r)esidual <rms.max[$opt_r m/s]>]",
+ "[-h water <depth|filename>]",
+ "[max LADCP time-series -g)ap <length[$opt_g s]>]",
+ "[-i)nitial CTD 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 correction <residual/CTD_w_tt>]",
+ "[-o)utput bin <resolution[$opt_o m]>] [-k) require <min[$opt_k]> samples]",
+ "[e-x)ecute <perl-expr>]",
+ "<profile-id> [run-label]");
- if ($opt_V) {
- printf(STDERR "+-------------------------+\n");
- printf(STDERR "| LADCP_w Software V%s |\n",$VERSION);
- printf(STDERR "|(c) 2015- A.M. Thurnherr |\n");
- printf(STDERR "+-------------------------+\n");
- exit(0);
- }
+if ($opt_V) {
+ printf(STDERR "+-------------------------+\n");
+ printf(STDERR "| LADCP_w Software V%s |\n",$VERSION);
+ printf(STDERR "|(c) 2015- A.M. Thurnherr |\n");
+ printf(STDERR "+-------------------------+\n");
+ exit(0);
+}
- &antsUsageError() if ($opt_u && !defined($opt_i));
- &antsUsageError() unless (@ARGV==1 || @ARGV==2);
+&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
+&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
+&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';
- }
+$ID = $PROF = $STN = &antsCardArg(); # station id ($STN for compatibility)
+if (@ARGV) { # run label
+ $RUN = $ARGV[0];
+ shift;
+} else {
+ $RUN = 'profiles';
+}
#-----------------------------
# Handle Processing Parameters
#-----------------------------
-require "$WCALC/defaults.pl"; # load default/option parameters
-do "$processing_param_file"; # load processing parameters
+
+require "$WCALC/default_paths.pl"; # load default input/output paths
+croak("$processing_param_file: $@")
+ unless ($err = do "$processing_param_file"); # load processing parameters
$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 .= "-i $opt_i" if defined($opt_i);
+$processing_options .= " -i $opt_i" if defined($opt_i);
+$processing_options .= " -r $opt_r" if defined($opt_r);
+$processing_options .= ' -l' if defined($opt_l);
+$processing_options .= ' -q' if defined($opt_q);
if (defined($opt_x)) { # eval cmd-line expression to override anything
$processing_options .= " -x $opt_x";
@@ -393,7 +416,7 @@
unless numberp($length_of_timelag_windows[0]) && numberp($length_of_timelag_windows[1]);
$processing_options .= " -w $opt_w";
-croak("$0: \$out_basename undefined\n") # all plotting routines use this
+croak("$0: \$out_basename undefined\n") # plotting routines use this to label the plots
unless defined($out_basename);
&antsAddParams('processing_options',$processing_options);
&antsAddParams('out_basename',$out_basename);
@@ -544,7 +567,7 @@
}
}
}
- print(STDERR " $nee ensembles edited\n");
+ progress(" $nee ensembles edited\n");
close(BVM);
}
@@ -726,10 +749,14 @@
foreach my $of (@out_LADCP) {
progress("<$of> ");
- my($plot,$out) = ($of =~ /^([^\(]+)\(([^\)]+)\)$/); # plot_sub(out_file)
- if (defined($out)) {
- require "$WCALC/${plot}.pl";
- &{$plot}($out);
+ 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*\|/);
@@ -834,7 +861,7 @@
$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]);
@@ -976,7 +1003,7 @@
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 w_pkg calcs
+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
@@ -1017,11 +1044,13 @@
$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;
- if ($ens <= $LADCP_atbottom) {
- $dc_sumsq += &SQR($CTD{W}[$scan]); $dc_n++;
- } else {
- $uc_sumsq += &SQR($CTD{W}[$scan]); $uc_n++;
+ $LADCP{ENSEMBLE}[$ens]->{CTD_SCAN} = $scan; # calc rms package acceleration
+ 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) {
@@ -1041,9 +1070,9 @@
}
}
-my($dc_rms_wpkg) = ($dc_n > 0) ? sqrt($dc_sumsq / $dc_n) : nan; # rms package velocity
-my($uc_rms_wpkg) = ($uc_n > 0) ? sqrt($uc_sumsq / $uc_n) : nan;
-&antsAddParams('dc_rms_w_pkg',$dc_rms_wpkg,'uc_rms_w_pkg',$uc_rms_wpkg);
+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_rms_accel_pkg',$dc_rms_apkg,'uc_rms_accel_pkg',$uc_rms_apkg);
if ($nWsq > 0 && $nWsqI > 0) {
&antsAddParams('rms_w_reflr_err',sqrt($sumWsq/$nWsq),'rms_w_reflr_err_interior',sqrt($sumWsqI/$nWsqI));
@@ -1080,10 +1109,10 @@
# 2) PPI
#----------------------------------------------------------------------------
-error("$0: conflicting water-depth information provided by user\n")
- if defined($opt_h) && defined($water_depth);
-
-if (defined($opt_h)) {
+error("$0: conflicting water-depth information provided by user\n") # this can only happen if
+ if defined($opt_h) && defined($water_depth); # the user is setting $water_depth
+ # explicitly?!
+if (defined($opt_h)) { # deal with user-provided water-depth info
if (numberp($opt_h)) {
$water_depth = $opt_h;
} elsif (-f $opt_h) {
@@ -1096,10 +1125,11 @@
}
}
-die("assertion failed (water_depth = $water_depth)")
- unless (!defined($water_depth) || numberp($water_depth));
+#die("assertion failed (water_depth = $water_depth)")
+# unless (!defined($water_depth) || numberp($water_depth));
-if ($LADCP{ENSEMBLE}[$LADCP_atbottom]->{XDUCER_FACING_DOWN} && !defined($water_depth)) {
+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});
@@ -1110,35 +1140,44 @@
# 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)) {
- warning(1,"using water_depth from ADCP BT data\n");
+ 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) {
+ if ($SS_use_BT) { # water depth from BT data
&antsAddParams('water_depth_from','BT_data');
$water_depth = $water_depth_BT;
$sig_water_depth = $sig_water_depth_BT;
- } else {
+ } else { # water depth from WT data
&antsAddParams('water_depth_from','echo_amplitudes');
}
}
-if (defined($water_depth)) {
+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\n",$water_depth,$sig_water_depth);
+ 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});
+ 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)) { # 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
+if ($LADCP{ENSEMBLE}[$LADCP_atbottom]->{XDUCER_FACING_DOWN}) { # DOWNLOOKER
&antsAddParams('ADCP_orientation','DL');
if (defined($water_depth)) {
@@ -1311,9 +1350,9 @@
my($bi) = $bindepth[$bin]/$opt_o;
push(@{$DNCAST{ENSEMBLE}[$bi]},$ens);
push(@{$DNCAST{ELAPSED}[$bi]},$CTD{ELAPSED}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]);
- push(@{$DNCAST{CTD_W}[$bi]},$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]);
- push(@{$DNCAST{CTD_W_tt}[$bi]},$CTD{W_tt}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]);
- push(@{$DNCAST{BIN}[$bi]},$bin);
+# push(@{$DNCAST{CTD_W}[$bi]},$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]);
+# push(@{$DNCAST{CTD_W_tt}[$bi]},$CTD{W_tt}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]);
+# push(@{$DNCAST{BIN}[$bi]},$bin);
push(@{$DNCAST{DEPTH}[$bi]},$bindepth[$bin]);
push(@{$DNCAST{W}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin]);
push(@{$DNCAST{W12}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin]);
@@ -1364,9 +1403,9 @@
my($bi) = $bindepth[$bin]/$opt_o;
push(@{$UPCAST{ENSEMBLE}[$bi]},$ens);
push(@{$UPCAST{ELAPSED}[$bi]},$CTD{ELAPSED}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]);
- push(@{$UPCAST{CTD_W}[$bi]},$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]);
- push(@{$UPCAST{CTD_W_tt}[$bi]},$CTD{W_tt}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]);
- push(@{$UPCAST{BIN}[$bi]},$bin);
+# push(@{$UPCAST{CTD_W}[$bi]},$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]);
+# push(@{$UPCAST{CTD_W_tt}[$bi]},$CTD{W_tt}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]);
+# push(@{$UPCAST{BIN}[$bi]},$bin);
push(@{$UPCAST{DEPTH}[$bi]},$bindepth[$bin]);
push(@{$UPCAST{W}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin]);
push(@{$UPCAST{W12}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin]);
@@ -1384,12 +1423,6 @@
$UPCAST{N_SAMP}[$bi] = @{$UPCAST{W}[$bi]};
}
-&antsAddParams('min_depth',$min_depth,'max_depth',$max_depth, # plot range limits
- 'min_ens',$LADCP{ENSEMBLE}[$firstGoodEns]->{NUMBER},
- 'max_ens',$LADCP{ENSEMBLE}[$realLastGoodEns]->{NUMBER},
- 'min_elapsed',$LADCP{ENSEMBLE}[$firstGoodEns]->{CTD_ELAPSED},
- 'max_elapsed',$LADCP{ENSEMBLE}[$realLastGoodEns]->{CTD_ELAPSED});
-
#------------------------------------------------------------------------------------------------------
# remove ping-coherent residuals
# - potential error sources:
@@ -1470,10 +1503,108 @@
} # unless ($opt_q);
+#----------------------------------------------------------------------
+# remove ensembles with large rms residuals
+#----------------------------------------------------------------------
+
+if (defined($opt_r)) {
+ progress("Removing ensembles with large residuals...\n");
+
+ $nerm = 0;
+ for ($ens=$firstGoodEns; $ens<=$lastGoodEns; $ens++) {
+ next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
+
+ my($sum) = my($n) = 0; # calculate rms residual
+ my(@bindepth) = calc_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;
+ my($res) = ($ens < $LADCP_atbottom) ?
+ $LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin] - $DNCAST{MEDIAN_W}[$bi] :
+ $LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin] - $UPCAST{MEDIAN_W}[$bi];
+ $sum += &SQR($res); $n++;
+ }
+
+ if ($n == 0 || sqrt($sum/$n) > $opt_r) { # ensemble is bad
+ undef($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
+ $nerm++;
+ }
+ }
+ progress("\t$nerm ensembles removed (%d%% of total)\n",round(100*$nerm/($lastGoodEns-$firstGoodEns+1)));
+
+ 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) = calc_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]);
+ push(@{$DNCAST{W34}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin]);
+ }
+ }
+ for (my($bi)=0; $bi<=$#{$DNCAST{ENSEMBLE}}; $bi++) { # bin data into profile
+ $DNCAST{MEAN_DEPTH}[$bi] = avg(@{$DNCAST{DEPTH}[$bi]});
+ $DNCAST{MEAN_ELAPSED}[$bi] = avg(@{$DNCAST{ELAPSED}[$bi]});
+ $DNCAST{MEDIAN_W}[$bi] = median(@{$DNCAST{W}[$bi]});
+ $DNCAST{MEDIAN_W12}[$bi] = median(@{$DNCAST{W12}[$bi]});
+ $DNCAST{MEDIAN_W34}[$bi] = median(@{$DNCAST{W34}[$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) = calc_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]);
+ push(@{$UPCAST{W34}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin]);
+ }
+ }
+ for (my($bi)=0; $bi<=$#{$UPCAST{ENSEMBLE}}; $bi++) {
+ $UPCAST{MEAN_DEPTH}[$bi] = avg(@{$UPCAST{DEPTH}[$bi]});
+ $UPCAST{MEAN_ELAPSED}[$bi] = avg(@{$UPCAST{ELAPSED}[$bi]});
+ $UPCAST{MEDIAN_W}[$bi] = median(@{$UPCAST{W}[$bi]});
+ $UPCAST{MEDIAN_W12}[$bi] = median(@{$UPCAST{W12}[$bi]});
+ $UPCAST{MEDIAN_W34}[$bi] = median(@{$UPCAST{W34}[$bi]});
+ $UPCAST{MAD_W}[$bi] = mad2($UPCAST{MEDIAN_W}[$bi],@{$UPCAST{W}[$bi]});
+ $UPCAST{N_SAMP}[$bi] = @{$UPCAST{W}[$bi]};
+ }
+} # if defined($opt_r)
+
#--------------------------------------------
# Calculate and Output Bin-Averaged Resiudals
#--------------------------------------------
+&antsAddParams('min_depth',$min_depth,'max_depth',$max_depth, # plot range limits
+ 'min_ens',$LADCP{ENSEMBLE}[$firstGoodEns]->{NUMBER},
+ 'max_ens',$LADCP{ENSEMBLE}[$realLastGoodEns]->{NUMBER},
+ 'min_elapsed',$LADCP{ENSEMBLE}[$firstGoodEns]->{CTD_ELAPSED},
+ 'max_elapsed',$LADCP{ENSEMBLE}[$realLastGoodEns]->{CTD_ELAPSED});
+
if (@out_BR) {
progress("Binning residuals...");
@@ -1550,10 +1681,10 @@
foreach my $of (@out_BR) {
progress("<$of> ");
- my($plot,$out) = ($of =~ /^([^\(]+)\(([^\)]+)\)$/); # plot_sub(out_file)
- if (defined($out)) {
- require "$WCALC/${plot}.pl";
- &{$plot}($out);
+ my($sub,$arg) = ($of =~ /^([^\(]+)\(([^\)]+)\)$/); # plot_sub(out_file)
+ if (defined($arg)) {
+ require "$WCALC/${sub}.pl";
+ &{$sub}($arg);
next;
}
$of = ">$of" unless ($of =~ /^$|^\s*\|/);
@@ -1607,10 +1738,10 @@
foreach my $of (@out_wsamp) {
progress("<$of> ");
- my($plot,$out) = ($of =~ /^([^\(]+)\(([^\)]+)\)$/); # plot_sub(out_file)
- if (defined($out)) {
- require "$WCALC/${plot}.pl";
- &{$plot}($out);
+ my($sub,$arg) = ($of =~ /^([^\(]+)\(([^\)]+)\)$/); # plot_sub(out_file)
+ if (defined($arg)) {
+ require "$WCALC/${sub}.pl";
+ &{$sub}($arg);
next;
}
@@ -1704,10 +1835,10 @@
foreach my $of (@out_profile) {
progress("<$of> ");
- my($plot,$out) = ($of =~ /^([^\(]+)\(([^\)]+)\)$/); # plot_sub(out_file)
- if (defined($out)) {
- require "$WCALC/${plot}.pl";
- &{$plot}($out);
+ my($sub,$arg) = ($of =~ /^([^\(]+)\(([^\)]+)\)$/); # plot_sub(out_file)
+ if (defined($arg)) {
+ require "$WCALC/${sub}.pl";
+ &{$sub}($arg);
next;
}
@@ -1748,10 +1879,10 @@
foreach my $of (@out_timeseries) {
progress("<$of> ");
- my($plot,$out) = ($of =~ /^([^\(]+)\(([^\)]+)\)$/); # plot_sub(out_file)
- if (defined($out)) {
- require "$WCALC/${plot}.pl";
- &{$plot}($out);
+ my($sub,$arg) = ($of =~ /^([^\(]+)\(([^\)]+)\)$/); # plot_sub(out_file)
+ if (defined($arg)) {
+ require "$WCALC/${sub}.pl";
+ &{$sub}($arg);
next;
}
$of = ">$of" unless ($of =~ /^$|^\s*\|/);
@@ -1785,7 +1916,7 @@
progress("\n");
}
-system("{ ./PostProcess.sh $out_basename $RUN $data_subdir $plot_subdir $log_subdir; }&")
- if (-x 'PostProcess.sh');
+system("{ ./LADCP_w.PostProcess $out_basename $RUN $data_dir $plot_dir $log_dir; }&")
+ if (-x 'LADCP_w.PostProcess');
exit(0);
--- a/LADCP_wspec Thu Mar 17 07:50:24 2016 -0400
+++ b/LADCP_wspec Tue Mar 29 15:03:23 2016 -0400
@@ -2,28 +2,32 @@
#======================================================================
# L A D C P _ W S P E C
# doc: Thu Jun 11 12:02:49 2015
-# dlm: Tue Mar 1 21:27:41 2016
+# dlm: Mon Mar 28 13:28:39 2016
# (c) 2012 A.M. Thurnherr
-# uE-Info: 26 29 NIL 0 0 72 10 2 4 NIL ofnI
+# uE-Info: 198 0 NIL 0 0 72 10 2 4 NIL ofnI
#======================================================================
$antsSummary = 'calculate VKE window spectra from LADCP profiles';
# HISTORY:
-# Jun 11, 2015: - adapted from [binpgrams]
-# Jun 12, 2015: - renamed %PARAM prefixes
-# Jun 15, 2015: - BUG: de-meaning did not respect _gap variables
-# - added %output_depth_resolution
-# - reversed semantics of -d/-u
-# Jun 16, 2015: - reversed semantics of -t
-# - re-added pwrdens.0 to make output consistent with [binpgrams]
-# Oct 12, 2015: - require ANTSlibs V6.2 for release
-# Oct 13, 2015: - adapted to [version.pl]
+# Jun 11, 2015: - adapted from [binpgrams]
+# Jun 12, 2015: - renamed %PARAM prefixes
+# Jun 15, 2015: - BUG: de-meaning did not respect _gap variables
+# - added %output_depth_resolution
+# - reversed semantics of -d/-u
+# Jun 16, 2015: - reversed semantics of -t
+# - re-added pwrdens.0 to make output consistent with [binpgrams]
+# Oct 12, 2015: - require ANTSlibs V6.2 for release
+# Oct 13, 2015: - adapted to [version.pl]
# Jan 25, 2016: - added software version %PARAM
-# Mar 1, 2016: - made trailing message much less frequent
-# - BUG: program croaked gracelessly or entered infinite loop
-# when presented with insufficient (no valid) input
-# data
+# Mar 1, 2016: - made trailing message much less frequent
+# - BUG: program croaked gracelessly or entered infinite loop
+# when presented with insufficient (no valid) input
+# data
+# Mar 27, 2016: - added -z)ap
+# Mar 28, 2016L - removed -z
+# - renamed nsamp to nspec
+# - added w.nsamp.avg
($ANTSBIN) = (`which list` =~ m{^(.*)/[^/]*$});
($ANTSLIB) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
@@ -46,66 +50,64 @@
#----------------------------------------------------------------------
&antsUsage('bc:dg:o:s:tuw:',0,
- '[poly-o)rder <n[0]> to de-mean data; -1 to disable>] [suppress cosine-t)aper]',
- '[-d)own/-u)pcast-only] [exclude -b)ottom window]',
- '[shortwave -c)utoff <kz or lambda>]',
- '[-s)urface <layer depth to exclude[150m]>',
- '[-g)ap <max depth layer to fill with interpolation[40m]>]',
- '[-w)indow <power-of-two input-records[32]>]',
- '[LADCP-profile(s)]');
+ '[poly-o)rder <n[0]> to de-mean data; -1 to disable>] [suppress cosine-t)aper]',
+ '[-d)own/-u)pcast-only] [exclude -b)ottom window]',
+ '[shortwave -c)utoff <kz or lambda>]',
+ '[-s)urface <layer depth to exclude[150m]>',
+ '[-g)ap <max depth layer to fill with interpolation[40m]>]',
+ '[-w)indow <power-of-two input-records[32]>]',
+ '[LADCP-profile(s)]');
-&antsIntOpt(\$opt_o,0); # polynomial order to remove
-if ($opt_o >= 0) { # init model
- &modelUsage();
- matrix(\@covar,1,$modelNFit,1,$modelNFit);
- vector(\@afunc,1,$modelNFit);
- &antsAddParams('wspec::demean_poly_order',$opt_o);
+&antsIntOpt(\$opt_o,0); # polynomial order to remove
+if ($opt_o >= 0) { # init model
+ &modelUsage();
+ matrix(\@covar,1,$modelNFit,1,$modelNFit);
+ vector(\@afunc,1,$modelNFit);
+ &antsAddParams('wspec::demean_poly_order',$opt_o);
}
croak("$0: cannot ignore both down- and upcast\n")
- unless ($opt_d+$opt_u < 2);
+ unless ($opt_d+$opt_u < 2);
if ($opt_d) {
- &antsAddParams('wspec::input_data','dc');
+ &antsAddParams('wspec::input_data','dc');
} elsif ($opt_u) {
- &antsAddParams('wspec::input_data','uc');
+ &antsAddParams('wspec::input_data','uc');
} else {
- &antsAddParams('wspec::input_data','dc/uc');
+ &antsAddParams('wspec::input_data','dc/uc');
}
&antsAddParams('wspec::cos_taper_applied',$opt_t ? 'no' : 'yes');
&antsAddParams('wspec::btm_window_included',$opt_b ? 'no' : 'yes');
-if (defined($opt_c)) { # shortwave cutoff
- $kzlim = ($opt_c < 1) ? $opt_c : 2*$PI/$opt_c;
- &antsAddParams('wspec::shortwave_cutoff',$kzlim);
+if (defined($opt_c)) { # shortwave cutoff
+ $kzlim = ($opt_c < 1) ? $opt_c : 2*$PI/$opt_c;
+ &antsAddParams('wspec::shortwave_cutoff',$kzlim);
} else {
- $kzlim = 9e99;
+ $kzlim = 9e99;
}
-&antsCardOpt($opt_w); # window size
+&antsCardOpt($opt_w); # window size
-&antsCardOpt(\$opt_g,40); # gap length [m]
+&antsCardOpt(\$opt_g,40); # gap length [m]
&antsAddParams('wspec::min_gap_thickness',$opt_g);
-&antsCardOpt(\$opt_s,150); # surface layer
+&antsCardOpt(\$opt_s,150); # surface layer
&antsAddParams('wspec::surface_layer',$opt_s);
-&ISUsage; # interpolation model
+&ISUsage; # interpolation model
#----------------------------------------------------------------------
# Read Data & Define Layout
#----------------------------------------------------------------------
croak("LADCP_VKE: spectral-input mode must be selected manually (-f)\n")
- if defined(fnrNoErr('pwrdens.0')) && !defined(fnrNoErr('dc_w'));
+ if defined(fnrNoErr('pwrdens.0')) && !defined(fnrNoErr('dc_w'));
-$zfnr = fnr('depth'); # required fields
-unless (defined(fnrNoErr('dc_w'))) { #
-}
-$dcwfnr = fnr('dc_w');
-$ucwfnr = fnr('uc_w');
-$habfnr = fnrNoErr('hab'); # optional fields
+$zfnr = fnr('depth'); # required fields
+$dcwfnr = fnr('dc_w'); $dcsfnr = fnr('dc_w.nsamp');
+$ucwfnr = fnr('uc_w'); $ucsfnr = fnr('uc_w.nsamp');
+$habfnr = fnrNoErr('hab'); # optional fields
-&antsInstallBufFull(0); # read entire file
+&antsInstallBufFull(0); # read entire file
&antsIn();
while (@ants_ && $ants_[0][$zfnr] < $opt_s) { # remove surface layer
@@ -141,7 +143,7 @@
$resolution_bandwidth *= 2*$PI;
&antsAddParams('wspec::resolution_bandwidth',$resolution_bandwidth);
-push(@antsNewLayout,'widx','depth','depth.min','depth.max','hab','nsamp');
+push(@antsNewLayout,'widx','depth','depth.min','depth.max','hab','nspec','w.nsamp.avg');
for (my($i)=0; $i<$opt_w/2+1; $i++) {
my($kz) = 2*$PI*$i/$zrange;
last if ($kz > $kzlim);
@@ -193,6 +195,16 @@
return avg(@vals);
}
+sub medianF($$) # average field over window
+{
+ my($f,$r) = @_;
+ my(@vals);
+
+ push(@vals,$ants_[$r++][$f])
+ while (@vals < $opt_w);
+ return median(@vals);
+}
+
unless ($opt_t) { # compile taper function only if needed
sub cosTaperWeight($)
{
@@ -238,10 +250,11 @@
push(@out,$ants_[$fromR][$zfnr]);
push(@out,$ants_[$fromR+$opt_w-1][$zfnr]);
}
- push(@out,defined($habfnr) ?
+ push(@out,defined($habfnr) ? # hab
avgF($habfnr,$fromR) : nan);
- push(@out,nan);
- for ($i=0; $i<=$opt_w/2; $i++) {
+ push(@out,nan); # nspec
+ push(@out,nan); # w.nsamp.avg
+ for ($i=0; $i<=$opt_w/2; $i++) { # power
push(@out,nan);
}
&antsOut(@out); # output nan record and go to next window
@@ -335,8 +348,12 @@
}
push(@out,defined($habfnr) ? # hab
avgF($habfnr,$fromR) : nan);
- my($nsamp) = !$dc_gap + !$uc_gap; # nsamp
- push(@out,$nsamp);
+ my($nspec) = !$dc_gap + !$uc_gap; # nspec
+ push(@out,$nspec);
+ my($nsamp_sum) = my($nsn) = 0; # w.nsamp.avg
+ $nsamp_sum+=medianF($dcsfnr,$fromR),$nsn++ unless ($dc_gap); # median to avoid biasing by short bottle stops
+ $nsamp_sum+=medianF($ucsfnr,$fromR),$nsn++ unless ($uc_gap);
+ push(@out,$nsamp_sum/$nsn);
my($totP) = 0; # power
my($i);
@@ -344,7 +361,7 @@
my($sumP) = 0;
$sumP += $dc_pwr[$i] * $taper_correction unless ($dc_gap);
$sumP += $uc_pwr[$i] * $taper_correction unless ($uc_gap);
- push(@out,$sumP/$nsamp/$resolution_bandwidth)
+ push(@out,$sumP/$nspec/$resolution_bandwidth)
unless (antsParam("k.$i") > $kzlim);
$totP += $sumP;
}
--- a/acoustic_backscatter.pl Thu Mar 17 07:50:24 2016 -0400
+++ b/acoustic_backscatter.pl Tue Mar 29 15:03:23 2016 -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: Tue Jan 26 19:24:42 2016
+# dlm: Sat Mar 26 06:10:57 2016
# (c) 2010 A.M. Thurnherr
-# uE-Info: 171 0 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 29 92 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -25,7 +25,8 @@
# Apr 21, 2015: - added debug statements
# May 14, 2015: - BUG: code did not work for partial-depth casts
# Jun 18, 2015: - removed assertion marked by ##???, which bombed on P16N1#41 DL
-# Jan 26, 2016: - adeed %PARAMs
+# Jan 26, 2016: - added %PARAMs
+# Mar 26, 2016: - BUG: nSv was declared local to this scope even though it is used outside
#----------------------------------------------------------------------
# Volume Scattering Coefficient, following Deines (IEEE 1999)
@@ -68,7 +69,7 @@
# nSv[$depth][$bin] number of samples in bin
#----------------------------------------------------------------------
-my(@sSv,@nSv);
+# my(@sSv,@nSv); BAD: nSv is referenced in [LADCP_w_ocean]
sub calc_backscatter_profs($$)
{
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/default_paths.pl Tue Mar 29 15:03:23 2016 -0400
@@ -0,0 +1,129 @@
+#======================================================================
+# D E F A U L T _ P A T H S . P L
+# doc: Tue Mar 29 07:09:52 2016
+# dlm: Tue Mar 29 13:47:24 2016
+# (c) 2016 A.M. Thurnherr
+# uE-Info: 11 0 NIL 0 0 72 0 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+# Mar 29, 2016: - split from [defaults.pl]
+
+#======================================================================
+# ProcessingParams file selection
+#======================================================================
+
+if (-r "ProcessingParams.$RUN") {
+ $processing_param_file = "ProcessingParams.$RUN";
+} elsif (-r "ProcessingParams.default") {
+ $processing_param_file = "ProcessingParams.default";
+} elsif (-r "ProcessingParams") {
+ $processing_param_file = "ProcessingParams";
+} else {
+ error("$0: cannot find either <ProcessingParams.$RUN> or <ProcessingParams[.default]>\n");
+}
+
+#======================================================================
+# Output
+#======================================================================
+
+# The "base name" of all output files (usually 0-padded 3-digits)
+
+$out_basename = sprintf('%03d',$PROF);
+
+
+# Output subdirectories
+# these are automatically created as long as they don't contain a "/"
+
+$data_dir = $plot_dir = $log_dir = $RUN;
+unless (-d $data_dir) {
+ unless ($data_dir =~ m{/}) {
+ warning(0,"Creating data sub-directory ./$data_dir\n");
+ mkdir($data_dir);
+ }
+ error("$data_dir: no such directory\n") unless (-d $data_dir);
+}
+unless (-d $plot_dir) {
+ unless ($plot_dir =~ m{/}) {
+ warning(0,"Creating plot sub-directory ./$plot_dir\n");
+ mkdir($plot_dir);
+ }
+ error("$plot_dir: no such directory\n") unless (-d $plot_dir);
+}
+unless (-d $log_dir) {
+ unless ($log_dir =~ m{/}) {
+ warning(0,"Creating log-file sub-directory ./$log_dir\n");
+ mkdir($log_dir);
+ }
+ error("$log_dir: no such directory\n") unless (-d $log_dir);
+}
+
+
+#----------------------------------------------------------------------
+# Processing log (diagnostic messages) output
+#----------------------------------------------------------------------
+
+$out_log = "$log_dir/$out_basename.log";
+
+
+#----------------------------------------------------------------------
+# Vertical-velocity profile output and plots:
+#
+# Data:
+# *.wprof vertical velocity profiles
+#
+# Plots:
+# *_wprof.ps vertical velocity profiles (main output plot)
+#----------------------------------------------------------------------
+
+@out_profile = ("plot_wprof($plot_dir/${out_basename}_wprof.ps)",
+ "$data_dir/$out_basename.wprof");
+
+
+#----------------------------------------------------------------------
+# Vertical-velocity sample data output and plots:
+#
+# Data (in $data_dir):
+# *.wsamp w sample data
+# residuals/<prof>/<ens>.rprof OPTIONAL: per-ensemble residuals
+#
+# Plots (in $plot_dir):
+# *_wsamp.ps vertical velocity time-depth plot
+# *_residuals.ps residual vertical velocity time-depth plot
+# *_backscatter.ps volume scattering coefficient time-depth plot
+# *_correlation.ps OPTIONAL: correlation time-depth plot
+#----------------------------------------------------------------------
+
+push(@out_wsamp,"$data_dir/$out_basename.wsamp");
+#push(@out_wsamp,sprintf('dump_residual_profiles(%s/residuals/%03d)',$data_dir,$PROF));
+
+push(@out_wsamp,"plot_residuals($plot_dir/${out_basename}_residuals.ps)");
+push(@out_wsamp,"plot_backscatter($plot_dir/${out_basename}_backscatter.ps)");
+push(@out_wsamp,"plot_wsamp($plot_dir/${out_basename}_wsamp.ps)");
+#push(@out_wsamp,"plot_correlation($plot_dir/${out_basename}_correlation.ps)");
+
+
+#----------------------------------------------------------------------
+# Time-series output
+#
+# *.tis combined CTD/LADCP time-series data, including
+# package- and LADCP reference layer w
+#----------------------------------------------------------------------
+
+@out_timeseries = ("$data_dir/$out_basename.tis");
+
+#----------------------------------------------------------------------
+# Per-bin vertical-velocity residuals (plot only)
+#----------------------------------------------------------------------
+
+@out_BR = ("plot_mean_residuals($plot_dir/${out_basename}_mean_residuals.ps)");
+
+
+#----------------------------------------------------------------------
+# Time-lagging correlation statistics (plot only)
+#----------------------------------------------------------------------
+
+@out_TL = ("plot_time_lags($plot_dir/${out_basename}_time_lags.ps)");
+
+1; # return true
+
--- a/defaults.pl Thu Mar 17 07:50:24 2016 -0400
+++ b/defaults.pl Tue Mar 29 15:03:23 2016 -0400
@@ -1,9 +1,9 @@
#======================================================================
# D E F A U L T S . P L
# doc: Tue Oct 11 17:11:21 2011
-# dlm: Sun Mar 13 10:50:09 2016
+# dlm: Tue Mar 29 07:23:24 2016
# (c) 2011 A.M. Thurnherr
-# uE-Info: 130 13 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 74 39 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -68,6 +68,25 @@
# - changed outGrid_firstBin default to '*', also lastBin
# Jan 27, 2016: - added documentation
# Mar 16, 2016: - added auto creation of output directory
+# Mar 18, 2016: - added comments about -l
+# Mar 19, 2016: - improved docu
+# Mar 29, 2016: - moved out dir creation to [LADCP_w_ocean]
+# - added opt_r support
+
+#======================================================================
+# Output Log Files
+# - there are 4 verbosity levels, selected by -v
+# 0 : errors
+# 1*: UNIX-like (warnings and info messages that are not produced for every cast; *DEFAULT)
+# 2 : progress messages and useful information
+# >2: debug messges
+# - the most useful ones of these are 1 & 2
+# - verbosity level can be set with the VERB shell variable
+#======================================================================
+
+&antsCardOpt(\$opt_v,$ENV{VERB});
+$opt_v = 1 unless numberp($opt_v);
+
#======================================================================
# Data Input
@@ -103,61 +122,24 @@
$opt_b = '2,*' unless defined($opt_b);
#======================================================================
-# Logging and Output
+# Data Editing
#======================================================================
-# - there are 4 verbosity levels, selected by -v
-# 0 : errors
-# 1*: UNIX-like (warnings and info messages that are not produced for every cast; *DEFAULT)
-# 2 : progress messages and useful information
-# >2: debug messges
-# - the most useful ones of these are 1 & 2
+# The following sets the max allowable rms residual w per ensemble;
+# data from ensembles with larger rms residuals are discarded.
-&antsCardOpt(\$opt_v,$ENV{VERB});
-$opt_v = 1 unless numberp($opt_v);
-
-
-# The "base name" of all output files (usually 0-padded 3-digits)
-
-$out_basename = sprintf('%03d',$PROF);
+&antsFloatOpt(\$opt_r,0.04);
-# Output subdirectories
-
-unless (-d $RUN) {
- warning(0,"Creating output directory ./$RUN\n");
- mkdir($RUN);
- error("./$RUN: no such directory\n") unless (-d $RUN);
-}
-$data_subdir = $plot_subdir = $log_subdir = $RUN;
-
-
-# The -k option defines the minimum number of w samples required in each
-# vertical-velocity bin. The following sets the default value.
-
-&antsCardOpt(\$opt_k,20);
-
+# By default, ensembles with uncertain time-lagging are discarded.
+# 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
+# most of the time --- setting $opt_l = 1 disables the lime-lagging
+# filter for those cases.
-# The -o option sets the output grid resolution in meters. The following
-# sets the default value.
-
-&antsFloatOpt(\$opt_o,20);
-
+# $opt_l = 1;
-# The following variables limit the bins used to grid w_oean
-# - in contrast to -b, the other bins are still used e.g. for BT
-# - values recorded in %outgrid_firstbin, %outgrid_lastbin
-# - values beyond range are:
-# - greyed out in *_mean_residuals.ps
-# - not used in *_w.ps, *_residuals.ps
-
-$outGrid_firstBin = '*'; # use $LADCP_firstBin (-b)
-$outGrid_lastBin = '*'; # use $LADCP_lastBin (-b)
-
-
-#======================================================================
-# Data Editing
-#======================================================================
# The following sets the default correlation limit; measurements with
# correlations below this limit are discarded.
@@ -415,80 +397,36 @@
$BT_max_w_error = 0.03;
+
#======================================================================
-# ProcessingParams file selection
+# Gridded Velocity Profile Output
#======================================================================
-if (-r "ProcessingParams.$RUN") {
- $processing_param_file = "ProcessingParams.$RUN";
-} elsif (-r "ProcessingParams.default") {
- $processing_param_file = "ProcessingParams.default";
-} elsif (-r "ProcessingParams") {
- $processing_param_file = "ProcessingParams";
-} else {
- error("$0: cannot find either <ProcessingParams.$RUN> or <ProcessingParams[.default]>\n");
-}
-
-#----------------------------------------------------------------------
-# Processing log (diagnostic messages) output
-#----------------------------------------------------------------------
-$out_log = "$log_subdir/$out_basename.log";
-
+# The -k option defines the minimum number of w samples required in each
+# vertical-velocity bin. The following sets the default value.
-#----------------------------------------------------------------------
-# Vertical-velocity profile output and plots:
-# Data:
-# *.wprof vertical velocity profiles
-# Standard Plots:
-# *_wprof.ps vertical velocity profiles (main output plot)
-#----------------------------------------------------------------------
-
-@out_profile = ("plot_wprof($plot_subdir/${out_basename}_wprof.ps)",
- "$data_subdir/$out_basename.wprof");
+&antsCardOpt(\$opt_k,20);
-#----------------------------------------------------------------------
-# Vertical-velocity sample data output and plots:
-# Data:
-# *.wsamp w sample data
-# Standard Plots:
-# *_wsamp.ps vertical velocity time-depth plot
-# *_residuals.ps residual vertical velocity time-depth plot
-# *_backscatter.ps volume scattering coefficient time-depth plot
-# Optional Plots:
-# *_correlation.ps correlation time-depth plot
-#----------------------------------------------------------------------
+# The -o option sets the output grid resolution in meters. The following
+# sets the default value.
-push(@out_wsamp,"$data_subdir/$out_basename.wsamp");
-
-push(@out_wsamp,"plot_residuals($plot_subdir/${out_basename}_residuals.ps)");
-push(@out_wsamp,"plot_backscatter($plot_subdir/${out_basename}_backscatter.ps)");
-push(@out_wsamp,"plot_wsamp($plot_subdir/${out_basename}_wsamp.ps)");
-
-#push(@out_wsamp,"plot_correlation($plot_subdir/${out_basename}_correlation.ps)");
+&antsFloatOpt(\$opt_o,20);
-#----------------------------------------------------------------------
-# Time-series output
-# Data:
-# *.tis combined CTD/LADCP time-series data, including
-# package- and LADCP reference layer w
-#----------------------------------------------------------------------
+# The following variables limit the bins used to grid w_oean
+# - in contrast to -b, the other bins are still used e.g. for BT
+# - values recorded in %outgrid_firstbin, %outgrid_lastbin
+# - values beyond range are:
+# - greyed out in *_mean_residuals.ps
+# - not used in *_w.ps, *_residuals.ps
-@out_timeseries = ("$data_subdir/$out_basename.tis");
-
-#----------------------------------------------------------------------
-# Per-bin vertical-velocity residuals (plot only)
-#----------------------------------------------------------------------
-
-@out_BR = ("plot_mean_residuals($plot_subdir/${out_basename}_mean_residuals.ps)");
+$outGrid_firstBin = '*'; # use $LADCP_firstBin (-b)
+$outGrid_lastBin = '*'; # use $LADCP_lastBin (-b)
-#----------------------------------------------------------------------
-# Time-lagging correlation statistics (plot only)
-#----------------------------------------------------------------------
-@out_TL = ("plot_time_lags($plot_subdir/${out_basename}_time_lags.ps)");
1; # return true
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/dump_residual_profiles.pl Tue Mar 29 15:03:23 2016 -0400
@@ -0,0 +1,60 @@
+#======================================================================
+# D U M P _ R E S I D U A L _ P R O F I L E S . P L
+# doc: Thu Mar 24 07:55:07 2016
+# dlm: Tue Mar 29 13:43:56 2016
+# (c) 2016 A.M. Thurnherr
+# uE-Info: 11 30 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+# Mar 24, 2016: - created from [plot_residuals.pl]
+# Mar 29, 2016: - cleaned up
+
+sub dump_residual_profiles($)
+{
+ my($odir) = @_;
+
+ unless (-d $odir) {
+ warning(0,"Creating residual-profile output directory ./$odir\n");
+ my(@dirs) = split('/',$odir);
+ my($path) = '.';
+ foreach my $d (@dirs) {
+ $path .= "/$d";
+ mkdir($path);
+ }
+ }
+
+ return unless ($P{max_depth});
+
+ @antsNewLayout = ('bin','depth','residual');
+
+ for ($ens=$firstGoodEns; $ens<=$lastGoodEns; $ens++) {
+ next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
+
+ my($of) = sprintf('>%s/%04d.rprof',$odir,$LADCP{ENSEMBLE}[$ens]->{NUMBER});
+ open(STDOUT,$of) || error("$of: $!\n");
+ undef($antsActiveHeader) unless ($ANTS_TOOLS_AVAILABLE);
+
+ &antsAddParams('ensemble', $LADCP{ENSEMBLE}[$ens]->{NUMBER},
+ 'elapsed', $LADCP{ENSEMBLE}[$ens]->{ELAPSED},
+ 'CTD_depth', $LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH},
+ 'CTD_w', $CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
+ 'CTD_accel', $CTD{W_t}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
+ 'ADCP_tilt', $LADCP{ENSEMBLE}[$ens]->{TILT});
+
+ my(@bindepth) = calc_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;
+ my($res) = ($ens < $LADCP_atbottom) ?
+ $LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin] - $DNCAST{MEDIAN_W}[$bi] :
+ $LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin] - $UPCAST{MEDIAN_W}[$bi];
+ &antsOut($bin,$bindepth[$bin],$res);
+
+ }
+ &antsOut('EOF'); open(STDOUT,'>&2');
+ }
+}
+
+1; # return true on require
--- a/plot_wprof.pl Thu Mar 17 07:50:24 2016 -0400
+++ b/plot_wprof.pl Tue Mar 29 15:03:23 2016 -0400
@@ -1,9 +1,9 @@
#======================================================================
# P L O T _ W P R O F . P L
# doc: Sun Jul 26 11:08:50 2015
-# dlm: Thu Mar 17 06:02:11 2016
+# dlm: Thu Mar 17 12:44:45 2016
# (c) 2015 A.M. Thurnherr
-# uE-Info: 144 24 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 145 24 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -135,13 +135,15 @@
printf(GMT "0.64 1.055 bin setup\n 0.77 1.055 : %.1fm/%1.fm/%1.fm\n",
$LADCP{BLANKING_DISTANCE},$LADCP{TRANSMITTED_PULSE_LENGTH},$LADCP{BIN_LENGTH});
print(GMT "0.64 1.090 rms tilt\n 0.77 1.096 :\n");
- print(GMT "0.64 1.130 rms w\@-pkg\@-\n 0.77 1.1315 :\n");
+ print(GMT "0.64 1.130 rms a\@-pkg\@-\n 0.77 1.1315 :\n");
GMT_pstext('-F+f9,Helvetica,coral+jTL -N');
- printf(GMT "0.788 1.090 %.1f\\260\n",$P{dc_rms_tilt});
- printf(GMT "0.788 1.125 %.1fm/s\n",$P{dc_rms_w_pkg});
+# printf(GMT "0.788 1.090 %.1f\\260\n",$P{dc_rms_tilt});
+ printf(GMT "0.808 1.090 %.1f\\260\n",$P{dc_rms_tilt});
+ printf(GMT "0.78 1.125 %.1fm\@+2\@+/s\n",$P{dc_rms_accel_pkg});
GMT_pstext('-F+f9,Helvetica,SeaGreen+jTL -N');
- printf(GMT "0.89 1.090 %.1f\\260\n",$P{uc_rms_tilt});
- printf(GMT "0.89 1.125 %.1fm/s\n",$P{uc_rms_w_pkg});
+# printf(GMT "0.89 1.090 %.1f\\260\n",$P{uc_rms_tilt});
+ printf(GMT "0.91 1.090 %.1f\\260\n",$P{uc_rms_tilt});
+ printf(GMT "0.89 1.125 %.1fm\@+2\@+/s\n",$P{uc_rms_accel_pkg});
my($depth_tics) = ($plot_wprof_ymax < 1000 ) ? 'f10a100' : 'f100a500'; # AXES
setR1();
--- a/version.pl Thu Mar 17 07:50:24 2016 -0400
+++ b/version.pl Tue Mar 29 15:03:23 2016 -0400
@@ -1,9 +1,9 @@
#======================================================================
# V E R S I O N . P L
# doc: Tue Oct 13 10:40:57 2015
-# dlm: Wed Mar 16 14:07:05 2016
+# dlm: Tue Mar 29 15:02:36 2016
# (c) 2015 A.M. Thurnherr
-# uE-Info: 16 20 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 15 58 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -11,10 +11,12 @@
# Jan 4, 2016: - added $ADCP_tools_minVersion
# Mar 8, 2016: - updated antsMinLibVersion to 6.3
# Mar 16, 2016: - updated antsMinLibVersion to 6.4 (gmt5)
+# Mar 29, 2016: - updated antsMinLibVersion to 6.6 (libSBE bugs)
+# - update tools to 1.5 (obsolete getopts)
#$VERSION = '1.1'; # Jan 4, 2016
-$VERSION = '1.2beta5'; # Jan 5, 2016
+$VERSION = '1.2beta6'; # Jan 5, 2016
-$antsMinLibVersion = 6.4;
-$ADCP_tools_minVersion = 1.4;
+$antsMinLibVersion = 6.6;
+$ADCP_tools_minVersion = 1.5;