--- a/HISTORY Wed May 25 12:14:29 2016 -0400
+++ b/HISTORY Fri Aug 05 11:02:51 2016 -0400
@@ -1,9 +1,9 @@
======================================================================
H I S T O R Y
doc: Mon Oct 12 16:09:24 2015
- dlm: Wed May 25 12:14:13 2016
+ dlm: Wed Jun 8 22:07:57 2016
(c) 2015 A.M. Thurnherr
- uE-Info: 172 19 NIL 0 0 72 3 2 4 NIL ofnI
+ uE-Info: 182 43 NIL 0 0 72 3 2 4 NIL ofnI
======================================================================
----------------------------------------------------------------------
@@ -170,3 +170,13 @@
May 25, 2016:
- updated version to V1.3beta1 [version.pl] [.hg/hgrc]
- published
+
+ May 26, 2016: V1.3beta2
+ - updated version to V1.3beta2 [version.pl] [.hg/hgrc]
+ - udated ANTS_tools lib to V1.7 (beam interpolation)
+ - adapted [LADCP_w_ocean] to beam interpolation
+
+...
+
+ Jun 8, 2016:
+ - removed plot_attitude_biases_w.pl
--- a/LADCP_w_CTD Wed May 25 12:14:29 2016 -0400
+++ b/LADCP_w_CTD Fri Aug 05 11:02:51 2016 -0400
@@ -2,9 +2,9 @@
#======================================================================
# L A D C P _ W _ C T D
# doc: Mon Nov 3 17:34:19 2014
-# dlm: Thu Mar 31 05:52:01 2016
+# dlm: Thu May 26 20:42:15 2016
# (c) 2014 A.M. Thurnherr
-# uE-Info: 62 42 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 639 69 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
$antsSummary = 'pre-process SBE 9plus CTD data for LADCP_w';
@@ -60,6 +60,8 @@
# - added $libSBE_quiet flag
# Mar 30, 2016: - added station number of ASCII format
# Mar 31, 2016: - changed version %PARAM
+# May 26, 2016: - renamed w[.raw] to CTD_w[.raw]
+# - added winch_w
($ANTS) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
($WCALC) = ($0 =~ m{^(.*)/[^/]*$});
@@ -68,24 +70,26 @@
require "$WCALC/version.pl";
require "$ANTS/ants.pl";
require "$ANTS/fft.pl";
+require "$ANTS/libstats.pl";
require "$ANTS/libGMT.pl";
require "$ANTS/libSBE.pl";
require "$ANTS/libEOS83.pl";
&antsAddParams('LADCP_w_CTD::version',$VERSION);
$antsParseHeader = 0; # usage
-&antsUsage('ai:l:orp:s:v:',1,
+&antsUsage('ai:l:orp:s:v:w:',1,
'[-v)erbosity <level[0]>]',
'[use -a)lternate sensor pair]',
'[-r)etain all data (no editing)] [allow infinite -o)utliers]',
'[-s)ampling <rate[6Hz]>]',
- '[-l)owpass w <cutoff[2s]>]',
+ '[-l)owpass w_CTD <cutoff[2s]>] [-w)inch-speed <granularity[10s]>]',
'[profile -i)d <id>]',
'[-p)lot_basenames <[%03d_w_CTD.ps],[%03d_sspd.ps]>]',
'<SBE CNV file>');
-&antsFloatOpt(\$opt_l,2); # defaults
-&antsCardOpt(\$opt_s,6);
+&antsFloatOpt(\$opt_l,2); # default low-pass cutoff for w_CTD
+&antsCardOpt(\$opt_s,6); # default output samplint rate (Hz)
+&antsFloatOpt(\$opt_w,10); # winch velocity granularity
&antsCardOpt(\$opt_v,$ENV{VERB}); # support VERB env variable
$CNVfile = $ARGV[0]; # open CNV file
@@ -579,7 +583,36 @@
@w_lp = @w;
}
-#----------------------------------------------------------------------
+#----------------------------------------
+# Estimate winch speed
+#----------------------------------------
+
+print(STDERR "Estimating winch velocity...") if ($opt_v);
+&antsAddParams('winch_velocity_granularity',$opt_w);
+
+my($from_r) = 0; my($to_r); # step 1: bin average in time
+for (my($from_r)=my($to_r)=0; $from_r<@elapsed; $from_r=$to_r+1) {
+ my($sumw) = $w_lp[2*$from_r]/@w_lp; my($n) = 1;
+ for ($to_r=$from_r+1; $to_r<@elapsed && $elapsed[$to_r]-$elapsed[$from_r]<$opt_w; $to_r++) {
+ $sumw += $w_lp[2*$to_r]/@w_lp; $n++;
+ }
+ $winch[$from_r] = $sumw/$n;
+}
+
+my($pwinch) = $winch[0];
+for (my($to_r),my($from_r)=0; $from_r<@elapsed; ) { # step 2: fill after median filtering
+ for ($to_r=$from_r+1; $to_r<@elapsed && !defined($winch[$to_r]); $to_r++) {}
+ my($nwinch) = $to_r<@elapsed ? $winch[$to_r] : $winch[$from_r];
+ my($winch) = median($pwinch,$winch[$from_r],$nwinch);
+ $pwinch = $winch[$from_r]; $winch[$from_r] = $winch;
+ while (++$from_r < $to_r) { $winch[$from_r] = $winch[$from_r-1]; }
+}
+
+printf(STDERR "\n") if ($opt_v);
+
+#----------------------------------------
+# Plot Data
+#----------------------------------------
if (defined($opt_p)) {
print(STDERR "Plotting data...\n") if ($opt_v);
@@ -597,10 +630,15 @@
printf(GMT "%f %f\n",$elapsed[$r]/60,$w_lp[2*$r]/@w_lp);
GMT_psxy('-W1,SeaGreen') if ($r == $atBtm);
}
+ GMT_psxy('-W1,magenta');
+ for ($r=0; $r<@w; $r++) {
+ printf(GMT "%f %f\n",$elapsed[$r]/60,$winch[$r]);
+ }
GMT_psbasemap('-Bg60a30f5:"Elapsed Time [min]":/g1a1f0.1:"Vertical Package Velocity [ms@+-1@+]":WeSn');
GMT_unitcoords();
- GMT_pstext('-F+f14,Helvetica,coral+jTR'); print(GMT "0.98 0.98 dc\n");
- GMT_pstext('-F+f14,Helvetica,SeaGreen+jTR'); print(GMT "0.98 0.94 uc\n");
+ GMT_pstext('-F+f14,Helvetica,coral+jBR'); print(GMT "0.98 0.96 dc\n");
+ GMT_pstext('-F+f14,Helvetica,SeaGreen+jBR'); print(GMT "0.98 0.93 uc\n");
+ GMT_pstext('-F+f14,Helvetica,magenta+jBR'); print(GMT "0.98 0.89 winch\n");
GMT_pstext('-F+f14,Helvetica,blue+jTL -N');
if (defined($outfile)) { printf(GMT "0.01 1.06 $outfile\n",$id); }
else { printf(GMT "0.01 1.06 %03d\n",$id); }
@@ -618,8 +656,8 @@
GMT_psxy('-W1.5,SeaGreen') if ($r == $atBtm);
}
GMT_unitcoords();
- GMT_pstext('-F+f14,Helvetica,coral+jTR'); print(GMT "0.98 0.02 dc\n");
- GMT_pstext('-F+f14,Helvetica,SeaGreen+jTR'); print(GMT "0.98 0.06 uc\n");
+ GMT_pstext('-F+f14,Helvetica,coral+jTR'); print(GMT "0.98 0.02 dc\n");
+ GMT_pstext('-F+f14,Helvetica,SeaGreen+jTR'); print(GMT "0.98 0.06 uc\n");
GMT_pstext('-F+f14,Helvetica,blue+jTL -N');
if (defined($outfile)) { printf(GMT "0.01 -0.06 $outfile\n",$id); }
else { printf(GMT "0.01 -0.06 %03d\n",$id); }
@@ -630,10 +668,11 @@
print(STDERR "Writing output...\n") if ($opt_v);
-@antsNewLayout = ('elapsed','press','temp','cond','depth','salin','sspd','w.raw','w','lat','lon');
+@antsNewLayout = ('elapsed','press','temp','cond','depth','salin','sspd','w_CTD.raw','w_CTD','w_winch','lat','lon');
for ($r=0; $r<@w; $r++) {
&antsOut($elapsed[$r],$press[$r],$temp[$r],$cond[$r],$depth[$r],$salin[$r],
- $sspd[$r],$w[$r],$w_lp[2*$r]/@w_lp,$lat[$r],$lon[$r]);
+ $sspd[$r],$w[$r],$w_lp[2*$r]/@w_lp,$winch[$r],
+ $lat[$r],$lon[$r]);
}
exit(0); # don't flush @ants_
--- a/LADCP_w_ocean Wed May 25 12:14:29 2016 -0400
+++ b/LADCP_w_ocean Fri Aug 05 11:02:51 2016 -0400
@@ -2,19 +2,22 @@
#======================================================================
# L A D C P _ W _ O C E A N
# doc: Fri Dec 17 18:11:13 2010
-# dlm: Tue May 24 16:33:29 2016
+# dlm: Sat Jun 11 18:34:12 2016
# (c) 2010 A.M. Thurnherr
-# uE-Info: 1868 26 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 272 65 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# TODO:
-# ! plots:
+# ! use instrument tilt in sidelobe editing
+# ! use instrument tilt in e.g. Sv
+#
+# - plots:
# - avoid over-plotting axis labels
# - allow for different nsamp magnitudes
# - add "dc" "uc" labels
# - add seabed to eps.VKE profile
# - add seabed to LADCP_w_postproc .wprof output
-# ! 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
@@ -251,13 +254,30 @@
# - added w12, w34 for earth-coordinates
# May 19, 2016: - updated to ADCP_tools V1.6
# May 24, 2016: - calc_binDepths() -> binDepths()
+# May 26, 2016: - added -d)isable bin interpolation
+# - adapted to w_CTD, w_winch
+# - re-enabled disabled L0 warning
+# - suppressing L0 warnings on screen
+# May 27, 2016: - cosmetics
+# - added reflr_w to .wsamp output
+# May 29, 2016: - BUG: do ProcessingParams did not handle all errors correctly
+# Jun 1, 2016: - made empirical accoustic backscatter correction optional
+# - split [default_paths.pl] to same plus [default_output.pl]
+# Jun 2, 2016: - added applyTiltCorrection()
+# - BUG: PPI-editing was no longer enabled automatically
+# for 150kHz instruments
+# Jun 3, 2016: - BUG: earth-coordinate data did not have gimbal-pitch
+# calculated
+# Jun 6, 2016: - removed applyTiltCorrection()
+# Jun 11, 2016: - BUG: w12, w34 were wrong for Earth-coord data
# HISTORY END
# CTD REQUIREMENTS
# - elapsed elapsed seconds; see note below
# - depth
# - ss sound speed
-# - w ddepth/dt
+# - w[_CTD] ddepth/dt
+# - w_winch OPTIONAL; winch-speed estimate
# - temp OPTIONAL; used for backscatter calculation
# - %lat/%lon OPTIONAL
@@ -302,9 +322,9 @@
# VELOCITY AMBIGUITY ERRORS
# - quite extensive tests with DIMES US2 station 146, which has a lot of
-# ambiguity velocity errors, reveal that -m 1 catches those errors
+# ambiguity velocity errors, reveal that $w_max_lim catches those errors
# quite nicely
-# - even when the errors are not filtered with -m 1, they do not
+# - even when the errors are not filtered with $w_max_lim, they do not
# affect the w profiles, as long as the median bin values are used
($ANTS) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
@@ -344,13 +364,13 @@
require "$WCALC/defaults.pl"; # load default/option parameters
$antsParseHeader = 0;
-&antsUsage('3:4a:b:c:e:g:h:i:k:lm:n:o:p:qr:s:t:uv:Vw:x:',0,
+&antsUsage('3:4a:b:c:de:g:h:i:k:lm:n:o:p:qr:s:t:uv:Vw:x:',0,
"[print software -V)ersion] [-v)erbosity <level[$opt_v]>]",
"[-q)uick (no single-ping denoising)]",
- "[require -4)-beam solutions] [apply beamvel-m)ask <file> if it exists]",
+ "[require -4)-beam solutions] [-d)isable bin interpolation] [apply beamvel-m)ask <file> if it exists]",
"[valid LADCP -b)ins <bin,bin[$opt_b]>",
"[-c)orrelation <min[$opt_c counts]>] [-t)ilt <max[$opt_t deg]> [-e)rr-vel <max[$opt_e m/s]>]",
- "[-r)esidual <rms.max[$opt_r m/s]>]",
+ "[-r)esidual <rms.max[,delta.max][$opt_r m/s]>]",
"[-h water <depth|filename>]",
"[max LADCP time-series -g)ap <length[$opt_g s]>]",
"[-i)nitial CTD time offset <guestimate> [-u)se as final]]",
@@ -386,19 +406,30 @@
$RUN = 'profiles';
}
-#-----------------------------
+#----------------------------------------------------------------
# Handle Processing Parameters
-#-----------------------------
+# - paths need to be read first to define ProcessingParams file
+# and output directories
+# - processing params file is read next to allow setting
+# plotting_level
+# - default_output.pl is read finally to set any missing output
+# variables to default values
+#----------------------------------------------------------------
-require "$WCALC/default_paths.pl"; # load default input/output paths
-croak("$processing_param_file: $@")
- unless ($err = do "$processing_param_file"); # load processing parameters
+require "$WCALC/default_paths.pl"; # define default input/output paths
+my($retval) = do $processing_param_file; # load processing parameters
+if (!defined($retval) && $@ ne '') {
+ croak("$processing_param_file: $@\n");
+} elsif (!defined($retval) && $! != 0) {
+ croak("$processing_param_file: $!\n");
+}
+require "$WCALC/default_output.pl"; # set default output plots and files
$processing_options = "-k $opt_k -o $opt_o -c $opt_c -t $opt_t -e $opt_e -g $opt_g -3 $opt_3";
$processing_options .= " -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);
+$processing_options .= ' -d' if defined($opt_d);
if (defined($opt_x)) { # eval cmd-line expression to override anything
$processing_options .= " -x $opt_x";
@@ -410,6 +441,12 @@
$RDI_Coords::minValidVels = 4;
}
+if (defined($opt_r)) { # residuals filters
+ $processing_options .= " -r $opt_r";
+ ($residuals_rms_max,$residuals_delta_max) = split(',',$opt_r);
+ croak("$0: cannot decode -r $opt_r\n") unless ($residuals_rms_max > 0);
+}
+
($LADCP_firstBin,$LADCP_lastBin) = split(',',$opt_b); # select valid bins
croak("$0: cannot decode -b $opt_b\n")
unless (numberp($LADCP_firstBin) &&
@@ -430,6 +467,7 @@
croak("$0: \$out_basename undefined\n") # plotting routines use this to label the plots
unless defined($out_basename);
+&antsAddParams('RDI_Coords::binMapping',$RDI_Coords::binMapping);
&antsAddParams('processing_options',$processing_options);
&antsAddParams('out_basename',$out_basename);
&antsAddParams('profile_id',$PROF,'run_label',$RUN);
@@ -467,7 +505,7 @@
printf(LOGF @msg);
print(LOGF "\n");
}
- return if ($opt_v == 0);
+ return if ($opt_v == 0) || ($lvl == 0);
if ($opt_v == 1) {
print(STDERR "$LADCP_file: ") if defined($LADCP_file);
@@ -502,7 +540,11 @@
progress("Reading LADCP data from <$LADCP_file>...\n");
readData($LADCP_file,\%LADCP);
-progress("\t%d ensembles\n",scalar(@{$LADCP{ENSEMBLE}}));
+if ($LADCP{BEAM_COORDINATES}) {
+ progress("\t%d ensembles (beam coordinates)\n",scalar(@{$LADCP{ENSEMBLE}}));
+} else {
+ progress("\t%d ensembles (Earth coordinates)\n",scalar(@{$LADCP{ENSEMBLE}}));
+}
error("$LADCP_file: cannot process multi-ping ensembles\n")
unless ($LADCP{PINGS_PER_ENSEMBLE} == 1);
@@ -548,6 +590,7 @@
#------------------------------------------------------------
# Edit beam-velocity data
# 0) beam-vel mask on -m
+# mask file has three columns: from_ens to_ens ignore_beam
# 1) correlation threshold
#------------------------------------------------------------
@@ -596,8 +639,7 @@
}
error("$LADCP_file: no valid data\n") unless ($nvv > 0);
progress("\tcorrelation threshold (-c %d counts): %d velocites removed (%d%% of total)\n",$opt_c,$cte,round(100*$cte/$nvv));
-# progress("\tattitude threshold (-t %d deg): %d velocites removed (%d%% of total)\n",$opt_t,$pte,round(100*$pte/$nvv));
-} else {
+} else { # if BEAM_COORDINATES
progress("Editing velocity data...\n");
error("$LADCP_file: cannot apply beamvel-mask $opt_m to earth-coordinate data\n")
if defined($opt_m);
@@ -608,7 +650,6 @@
}
error("$LADCP_file: no valid data\n") unless ($nvv > 0);
progress("\tcorrelation threshold (-c %d counts): %d velocites removed (%d%% of total)\n",$opt_c,$cte,round(100*$cte/$nvv));
-# progress("\tattitude threshold (-t %d deg): %d velocites removed (%d%% of total)\n",$opt_t,$pte,round(100*$pte/$nvv));
}
#-------------------------------------------------------------------
@@ -630,29 +671,40 @@
}
$nvw = 0;
for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
+ $LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH} =
+ gimbal_pitch($LADCP{ENSEMBLE}[$ens]->{PITCH},$LADCP{ENSEMBLE}[$ens]->{ROLL});
+
for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
if ($bad_beam) {
undef($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$bad_beam-1]);
undef($LADCP{ENSEMBLE}[$ens]->{BT_VELOCITY}[$bin][$bad_beam-1]);
}
- ($LADCP{ENSEMBLE}[$ens]->{U}[$bin],
- $LADCP{ENSEMBLE}[$ens]->{V}[$bin],
- $LADCP{ENSEMBLE}[$ens]->{W}[$bin],
- $LADCP{ENSEMBLE}[$ens]->{ERRVEL}[$bin]) =
- velBeamToEarth(\%LADCP,$ens,@{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]});
+ ($LADCP{ENSEMBLE}[$ens]->{INTERP_U}[$bin],
+ $LADCP{ENSEMBLE}[$ens]->{INTERP_V}[$bin],
+ $LADCP{ENSEMBLE}[$ens]->{INTERP_W}[$bin],
+ $LADCP{ENSEMBLE}[$ens]->{INTERP_ERRVEL}[$bin]) = earthVels(\%LADCP,$ens,$bin);
+ ($LADCP{ENSEMBLE}[$ens]->{V12}[$bin],$LADCP{ENSEMBLE}[$ens]->{W12}[$bin],
+ $LADCP{ENSEMBLE}[$ens]->{V34}[$bin],$LADCP{ENSEMBLE}[$ens]->{W34}[$bin]) = BPEarthVels(\%LADCP,$ens,$bin);
+ }
+
+ for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
+ $LADCP{ENSEMBLE}[$ens]->{U}[$bin] = $LADCP{ENSEMBLE}[$ens]->{INTERP_U}[$bin];
+ $LADCP{ENSEMBLE}[$ens]->{V}[$bin] = $LADCP{ENSEMBLE}[$ens]->{INTERP_V}[$bin];
+ $LADCP{ENSEMBLE}[$ens]->{W}[$bin] = $LADCP{ENSEMBLE}[$ens]->{INTERP_W}[$bin];
+ $LADCP{ENSEMBLE}[$ens]->{ERRVEL}[$bin] = $LADCP{ENSEMBLE}[$ens]->{INTERP_ERRVEL}[$bin];
+ undef($LADCP{ENSEMBLE}[$ens]->{INTERP_U}[$bin]);
+ undef($LADCP{ENSEMBLE}[$ens]->{INTERP_V}[$bin]);
+ undef($LADCP{ENSEMBLE}[$ens]->{INTERP_W}[$bin]);
+ undef($LADCP{ENSEMBLE}[$ens]->{INTERP_ERRVEL}[$bin]);
+
if (defined($LADCP{ENSEMBLE}[$ens]->{W}[$bin])) {
$per_bin_nsamp[$bin]++;
$nvw++;
}
- $LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH} =
- gimbal_pitch($LADCP{ENSEMBLE}[$ens]->{PITCH},$LADCP{ENSEMBLE}[$ens]->{ROLL});
+
$LADCP{ENSEMBLE}[$ens]->{U_UNEDITED}[$bin] = $LADCP{ENSEMBLE}[$ens]->{U}[$bin];
$LADCP{ENSEMBLE}[$ens]->{V_UNEDITED}[$bin] = $LADCP{ENSEMBLE}[$ens]->{V}[$bin];
$LADCP{ENSEMBLE}[$ens]->{W_UNEDITED}[$bin] = $LADCP{ENSEMBLE}[$ens]->{W}[$bin];
-
- ($LADCP{ENSEMBLE}[$ens]->{V12}[$bin],$LADCP{ENSEMBLE}[$ens]->{W12}[$bin],
- $LADCP{ENSEMBLE}[$ens]->{V34}[$bin],$LADCP{ENSEMBLE}[$ens]->{W34}[$bin]) =
- velBeamToBPEarth(\%LADCP,$ens,@{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]});
}
}
progress("\t$nvw valid velocities in bins $LADCP_firstBin-$LADCP_lastBin\n");
@@ -665,6 +717,8 @@
progress("Counting valid vertical velocities...\n");
$nvw = 0;
for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
+ $LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH} =
+ gimbal_pitch($LADCP{ENSEMBLE}[$ens]->{PITCH},$LADCP{ENSEMBLE}[$ens]->{ROLL});
for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
($LADCP{ENSEMBLE}[$ens]->{U}[$bin],
$LADCP{ENSEMBLE}[$ens]->{V}[$bin],
@@ -681,9 +735,9 @@
$LADCP{ENSEMBLE}[$ens]->{V12}[$bin] = $LADCP{ENSEMBLE}[$ens]->{V34}[$bin] = nan;
# for 3-beam solutions, the following code sets w12 and w34 equal to w
- my($delta) = $LADCP{ENSEMBLE}[$ens]->{ERRVEL}[$bin] * sqrt(2)/2 * tan(rad($LADCP{BEAM_ANGLE}));
- $LADCP{ENSEMBLE}[$ens]->{W12}[$bin] = $LADCP{ENSEMBLE}[$ens]->{W}[$bin] - $delta;
- $LADCP{ENSEMBLE}[$ens]->{W34}[$bin] = $LADCP{ENSEMBLE}[$ens]->{W}[$bin] + $delta;
+ my($delta) = $LADCP{ENSEMBLE}[$ens]->{ERRVEL}[$bin] / sqrt(2) / tan(rad($LADCP{BEAM_ANGLE}));
+ $LADCP{ENSEMBLE}[$ens]->{W12}[$bin] = $LADCP{ENSEMBLE}[$ens]->{W}[$bin] + $delta/2;
+ $LADCP{ENSEMBLE}[$ens]->{W34}[$bin] = $LADCP{ENSEMBLE}[$ens]->{W}[$bin] - $delta/2;
}
}
progress("\t$nvw valid velocities in bins $LADCP_firstBin-$LADCP_lastBin\n");
@@ -844,8 +898,11 @@
&antsAddParams('lat',$P{lat}) if defined($P{lat});
&antsAddParams('lon',$P{lon}) if defined($P{lon});
-($CTD_elapsed,$CTD_depth,$CTD_svel,$CTD_w) = &fnr('elapsed','depth','ss','w');
-$CTD_temp = &fnrNoErr('temp');
+($CTD_elapsed,$CTD_depth,$CTD_svel) = &fnr('elapsed','depth','ss');
+$CTD_w = &fnrNoErr('w_CTD');
+$CTD_w = &fnr('w') unless defined($CTD_w);
+$CTD_winchvel = &fnrNoErr('w_winch');
+$CTD_temp = &fnrNoErr('temp');
warning(0,"no CTD temperature --- using ADCP temperature instead => Sv degraded!\n",$s)
unless defined($CTD_temp);
@@ -857,7 +914,8 @@
push(@{$CTD{DEPTH}}, $ants_[0][$CTD_depth]);
push(@{$CTD{SVEL}}, $ants_[0][$CTD_svel]);
push(@{$CTD{W}}, $ants_[0][$CTD_w]);
- push(@{$CTD{TEMP}}, $ants_[0][$CTD_temp]) if defined($CTD_temp);
+ push(@{$CTD{TEMP}}, $ants_[0][$CTD_temp]) if defined($CTD_temp);
+ push(@{$CTD{W_WINCH}},$ants_[0][$CTD_winchvel]) if defined($CTD_winchvel);
if ($ants_[0][$CTD_depth] > $CTD_maxdepth) {
$CTD_maxdepth = $ants_[0][$CTD_depth];
$CTD_atbottom = $#{$CTD{DEPTH}};
@@ -1091,10 +1149,9 @@
&antsAddParams('rms_w_reflr_err',sqrt($sumWsq/$nWsq),'rms_w_reflr_err_interior',sqrt($sumWsqI/$nWsqI));
progress("\t%.2f cm/s rms reference-layer w_ocean, %.2f cm/s away from boundaries\n",
100*sqrt($sumWsq/$nWsq),100*sqrt($sumWsqI/$nWsqI));
-# if (sqrt($sumWsqI/$nWsqI) > 0.05) {
-# warning(0,"%.2f cm/s (large) reference-layer w_ocean away from boundaries\n",100*sqrt($sumWsqI/$nWsqI));
-# } els
- if (sqrt($sumWsqI/$nWsqI) > 0.15) {
+ if (sqrt($sumWsqI/$nWsqI) > 0.05) {
+ warning(0,"%.2f cm/s (large) reference-layer w_ocean away from boundaries\n",100*sqrt($sumWsqI/$nWsqI));
+ } elsif (sqrt($sumWsqI/$nWsqI) > 0.15) {
warning(2,"%.2f cm/s (implausibly large) reference-layer w_ocean away from boundaries\n",100*sqrt($sumWsqI/$nWsqI));
}
} elsif ($nWsq > 0) {
@@ -1111,8 +1168,9 @@
progress("Calculating volume-scattering coefficients (Sv)...\n");
calc_backscatter_profs($firstGoodEns,$lastGoodEns);
-if (@nSv) {
+if (defined($Sv_ref_bin) && @nSv) {
progress("\tapplying empirical Sv correction\n");
+ &antsAddParams('Sv_correction','empirical');
correct_backscatter($firstGoodEns,$lastGoodEns);
}
@@ -1207,6 +1265,7 @@
&antsAddParams('sidelobe_editing','seabed');
}
+ $PPI_editing |= eval($PPI_seabed_editing_required);
if ($PPI_editing) {
&antsAddParams('PPI_editing','seabed');
&antsAddParams('PPI_extend_upper_limit',$PPI_extend_upper_limit)
@@ -1267,7 +1326,7 @@
#----------------------------------------------------------------------
# Data Editing after LADCP and CTD data have been merged
# 1) surface layer editing
-# 2) user-supplied $edit_data_hook
+# 2) Execute user-supplied $edit_data_hook
#----------------------------------------------------------------------
progress("Removing data from instrument at surface...\n");
@@ -1321,7 +1380,7 @@
progress("Creating binned profiles at ${opt_o}m resolution...\n");
-&antsAddParams('outgrid_dz',$opt_o,'outgrid_minsamp',$opt_k); # used by LADCP_w_regrid
+&antsAddParams('outgrid_dz',$opt_o,'outgrid_minsamp',$opt_k); # used by LADCP_w_postproc
&antsAddParams('outgrid_firstbin',$outGrid_firstBin,'outgrid_lastbin',$outGrid_lastBin);
my($min_depth) = 9e99;
@@ -1471,7 +1530,7 @@
#------------------------------------------------------------------------------------------------------
unless ($opt_q) {
- progress("Removing ping-coherent residuals...\n");
+ progress("Subtracting ping-coherent residuals...\n");
for ($ens=$firstGoodEns; $ens<=$lastGoodEns; $ens++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
@@ -1539,30 +1598,17 @@
#----------------------------------------------------------------------
if (defined($opt_r)) {
- progress("Removing ensembles with large residuals...\n");
+ progress("Applying residuals filters...\n");
- $nerm = 0;
- for ($ens=$firstGoodEns; $ens<=$lastGoodEns; $ens++) {
- next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
+ progress("\trms residuals > $residuals_rms_max: ");
+ my($nerm) = editResiduals_rmsMax($firstGoodEns,$lastGoodEns,$residuals_rms_max);
+ progress("$nerm ensembles removed (%d%% of total)\n",round(100*$nerm/($lastGoodEns-$firstGoodEns+1)));
- my($sum) = my($n) = 0; # calculate rms residual
- 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;
- 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++;
- }
+ if (defined($residuals_delta_max)) { # filter based on difference between
+ progress("\tbeampair residuals difference > $residuals_delta_max [m/s]: ");
+ my($nvrm) = editResiduals_deltaMax($firstGoodEns,$lastGoodEns,$residuals_delta_max);
+ progress("$nvrm velocities removed (%d%% of total in bins $LADCP_firstBin-$LADCP_lastBin)\n",round(100*$nvrm/$nvw));
}
- progress("\t$nerm ensembles removed (%d%% of total)\n",round(100*$nerm/($lastGoodEns-$firstGoodEns+1)));
progress("\tre-binning profile data...\n");
undef(%DNCAST); undef(%UPCAST);
@@ -1807,7 +1853,7 @@
progress("Writing vertical-velocity data to ");
@antsNewLayout = ('ensemble','bin','elapsed','depth','CTD_depth','downcast',
'w','w12','w34','v12','v34','residual','residual12','residual34',
- 'CTD_w','CTD_w_t','CTD_w_tt','LADCP_w',
+ 'CTD_w','CTD_w_t','CTD_w_tt','LADCP_w','LADCP_reflr_w','winch_w',
'errvel','correlation','echo_amplitude','Sv',
'pitch','roll','tilt','heading','3_beam','svel');
@@ -1846,6 +1892,8 @@
$CTD{W_t}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$CTD{W_tt}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_W}[$bin],
+ $LADCP{ENSEMBLE}[$ens]->{REFLR_W},
+ $CTD{W_WINCH}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$LADCP{ENSEMBLE}[$ens]->{ERRVEL}[$bin],
median(@{$LADCP{ENSEMBLE}[$ens]->{CORRELATION}[$bin]}),
median(@{$LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin]}),
@@ -1884,6 +1932,8 @@
$CTD{W_t}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$CTD{W_tt}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_W}[$bin],
+ $LADCP{ENSEMBLE}[$ens]->{REFLR_W},
+ $CTD{W_WINCH}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$LADCP{ENSEMBLE}[$ens]->{ERRVEL}[$bin],
median(@{$LADCP{ENSEMBLE}[$ens]->{CORRELATION}[$bin]}),
median(@{$LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin]}),
--- a/LADCP_w_postproc Wed May 25 12:14:29 2016 -0400
+++ b/LADCP_w_postproc Fri Aug 05 11:02:51 2016 -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: Tue May 24 22:48:51 2016
+# dlm: Thu May 26 18:21:11 2016
# (c) 2015 A.M. Thurnherr
-# uE-Info: 602 50 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 101 35 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
$antsSummary = 'edit and re-grid LADCP vertical-velocity samples';
@@ -64,6 +64,7 @@
# Mar 31, 2016: - changed version %PARAM
# Apr 14, 2016: - added profile id to warning messages
# May 24, 2016: - improved plot
+# May 26, 2016: - added automatic directory creation on -d
($ANTS) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
($WCALC) = ($0 =~ m{^(.*)/[^/]*$});
@@ -88,20 +89,25 @@
'[-p)lot <[%03d_wprof.eps]> [-b)t <wprof>]]',
'<DL.wsamp file> [UL.wsamp file] (or only <UL.wsamp file>)');
-$dual_head = (@ARGV==1); # single or dual head
+$dual_head = (@ARGV==1); # single or dual head
-$id = defined($opt_i) ? $opt_i : &antsParam('profile_id'); # ensure profile id
+$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);
-if (defined($opt_d)) { # different output directory
- croak("$opt_d: no such directory\n")
- unless (-d "$opt_d");
+if (defined($opt_d)) { # select output directory
+ unless (-d $opt_d) {
+ unless ($opt_d =~ m{/}) {
+ print(STDERR "Warning: Creating output sub-directory ./$opt_d\n");
+ mkdir($opt_d);
+ }
+ croak("$opt_d: no such directory\n") unless (-d $opt_d);
+ }
}
-&antsCardOpt(\$opt_s,300); # surface layer depth limit
+&antsCardOpt(\$opt_s,300); # surface layer depth limit
-if (defined($opt_v)) { # bin ranges with valid data to use
+if (defined($opt_v)) { # bin ranges with valid data to use
($fvBin,$lvBin,$UL_fvBin,$UL_lvBin) = split(/,/,$opt_v);
croak("$0: cannot decode -v $opt_v\n")
unless (defined($lvBin) && (!$dual_head || defined($UL_lvBin)));
@@ -114,8 +120,8 @@
'DL_last_valid_bin',&antsRequireParam('outgrid_lastbin'));
}
-if (defined($opt_w)) { # vertical-velocity fields to use
- ($Ddwf,$Duwf,$Udwf,$Uuwf) = split(/,/,$opt_w); # DL dc, DL uc, ...
+if (defined($opt_w)) { # vertical-velocity fields to use
+ ($Ddwf,$Duwf,$Udwf,$Uuwf) = split(/,/,$opt_w); # DL dc, DL uc, ...
croak("$0: cannot decode -w $opt_w\n")
unless (defined($Duwf) && (!$dual_head || defined($Uuwf)));
&antsAddParams('DL_dc_w_field',$Ddwf,'DL_uc_w_field',$Duwf,
@@ -124,14 +130,14 @@
($Ddwf,$Duwf,$Udwf,$Uuwf) = ('w','w','w','w');
}
-if (defined($opt_o)) { # output grid resolution
+if (defined($opt_o)) { # output grid resolution
$opt_o_override = 1;
&antsCardOpt($opt_o);
} else {
$opt_o = &antsRequireParam('outgrid_dz');
}
-if (defined($opt_k)) { # minimum number of required samples
+if (defined($opt_k)) { # minimum number of required samples
$opt_k_override = 1;
&antsCardOpt($opt_k);
} else {
@@ -588,7 +594,7 @@
GMT_pstext('-F+f12,Helvetica,coral+jTL -Gred');
} elsif ($dc_R < 0.5) { GMT_pstext('-F+f12,Helvetica,coral+jTL -Gyellow'); }
else { GMT_pstext('-F+f12,Helvetica,coral+jTL -Gwhite'); }
- printf(GMT "%f %f r = %.1f\n",0.64,0.01,$dc_R);
+ printf(GMT "%f %f r = %.1f\n",0.61,0.01,$dc_R);
}
GMT_pstext('-F+f12,Helvetica,coral+jTR -Gwhite');
printf(GMT "%f %f [%.1f/%.1f cm/s @~s@~\@-e/r\@-]\n",
@@ -599,7 +605,7 @@
GMT_pstext('-F+f12,Helvetica,SeaGreen+jTL -Gred');
} elsif ($uc_R < 0.5) { GMT_pstext('-F+f12,Helvetica,SeaGreen+jTL -Gyellow'); }
else { GMT_pstext('-F+f12,Helvetica,SeaGreen+jTL -Gwhite'); }
- printf(GMT "%f %f r = %.1f\n",0.64,0.05,$uc_R);
+ printf(GMT "%f %f r = %.1f\n",0.61,0.05,$uc_R);
}
GMT_pstext('-F+f12,Helvetica,SeaGreen+jTR -Gwhite');
printf(GMT "%f %f [%.1f/%.1f cm/s @~s@~\@-e/r\@-]\n",
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/default_output.pl Fri Aug 05 11:02:51 2016 -0400
@@ -0,0 +1,111 @@
+#======================================================================
+# D E F A U L T _ O U T P U T . P L
+# doc: Wed Jun 1 19:21:19 2016
+# dlm: Wed Jun 1 19:26:01 2016
+# (c) 2016 A.M. Thurnherr
+# uE-Info: 106 1 NIL 0 0 72 0 2 4 NIL ofnI
+#======================================================================
+
+# NOTES:
+# - as this file is executed after ProcessingParams, only undefined
+# variables are set
+
+# HISTORY:
+# Jun 1, 2016: - created from [default_paths.pl]
+
+#----------------------------------------------------------------------
+# Processing log (diagnostic messages) output
+#----------------------------------------------------------------------
+
+$out_log = "$log_dir/$out_basename.log"
+ unless defined($out_log);
+
+
+#----------------------------------------------------------------------
+# Vertical-velocity profile output and plots:
+#
+# Data:
+# *.wprof vertical velocity profiles
+#
+# Plots:
+# *_wprof.ps vertical velocity profiles (main output plot)
+#----------------------------------------------------------------------
+
+unless (@out_profile) {
+ push(@out_profile,"$data_dir/$out_basename.wprof");
+ push(@out_profile,"plot_wprof($plot_dir/${out_basename}_wprof.ps)")
+ if ($plotting_level > 0);
+}
+
+#--------------------------------------------------------------------------------------------------
+# 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 1: vertical velocity time-depth plot
+# *_residuals.ps 1: residual vertical velocity time-depth plot
+# *_residual_profs.ps 1: residuals binned in depth
+# *_backscatter.ps 1: volume scattering coefficient time-depth plot
+# *_attitude_res.ps 2: residuals binned wrt. pitch/roll
+# *_residuals12.ps 2: beampair <1,2> residual vertical velocity time-depth plot
+# *_residuals34.ps 2: beampair <3,4> residual vertical velocity time-depth plot
+# *_attitude_res.ps 2: residuals binned wrt. package attitude
+# *_acceleration_res.ps 2: residuals binned wrt. package acceleration derivative
+# *_correlation.ps 3: correlation time-depth plot
+#--------------------------------------------------------------------------------------------------
+
+unless (@out_wsamp) {
+ push(@out_wsamp,"$data_dir/$out_basename.wsamp");
+ #push(@out_wsamp,sprintf('dump_residual_profiles(%s/residuals/%03d)',$data_dir,$PROF));
+
+ if ($plotting_level > 0) {
+ push(@out_wsamp,"plot_wsamp($plot_dir/${out_basename}_wsamp.ps)");
+ push(@out_wsamp,"plot_residuals($plot_dir/${out_basename}_residuals.ps)");
+ push(@out_wsamp,"plot_residual_profs($plot_dir/${out_basename}_residual_profs.ps)");
+ push(@out_wsamp,"plot_backscatter($plot_dir/${out_basename}_backscatter.ps)");
+ if ($plotting_level > 1) {
+ push(@out_wsamp,"plot_residuals12($plot_dir/${out_basename}_residuals12.ps)");
+ push(@out_wsamp,"plot_residuals34($plot_dir/${out_basename}_residuals34.ps)");
+ push(@out_wsamp,"plot_attitude_residuals($plot_dir/${out_basename}_attitude_res.ps)");
+ push(@out_wsamp,"plot_acceleration_residuals($plot_dir/${out_basename}_acceleration_res.ps)");
+ if ($plotting_level > 2) {
+ 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
+#----------------------------------------------------------------------
+
+push(@out_timeseries,"$data_dir/$out_basename.tis")
+ unless (@out_timeseries);
+
+#----------------------------------------------------------------------
+# Per-bin vertical-velocity residuals (plot only)
+#----------------------------------------------------------------------
+
+push(@out_BR,"plot_mean_residuals($plot_dir/${out_basename}_bin_residuals.ps)")
+ unless (@out_BR);
+
+
+#----------------------------------------------------------------------
+# Time-lagging correlation statistics (plot only)
+#----------------------------------------------------------------------
+
+unless (@out_TL) {
+ push(@out_TL,"plot_time_lags($plot_dir/${out_basename}_time_lags.ps)")
+ if ($plotting_level > 0);
+}
+
+#----------------------------------------------------------------------
+
+1; # return true
+
--- a/default_paths.pl Wed May 25 12:14:29 2016 -0400
+++ b/default_paths.pl Fri Aug 05 11:02:51 2016 -0400
@@ -1,9 +1,9 @@
#======================================================================
# D E F A U L T _ P A T H S . P L
# doc: Tue Mar 29 07:09:52 2016
-# dlm: Wed May 18 20:23:32 2016
+# dlm: Wed Jun 1 19:20:45 2016
# (c) 2016 A.M. Thurnherr
-# uE-Info: 13 44 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 16 57 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -11,6 +11,9 @@
# May 18, 2016: - added new attitude and acceleration residuals plots
# - renamed mean_residual to bin_residuals plots
# - added new residual_profs plot
+# Jun 1, 2016: - added residuals12 plots
+# - added support for $plotting_level
+# - exported stuff to [default_output.pl]
#======================================================================
# ProcessingParams file selection
@@ -46,13 +49,6 @@
}
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");
@@ -60,78 +56,17 @@
}
error("$log_dir: no such directory\n") unless (-d $log_dir);
}
-
-
-#----------------------------------------------------------------------
-# Processing log (diagnostic messages) output
-#----------------------------------------------------------------------
-
-$out_log = "$log_dir/$out_basename.log";
-
+if ($plotting_level > 0) {
+ 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);
+ }
+}
#----------------------------------------------------------------------
-# 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
-# *_attitude_res.ps residuals binned wrt. pitch/roll
-# *_res_profs.ps residuals binned in depth
-# *_acceleration_res.ps OPTIONAL: residuals binned wrt. package acceleration derivative
-# *_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_attitude_residuals($plot_dir/${out_basename}_attitude_res.ps)");
-push(@out_wsamp,"plot_residual_profs($plot_dir/${out_basename}_residual_profs.ps)");
-#push(@out_wsamp,"plot_acceleration_residuals($plot_dir/${out_basename}_acceleration_res.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}_bin_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 Wed May 25 12:14:29 2016 -0400
+++ b/defaults.pl Fri Aug 05 11:02:51 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: Tue Mar 29 07:23:24 2016
+# dlm: Mon Jun 6 22:07:19 2016
# (c) 2011 A.M. Thurnherr
-# uE-Info: 74 39 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 179 18 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -72,6 +72,10 @@
# Mar 19, 2016: - improved docu
# Mar 29, 2016: - moved out dir creation to [LADCP_w_ocean]
# - added opt_r support
+# May 26, 2016: - added RDI_Coords::binMapping
+# May 28, 2016: - added delta-residual filter (-r)
+# Jun 1, 2016: - added $plotting_level
+# Jun 2, 2016: - added $tilt_correction_*
#======================================================================
# Output Log Files
@@ -87,17 +91,35 @@
&antsCardOpt(\$opt_v,$ENV{VERB});
$opt_v = 1 unless numberp($opt_v);
+#======================================================================
+# Output Plots
+# - there are 3 plotting levels
+# 0 : suppress all plots
+# 1* : produce default plots; *DEFAULT
+# 2 : produce default and diagnostic plots
+# >2: produce all plots, including useless ones
+#======================================================================
+
+$plotting_level = 1;
+
#======================================================================
-# Data Input
+# Input Data
#======================================================================
# Set $opt_4 to 1 (or use the -4 option) to suppress 3-beam LADCP
-# solutions
+# solutions. This only has an effect for beam-coordinate data.
#$opt_4 = 1;
+# Set $RDI_Coords::binMapping to 'none' or use -d to disable
+# linear bin interpolation (which is better than bin mapping).
+# This only has an effect for beam-coordinate data.
+
+$RDI_Coords::binMapping = $opt_d ? 'none' : 'linterp';
+
+
# The following variables allow bias-correcting the attiude
# sensors.
# NB: heading is not used for vertical-velocity processing!
@@ -125,10 +147,15 @@
# Data Editing
#======================================================================
-# The following sets the max allowable rms residual w per ensemble;
-# data from ensembles with larger rms residuals are discarded.
+# The following sets the max allowable rms residual w per ensemble, as
+# well as the max allowable difference between the two beam-pair
+# residuals. Measurements that fail either of these tests are are
+# discarded. The limiting values were chosed by inspection of
+# log files and diagnostic plots of a few example profiles. In case of
+# the delta-residual limit, histograms of this parameter show a very
+# steep cutoff at 0.05 cm/s for both IWISE and 2016_I08S data.
-&antsFloatOpt(\$opt_r,0.04);
+&antsFloatOpt(\$opt_r,'0.06,0.06');
# By default, ensembles with uncertain time-lagging are discarded.
@@ -147,25 +174,18 @@
&antsFloatOpt(\$opt_c,70);
-# The following sets the default limit for instrument attitude
-# (pitch/roll).
-#
-# The default value was established with IWISE profiles 004, 005, 045
-# and 049, which all show considerabe tilt-related discrepancies between
-# the corresponding 2-beam solutions. The first and second pair of
-# profiles were collected with 8 and 6m bins, respectively, without
-# bin re-mapping. The original default of 15 degrees led to large
-# beam-pair differences. Based on diagnostic plots it appears that
-# only tilt angles smaller than 9 degrees or so are satisfactory.
-# In case of the IWISE data set, such a tight constraint causes
-# too many data gaps. The compromise of 12 degrees seems to work
-# quite well, based on the p0 vs epsilon correlation across 5
-# data sets.
+# Instrument Tilt
#
-# NB: if this default is changed, the usage message in [LADCP_w]
-# needs to be updated as well.
+# It is not fully clear what tilt angles are acceptable for
+# obtaining good vertical velocities. Up to 2016 (V1.3) the
+# default limit was 12 degrees based on an analysis of
+# the 2010 IWISE data with inaccurate 2-beam transformations.
+# As re-processing with a limit of 20 degrees improves the
+# agreement between DL and UL data (R = 0.77/0.67 => 0.79/0.74)
+# the limit was changed to 22 degrees, the same used
+# in Martin Visbeck's inversion code.
-&antsFloatOpt(\$opt_t,12);
+&antsFloatOpt(\$opt_t,22);
# The following sets the default error velocity limit; measurements
@@ -224,8 +244,10 @@
# Previous Ping Interference editing as described in [edit_data.pl]
# - enabled by default for WH150 data
-# - the variable defines a string with a perl expression, which is
-# evaluated once the data are loaded
+# - PPI_seabed_editing_required defines a string with a perl expression
+# that is evaluated once the data are loaded; if true, seabed PPI
+# editing is enabled
+# - to enable PPI editing without condition, set $PPI_editing = 1;
# - 2014 CLIVAR P16 #47 has a slight discontinuity at 4000m; this
# discontinuity is there without PPI filtering but gets slightly
# worse with PPI filtering. Setting $PPI_extend_upper_limit to
@@ -236,8 +258,9 @@
# set by the shortest acoustic path between the ADCP and the
# seabed.
-$PPI_editing_required = '($LADCP{BEAM_FREQUENCY} < 300)';
+$PPI_seabed_editing_required = '($LADCP{BEAM_FREQUENCY} < 300)';
+#$PPI_editing = 1; # uncomment to enable PPI always
#$PPI_extend_upper_limit = 1.03; # see comments above
@@ -311,7 +334,9 @@
# bin is chosen to construct a reference profile for Sv. The bin number
# is automatically increased if the selected bin does not contain valid
# data, i.e. the default value of 1 ensures that the closest valid bin
-# is used to construct the reference profile.
+# is used to construct the reference profile. The empirical correction
+# causes artifacts every 100m. To disable the empirical
+# correction, undefine the following variable.
$Sv_ref_bin = 1;
--- a/edit_data.pl Wed May 25 12:14:29 2016 -0400
+++ b/edit_data.pl Fri Aug 05 11:02:51 2016 -0400
@@ -1,9 +1,9 @@
#======================================================================
# E D I T _ D A T A . P L
# doc: Sat May 22 21:35:55 2010
-# dlm: Tue May 24 16:36:37 2016
+# dlm: Mon Jun 6 21:13:28 2016
# (c) 2010 A.M. Thurnherr
-# uE-Info: 374 18 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 543 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -35,18 +35,25 @@
# Jan 23, 2016: - added &editBadTimeLagging()
# May 18, 2016: - removed assumption of 1500m/s soundspeed setting
# May 24, 2016: - calc_binDepths() -> binDepths()
+# May 28, 2016: - added editResiduals_rmsMax, editResiduals_deltaMax
+# Jun 2, 2016: - added applyTiltCorrection()
+# - maded editCorr_Earthcoords() less conservative
+# - verified that removed velocities are counted correctly
+# Jun 3, 2016: - BUG: applyTiltCorrection() did not use gimbal_pitch
+# Jun 6, 2016: - removed applyTiltCorrection()
# NOTES:
-# - editCorr_Earthcoords() is overly conservative and removed most
-# or all 3-beam solutions
# - all bins must be edited (not just the ones between $LADCP_firstBin
# and $LADCP_lastBin to allow reflr calculations to use bins outside
# this range (ONLY FOR BEAM-COORD EDITS)
# - however, to make the stats work, only the edited velocities
-# inside the bin range are counted
+# inside the bin range are counted for those edit functions that
+# report their stats wrt $nvw (for those which use $nvv,
+# all velocities must be counted)
#======================================================================
# correctAttitude($ens,$pitch_bias,$roll_bias,$heading_bias)
+# - this is called before gimbal_pitch is calculated
#======================================================================
sub correctAttitude($$$$)
@@ -59,6 +66,10 @@
#======================================================================
# $vv = countValidVels($ens)
+#
+# NOTES:
+# - in case of Earth coords, this counts the velocity components
+# (including errvel)
#======================================================================
sub countValidBeamVels($)
@@ -79,7 +90,8 @@
# $removed = editCorr($ens,$threshold)
#
# NOTES:
-# - called before Earth vels have been calculated
+# - called before Earth vels are calculated
+# - count removed velocities in all bins
#======================================================================
sub editCorr($$)
@@ -102,12 +114,11 @@
# $removed = editCorr_Earthcoords($ens,$threshold)
#
# NOTES:
-# - if any of the 4 beam correlations is below the threshold,
+# - if any of the used correlations is below the threshold,
# the entire velocity is removed
-# - this implies that (most? all?) three-beam solutions will
-# be edited out, which is overly conserative
-# - a potentially better algorithm (used in LADCPproc) ignores the
-# lowest correlation in all 3-beam solutions
+# - for three-beam solutions two correlations must fail the
+# test
+# - count velocities in all bins
#======================================================================
sub editCorr_Earthcoords($$)
@@ -116,14 +127,14 @@
my($nrm) = 0;
for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
- my($beam);
- for ($beam=0; $beam<4; $beam++) {
- last unless ($LADCP{ENSEMBLE}[$ens]->{CORRELATION}[$bin][$beam] >= $lim);
+ my($nBad) = 0;
+ for (my($beam)=0; $beam<4; $beam++) {
+ $nBad++ unless ($LADCP{ENSEMBLE}[$ens]->{CORRELATION}[$bin][$beam] > $lim);
}
- if ($beam < 4) {
- for (my($c)=0; $c<4; $c++) {
- next unless defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$c]);
- undef($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$c]);
+ if ($nBad-$LADCP{ENSEMBLE}[$ens]->{THREE_BEAM}[$bin] > 0) {
+ for (my($beam)=0; $beam<4; $beam++) {
+ next unless defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$beam]);
+ undef($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$beam]);
$nrm++;
}
}
@@ -137,6 +148,7 @@
# NOTES:
# - called after Earth vels have been calculated
# - sets TILT field for each ensemble as a side-effect
+# - count all removed velocities
#======================================================================
sub editTilt($$)
@@ -209,6 +221,7 @@
#
# NOTES:
# - call after Earth vels have been calculated
+# - count only removed vels in selected bin range
#======================================================================
sub editTruncRange($$)
@@ -232,6 +245,7 @@
# - remove data from far bins
# - only bins in valid range are considered here, because
# $per_bin_nsamp is only defined for those
+# - only velocities from bins in valid range are counted
#======================================================================
sub editFarBins($$)
@@ -255,6 +269,8 @@
# 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
+# - all velocities are counted, even those outside valid bin range,
+# because the %age is not reported
#======================================================================
sub editSideLobes($$$)
@@ -290,6 +306,7 @@
# ($nvrm,$nerm) = editPPI($fromEns,$toEns,$water_depth)
#
# NOTES:
+# - only velocities in good bin range are removed/counted
# - for UL, water_depth == undef; for DL water_depth is always defined,
# or else editPPI is not called
# - when this code is executed a suitable UL or DL depth-average-soundspeed
@@ -422,6 +439,9 @@
#===============================================================================
# $nerm = editBadTimeLagging($fromEns,$toEns,$good_from_elapsed1,$good_to_elapsed1,...)
+#
+# NOTES:
+# - deleted velocities are not counted
#===============================================================================
sub editBadTimeLagging($$@)
@@ -458,5 +478,69 @@
return $nerm;
}
+#======================================================================
+# $nerm = editResiduals_rmsMax($fe,$te,$max_val)
+#
+# NOTES:
+# - removed velocities are not counted
+#======================================================================
+
+sub editResiduals_rmsMax($$$)
+{
+ my($fe,$te,$limit) = @_;
+ my($nerm) = 0;
+ for (my($ens)=$fe; $ens<=$te; $ens++) {
+ next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
+ my($sum) = my($n) = 0; # calculate rms residual
+ 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;
+ 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) > $limit) { # ensemble is bad
+ undef($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
+ $nerm++;
+ }
+ }
+ return $nerm;
+}
+
+#======================================================================
+# $nerm = editResiduals_deltaMax($fe,$te,$max_val)
+# - delta residual = delta beampair w => equal to scaled error velocity?
+# - sharp cutoff near 5cm/s for std parameters (0.1 m/s error velocity
+# filter) in several data sets
+# - samples with large residuals differences are clear outliers in
+# the residuals vs tilt plots => obvious to remove
+# - how are large delta res possible given the errvel limit???
+# - only valid bin range is edited/counted
+#======================================================================
+
+sub editResiduals_deltaMax($$$)
+{
+ my($fe,$te,$limit) = @_;
+ my($nvrm) = 0;
+ for (my($ens)=$fe; $ens<=$te; $ens++) {
+ next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
+ 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($Dr) = abs($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin] -
+ $LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin]);
+ if ($Dr > $limit) {
+ undef($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
+ $nvrm++;
+ }
+ }
+ }
+ return $nvrm;
+}
+
+#======================================================================
1;
--- a/plot_attitude_biases_w.pl Wed May 25 12:14:29 2016 -0400
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,238 +0,0 @@
-#======================================================================
-# P L O T _ A T T I T U D E _ B I A S E S _ W . P L
-# doc: Sun May 15 16:08:59 2016
-# dlm: Tue May 24 16:38:31 2016
-# (c) 2016 A.M. Thurnherr
-# uE-Info: 15 31 NIL 0 0 72 2 2 4 NIL ofnI
-#======================================================================
-
-# HISTORY:
-# May 15, 2016: - created from [plot_mean_residuals.pl]
-# May 16, 2016: - continued
-# May 17, 2016: - renamed from [plot_attitude_biases.pl]
-# May 18, 2016: - added version
-# - expunged $realLastGoodEns
-# May 19, 2016: - added notes about wrong beam plane
-# May 24, 2016: - calc_binDepths() -> binDepths()
-
-# IMPORTANT NOTE:
-# - the variables prefixed with p/r refer to beam-pairs 1,2 and 3,4 respectively,
-# i.e. the p variables correspond to the roll plane and the r variables
-# correspond to the pitch plane
-
-require "$ANTS/libGMT.pl";
-
-sub plot_attitude_biases_w($)
-{
- my($pfn) = @_; # plot file name
-
- my($xmin) = -round($opt_t); # full pitch range
- my($xmax) = round($opt_t);
- my($ymin) = -0.05;
- my($ymax) = 0.05;
-
- my($min_thin) = 200;
- my($min_fat) = 500;
- my($btstrp_ndraw) = 20;
- my($excluded_surf_layer) = 150;
-
- #--------------------------------------------------------
- # Bin-Average & Create Histogram from Beampair Velocities
- # - use 1-degree bins for simplicity
- # - use gimbal pitch (not measured pitch)
- #--------------------------------------------------------
-
- my(@pHistDC,@rHistDC,@pSumDC,@rSumDC,$pHistDC,$rHistDC,$pSumDC,$rSumDC);
- my(@pHistUC,@rHistUC,@pSumUC,@rSumUC,$pHistUC,$rHistUC,$pSumUC,$rSumUC);
- my(@pValsDC,@rValsDC,@pValsUC,@rValsUC,$mode);
- for (my($e)=$firstGoodEns; $e<=$lastGoodEns; $e++) {
- next unless numberp($LADCP{ENSEMBLE}[$e]->{CTD_DEPTH});
- my(@bindepth) = binDepths($e);
- for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
- next if ($bindepth[$bin] <= $excluded_surf_layer);
- next unless ($bin+1>=$outGrid_firstBin && $bin+1<=$outGrid_lastBin);
- next unless numberp($LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W12}[$bin]) &&
- numberp($LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W34}[$bin]);
- my($pi) = int($LADCP{ENSEMBLE}[$e]->{GIMBAL_PITCH}+$opt_t); # pitch/roll indices
- my($ri) = int($LADCP{ENSEMBLE}[$e]->{ROLL}+$opt_t);
- if ($e < $LADCP_atbottom) { # downcast
- push(@{$pValsDC[$pi]},$LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W12}[$bin]);
- push(@{$rValsDC[$ri]},$LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W34}[$bin]);
- $pSumDC += $LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W12}[$bin]; $pHistDC++;
- $rSumDC += $LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W34}[$bin]; $rHistDC++;
- $pSumDC[$pi] += $LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W12}[$bin]; $pHistDC[$pi]++;
- $rSumDC[$ri] += $LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W34}[$bin]; $rHistDC[$ri]++;
- $mode = $pHistDC[$pi] if ($pHistDC[$pi] > $mode);
- $mode = $rHistDC[$ri] if ($rHistDC[$ri] > $mode);
- } else { # upcast
- push(@{$pValsUC[$pi]},$LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W12}[$bin]);
- push(@{$rValsUC[$ri]},$LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W34}[$bin]);
- $pSumUC += $LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W12}[$bin]; $pHistUC++;
- $rSumUC += $LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W34}[$bin]; $rHistUC++;
- $pSumUC[$pi] += $LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W12}[$bin]; $pHistUC[$pi]++;
- $rSumUC[$ri] += $LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W34}[$bin]; $rHistUC[$ri]++;
- $mode = $pHistUC[$pi] if ($pHistUC[$pi] > $mode);
- $mode = $rHistUC[$ri] if ($rHistUC[$ri] > $mode);
- }
- }
- }
-
- #----------
- # Plot Data
- #----------
-
- my($R) = "-R$xmin/$xmax/$ymin/$ymax"; # begin plot
- GMT_begin($pfn,'-JX10/10',$R,"-P -Bg5f1a5:'Pitch/Roll [\260]':/g0.01f0.001a0.01:'Beam-Plane Vertical Velocity [m/s]':WeSn");
-
- # ZERO LINE
- GMT_psxy('-W4,CornflowerBlue');
- print(GMT "$xmin 0\n$xmax 0\n");
-
- # DC BEAMS 1,2
- GMT_psxy('-Ey0.2/2,coral');
- for (my($i)=0; $i<2*round($opt_t); $i++) { # error bars
- next unless ($pHistDC[$i] >= $min_fat);
- my($minLim,$maxLim) = &bootstrap($btstrp_ndraw,0.95,\&avg,@{$pValsDC[$i]}); # 95% bootstrap conf limits
- printf(GMT "%f %f %f\n",$i-round($opt_t)-0.3,($maxLim+$minLim)/2,($maxLim-$minLim)/2);
- }
- GMT_psxy('-Ey0.2/1,coral');
- for (my($i)=0; $i<2*round($opt_t); $i++) {
- next unless ($pHistDC[$i]>=$min_thin && $pHistDC[$i]<$min_fat);
- my($minLim,$maxLim) = &bootstrap($btstrp_ndraw,0.95,\&avg,@{$pValsDC[$i]}); # 95% bootstrap conf limits
- printf(GMT "%f %f %f\n",$i-round($opt_t)-0.3,($maxLim+$minLim)/2,($maxLim-$minLim)/2);
- }
- GMT_psxy('-Sc0.15 -Gcoral');
- for (my($i)=0; $i<2*round($opt_t); $i++) {
- next unless ($pHistDC[$i] >= $min_thin);
- printf(GMT "%f %f\n",$i-round($opt_t)-0.3,$pSumDC[$i]/$pHistDC[$i]); # errorbar center symbol
- printf(GMT "%f %f\n",$i-round($opt_t)-0.3,$ymin+0.03*$pHistDC[$i]/$mode); # histogram symbol
- }
- if ($pHistDC) {
- GMT_psxy('-W1,coral,8_2:0'); # average bias (horizontal line);
- printf(GMT ">\n%f %f\n%f %f\n",-$opt_t,$pSumDC/$pHistDC,$opt_t,$pSumDC/$pHistDC)
- }
- GMT_psxy('-W2,coral,8_2:0');
- for (my($i)=0; $i<2*round($opt_t); $i++) {
- next unless ($pHistDC[$i] >= $min_thin);
- printf(GMT ">\n%f %f\n%f %f\n%f %f\n>\n%f %f\n%f %f\n", # histogram bar
- $i-round($opt_t)-0.3-0.5,$ymin,
- $i-round($opt_t)-0.3-0.5,$ymin+0.03*$pHistDC[$i]/$mode,
- $i-round($opt_t)-0.3+0.5,$ymin+0.03*$pHistDC[$i]/$mode,
- $i-round($opt_t)-0.3+0.5,$ymin,
- $i-round($opt_t)-0.3+0.5,$ymin+0.03*$pHistDC[$i]/$mode);
- }
-
- # DC BEAMS 3,4
- GMT_psxy('-Ey0.2/2,coral');
- for (my($i)=0; $i<2*round($opt_t); $i++) {
- next unless ($rHistDC[$i] >= $min_fat);
- my($minLim,$maxLim) = &bootstrap($btstrp_ndraw,0.95,\&avg,@{$rValsDC[$i]});
- printf(GMT "%f %f %f\n",$i-round($opt_t)-0.1,($maxLim+$minLim)/2,($maxLim-$minLim)/2);
- }
- GMT_psxy('-Ey0.2/1,coral');
- for (my($i)=0; $i<2*round($opt_t); $i++) {
- next unless ($rHistDC[$i]>=$min_thin && $rHistDC[$i]<$min_fat);
- my($minLim,$maxLim) = &bootstrap($btstrp_ndraw,0.95,\&avg,@{$rValsDC[$i]});
- printf(GMT "%f %f %f\n",$i-round($opt_t)-0.1,($maxLim+$minLim)/2,($maxLim-$minLim)/2);
- }
- GMT_psxy('-Sx0.25 -W2,coral');
- for (my($i)=0; $i<2*round($opt_t); $i++) {
- next unless ($rHistDC[$i] >= $min_thin);
- printf(GMT "%f %f\n",$i-round($opt_t)-0.1,$rSumDC[$i]/$rHistDC[$i]);
- printf(GMT "%f %f\n",$i-round($opt_t)-0.1,$ymin+0.03*$rHistDC[$i]/$mode);
- }
- if ($rHistDC) {
- GMT_psxy('-W1,coral,2_2:0');
- printf(GMT ">\n%f %f\n%f %f\n",-$opt_t,$rSumDC/$rHistDC,$opt_t,$rSumDC/$rHistDC);
- }
- GMT_psxy('-W2,coral,2_2:0');
- for (my($i)=0; $i<2*round($opt_t); $i++) {
- next unless ($rHistDC[$i] >= $min_thin);
- printf(GMT ">\n%f %f\n%f %f\n%f %f\n>\n%f %f\n%f %f\n",
- $i-round($opt_t)-0.1-0.5,$ymin,
- $i-round($opt_t)-0.1-0.5,$ymin+0.03*$rHistDC[$i]/$mode,
- $i-round($opt_t)-0.1+0.5,$ymin+0.03*$rHistDC[$i]/$mode,
- $i-round($opt_t)-0.1+0.5,$ymin,
- $i-round($opt_t)-0.1+0.5,$ymin+0.03*$rHistDC[$i]/$mode);
- }
-
- # UC BEAMS 1,2
- GMT_psxy('-Ey0.2/2,SeaGreen');
- for (my($i)=0; $i<2*round($opt_t); $i++) {
- next unless ($pHistUC[$i] >= $min_fat);
- my($minLim,$maxLim) = &bootstrap($btstrp_ndraw,0.95,\&avg,@{$pValsUC[$i]});
- printf(GMT "%f %f %f\n",$i-round($opt_t)+0.1,($maxLim+$minLim)/2,($maxLim-$minLim)/2);
- }
- GMT_psxy('-Ey0.2/1,SeaGreen');
- for (my($i)=0; $i<2*round($opt_t); $i++) {
- next unless ($pHistUC[$i]>=$min_thin && $pHistUC[$i]<$min_fat);
- my($minLim,$maxLim) = &bootstrap($btstrp_ndraw,0.95,\&avg,@{$pValsUC[$i]});
- printf(GMT "%f %f %f\n",$i-round($opt_t)+0.1,($maxLim+$minLim)/2,($maxLim-$minLim)/2);
- }
- GMT_psxy('-Sc0.15 -GSeaGreen');
- for (my($i)=0; $i<2*round($opt_t); $i++) {
- next unless ($pHistUC[$i] >= $min_thin);
- printf(GMT "%f %f\n",$i-round($opt_t)+0.1,$pSumUC[$i]/$pHistUC[$i]);
- printf(GMT "%f %f\n",$i-round($opt_t)+0.1,$ymin+0.03*$pHistUC[$i]/$mode);
- }
- if ($pHistUC) {
- GMT_psxy('-W1,SeaGreen,8_2:0');
- printf(GMT ">\n%f %f\n%f %f\n",-$opt_t,$pSumUC/$pHistUC,$opt_t,$pSumUC/$pHistUC);
- }
- GMT_psxy('-W2,SeaGreen,8_2:0');
- for (my($i)=0; $i<2*round($opt_t); $i++) {
- next unless ($pHistUC[$i] >= $min_thin);
- printf(GMT ">\n%f %f\n%f %f\n%f %f\n>\n%f %f\n%f %f\n",
- $i-round($opt_t)+0.1-0.5,$ymin,
- $i-round($opt_t)+0.1-0.5,$ymin+0.03*$pHistUC[$i]/$mode,
- $i-round($opt_t)+0.1+0.5,$ymin+0.03*$pHistUC[$i]/$mode,
- $i-round($opt_t)+0.1+0.5,$ymin,
- $i-round($opt_t)+0.1+0.5,$ymin+0.03*$pHistUC[$i]/$mode);
- }
-
- # UC BEAMS 3,4
- GMT_psxy('-Ey0.2/2,SeaGreen');
- for (my($i)=0; $i<2*round($opt_t); $i++) {
- next unless ($rHistUC[$i] >= $min_fat);
- my($minLim,$maxLim) = &bootstrap($btstrp_ndraw,0.95,\&avg,@{$rValsUC[$i]});
- printf(GMT "%f %f %f\n",$i-round($opt_t)+0.3,($maxLim+$minLim)/2,($maxLim-$minLim)/2);
- }
- GMT_psxy('-Ey0.2/1,SeaGreen');
- for (my($i)=0; $i<2*round($opt_t); $i++) {
- next unless ($rHistUC[$i]>=$min_thin && $rHistUC[$i]<$min_fat);
- my($minLim,$maxLim) = &bootstrap($btstrp_ndraw,0.95,\&avg,@{$rValsUC[$i]});
- printf(GMT "%f %f %f\n",$i-round($opt_t)+0.3,($maxLim+$minLim)/2,($maxLim-$minLim)/2);
- }
- GMT_psxy('-Sx0.25 -W2,SeaGreen');
- for (my($i)=0; $i<2*round($opt_t); $i++) {
- next unless ($rHistUC[$i] >= $min_thin);
- printf(GMT "%f %f\n",$i-round($opt_t)+0.3,$rSumUC[$i]/$rHistUC[$i]);
- printf(GMT "%f %f\n",$i-round($opt_t)+0.3,$ymin+0.03*$rHistUC[$i]/$mode);
- }
- if ($rHistUC) {
- GMT_psxy('-W1,SeaGreen,2_2:0');
- printf(GMT ">\n%f %f\n%f %f\n",-$opt_t,$rSumUC/$rHistUC,$opt_t,$rSumUC/$rHistUC);
- }
- GMT_psxy('-W2,SeaGreen,2_2:0');
- for (my($i)=0; $i<2*round($opt_t); $i++) {
- next unless ($rHistUC[$i] >= $min_thin);
- printf(GMT ">\n%f %f\n%f %f\n%f %f\n>\n%f %f\n%f %f\n",
- $i-round($opt_t)+0.3-0.5,$ymin,
- $i-round($opt_t)+0.3-0.5,$ymin+0.03*$rHistUC[$i]/$mode,
- $i-round($opt_t)+0.3+0.5,$ymin+0.03*$rHistUC[$i]/$mode,
- $i-round($opt_t)+0.3+0.5,$ymin,
- $i-round($opt_t)+0.3+0.5,$ymin+0.03*$rHistUC[$i]/$mode);
- }
-
- GMT_unitcoords(); # LABELS
- GMT_pstext('-F+f9,Helvetica,orange+jTR -N -Gwhite');
- print(GMT "0.99 0.01 V$VERSION\n");
-
- GMT_pstext('-F+f14,Helvetica,blue+jBL -N'); # profile id
- print(GMT "0.0 1.03 $P{out_basename} $P{run_label}\n");
-
- GMT_setR($R); # FINISH PLOT
- GMT_end();
-}
-
-1; # return true on require
--- a/plot_attitude_residuals.pl Wed May 25 12:14:29 2016 -0400
+++ b/plot_attitude_residuals.pl Fri Aug 05 11:02:51 2016 -0400
@@ -1,9 +1,9 @@
#======================================================================
# P L O T _ A T T I T U D E _ R E S I D U A L S . P L
# doc: Sun May 15 16:08:59 2016
-# dlm: Tue May 24 16:39:01 2016
+# dlm: Wed Jun 8 22:12:19 2016
# (c) 2016 A.M. Thurnherr
-# uE-Info: 54 24 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 92 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -15,11 +15,12 @@
# - expunged $realLastGoodEns
# May 19, 2016: - added notes about beam planes
# May 24, 2016: - calc_binDepths() -> binDepths()
+# Jun 8, 2016: - adapted to correct tilt axis
# IMPORTANT NOTE:
# - the variables prefixed with p/r refer to beam-pairs 1,2 and 3,4 respectively,
# i.e. the p variables correspond to the roll plane and the r variables
-# correspond to the pitch plane
+# correspond to the pitch plane. This is confusing.
require "$ANTS/libGMT.pl";
@@ -59,8 +60,8 @@
next unless numberp($LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W12}[$bin]) &&
numberp($LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W34}[$bin]);
- my($pi) = int($LADCP{ENSEMBLE}[$e]->{GIMBAL_PITCH}+$opt_t); # pitch/roll indices
- my($ri) = int($LADCP{ENSEMBLE}[$e]->{ROLL}+$opt_t);
+ my($ri) = int($LADCP{ENSEMBLE}[$e]->{GIMBAL_PITCH}+$opt_t); # pitch/roll indices
+ my($pi) = int($LADCP{ENSEMBLE}[$e]->{ROLL}+$opt_t);
my($bi) = $bindepth[$bin]/$opt_o;
if ($e < $LADCP_atbottom) { # downcast
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/plot_residuals12.pl Fri Aug 05 11:02:51 2016 -0400
@@ -0,0 +1,87 @@
+#======================================================================
+# P L O T _ R E S I D U A L S 1 2 . P L
+# doc: Wed Jun 1 19:05:22 2016
+# dlm: Wed Jun 1 19:35:43 2016
+# (c) 2016 A.M. Thurnherr
+# uE-Info: 82 67 NIL 0 0 72 0 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+# Jun 1, 2016: - created from [plot_residuals.pl]
+
+require "$ANTS/libGMT.pl";
+
+sub plot_residuals12($)
+{
+ my($pfn) = @_;
+
+ return unless ($P{max_depth});
+
+ my($xmin) = $P{min_ens}-0.5;
+ my($xmax) = $P{max_ens}+0.5;
+ my($ymin) = round(antsParam('min_depth')-25,50);
+ my($ymax) = ($P{water_depth} > 0) ?
+ round($P{water_depth}+25,50) :
+ round($P{max_depth}+$P{ADCP_bin_length}+25,50);
+
+ my($ens_width) = 10 / ($P{max_ens} - $P{min_ens} + 1);
+ my($bin_length) = 10 * $P{ADCP_bin_length} /
+ ($P{max_depth}-$P{min_depth}+$P{ADCP_bin_length});
+
+ my($R) = "-R$xmin/$xmax/$ymin/$ymax";
+ GMT_begin($pfn,'-JX10/-10',$R,'-P');
+
+ my($C) = "-C$WCALC/residuals.cpt";
+ GMT_psxy("$C -Sr");
+ for ($ens=$firstGoodEns; $ens<$LADCP_atbottom; $ens++) { # downcast
+ next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
+ my(@bindepth) = binDepths($ens);
+ for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
+ next unless ($bin+1>=$outGrid_firstBin && $bin+1<=$outGrid_lastBin);
+ next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
+ next unless numberp($LADCP{ENSEMBLE}[$ens]->{W12}[$bin]);
+ my($bi) = $bindepth[$bin]/$opt_o;
+ printf(GMT "%d %f %f $ens_width $bin_length\n",
+ $LADCP{ENSEMBLE}[$ens]->{NUMBER},
+ $bindepth[$bin],
+ $LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin] - $DNCAST{MEDIAN_W}[$bi]);
+ }
+ }
+ for ($ens=$LADCP_atbottom; $ens<=$lastGoodEns; $ens++) { # upcast
+ next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
+ my(@bindepth) = binDepths($ens);
+ for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
+ next unless ($bin+1>=$outGrid_firstBin && $bin+1<=$outGrid_lastBin);
+ next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
+ next unless numberp($LADCP{ENSEMBLE}[$ens]->{W12}[$bin]);
+ my($bi) = $bindepth[$bin]/$opt_o;
+ printf(GMT "%d %f %f $ens_width $bin_length\n",
+ $LADCP{ENSEMBLE}[$ens]->{NUMBER},
+ $bindepth[$bin],
+ $LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin] - $UPCAST{MEDIAN_W}[$bi]);
+ }
+ }
+
+ if ($P{water_depth} > 0) { # SEABED
+ GMT_psxy('-G204/153/102');
+ print(GMT "$xmin $ymax\n$xmax $ymax\n$xmax $P{water_depth}\n $xmin $P{water_depth}\n");
+ }
+
+ GMT_unitcoords(); # LABELS
+ GMT_pstext('-F+f9,Helvetica,orange+jTR -N -Gwhite');
+ print(GMT "0.99 0.01 V$VERSION\n");
+ GMT_pstext('-F+f14,Helvetica,blue+jTL -N');
+ print(GMT "0.01 -0.06 $P{out_basename} $P{run_label}\n");
+
+ my($depth_tics) = ($ymax-$ymin < 1000) ? 'f10a100' : 'f100a500'; # AXES
+ my($ens_tics) = ($xmax-$xmin < 4000) ? 'f50a500' : 'f500a2000';
+ GMT_setR($R);
+ GMT_psbasemap("-B$ens_tics:'Ensemble [#]':/$depth_tics:'Depth [m]':WeSn");
+
+ GMT_setAnnotFontSize(7); # SCALE BAR
+ GMT_psscale("-Dn0.83/0.1+w3/0.4+e $C -B/:'w<1,2>\@-residual\@-':");
+
+ GMT_end(); # FINISH PLOT
+}
+
+1; # return true on require
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/plot_residuals34.pl Fri Aug 05 11:02:51 2016 -0400
@@ -0,0 +1,87 @@
+#======================================================================
+# P L O T _ R E S I D U A L S 3 4 . P L
+# doc: Wed Jun 1 19:05:22 2016
+# dlm: Wed Jun 1 19:35:58 2016
+# (c) 2016 A.M. Thurnherr
+# uE-Info: 82 67 NIL 0 0 72 0 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+# Jun 1, 2016: - created from [plot_residuals.pl]
+
+require "$ANTS/libGMT.pl";
+
+sub plot_residuals34($)
+{
+ my($pfn) = @_;
+
+ return unless ($P{max_depth});
+
+ my($xmin) = $P{min_ens}-0.5;
+ my($xmax) = $P{max_ens}+0.5;
+ my($ymin) = round(antsParam('min_depth')-25,50);
+ my($ymax) = ($P{water_depth} > 0) ?
+ round($P{water_depth}+25,50) :
+ round($P{max_depth}+$P{ADCP_bin_length}+25,50);
+
+ my($ens_width) = 10 / ($P{max_ens} - $P{min_ens} + 1);
+ my($bin_length) = 10 * $P{ADCP_bin_length} /
+ ($P{max_depth}-$P{min_depth}+$P{ADCP_bin_length});
+
+ my($R) = "-R$xmin/$xmax/$ymin/$ymax";
+ GMT_begin($pfn,'-JX10/-10',$R,'-P');
+
+ my($C) = "-C$WCALC/residuals.cpt";
+ GMT_psxy("$C -Sr");
+ for ($ens=$firstGoodEns; $ens<$LADCP_atbottom; $ens++) { # downcast
+ next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
+ my(@bindepth) = binDepths($ens);
+ for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
+ next unless ($bin+1>=$outGrid_firstBin && $bin+1<=$outGrid_lastBin);
+ next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
+ next unless numberp($LADCP{ENSEMBLE}[$ens]->{W34}[$bin]);
+ my($bi) = $bindepth[$bin]/$opt_o;
+ printf(GMT "%d %f %f $ens_width $bin_length\n",
+ $LADCP{ENSEMBLE}[$ens]->{NUMBER},
+ $bindepth[$bin],
+ $LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin] - $DNCAST{MEDIAN_W}[$bi]);
+ }
+ }
+ for ($ens=$LADCP_atbottom; $ens<=$lastGoodEns; $ens++) { # upcast
+ next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
+ my(@bindepth) = binDepths($ens);
+ for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
+ next unless ($bin+1>=$outGrid_firstBin && $bin+1<=$outGrid_lastBin);
+ next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
+ next unless numberp($LADCP{ENSEMBLE}[$ens]->{W34}[$bin]);
+ my($bi) = $bindepth[$bin]/$opt_o;
+ printf(GMT "%d %f %f $ens_width $bin_length\n",
+ $LADCP{ENSEMBLE}[$ens]->{NUMBER},
+ $bindepth[$bin],
+ $LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin] - $UPCAST{MEDIAN_W}[$bi]);
+ }
+ }
+
+ if ($P{water_depth} > 0) { # SEABED
+ GMT_psxy('-G204/153/102');
+ print(GMT "$xmin $ymax\n$xmax $ymax\n$xmax $P{water_depth}\n $xmin $P{water_depth}\n");
+ }
+
+ GMT_unitcoords(); # LABELS
+ GMT_pstext('-F+f9,Helvetica,orange+jTR -N -Gwhite');
+ print(GMT "0.99 0.01 V$VERSION\n");
+ GMT_pstext('-F+f14,Helvetica,blue+jTL -N');
+ print(GMT "0.01 -0.06 $P{out_basename} $P{run_label}\n");
+
+ my($depth_tics) = ($ymax-$ymin < 1000) ? 'f10a100' : 'f100a500'; # AXES
+ my($ens_tics) = ($xmax-$xmin < 4000) ? 'f50a500' : 'f500a2000';
+ GMT_setR($R);
+ GMT_psbasemap("-B$ens_tics:'Ensemble [#]':/$depth_tics:'Depth [m]':WeSn");
+
+ GMT_setAnnotFontSize(7); # SCALE BAR
+ GMT_psscale("-Dn0.83/0.1+w3/0.4+e $C -B/:'w<3,4>\@-residual\@-':");
+
+ GMT_end(); # FINISH PLOT
+}
+
+1; # return true on require
--- a/plot_wprof.pl Wed May 25 12:14:29 2016 -0400
+++ b/plot_wprof.pl Fri Aug 05 11:02:51 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: Tue May 24 22:31:14 2016
+# dlm: Thu May 26 11:29:32 2016
# (c) 2015 A.M. Thurnherr
-# uE-Info: 19 51 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 20 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -17,6 +17,7 @@
# May 24, 2016: - BUG: ymin did not work for nsamp
# - fixed for partial-depth profiles
# - suppress plotting of nsamp == 0
+# May 26, 2016: - added instrument coord system to plot labels
# Tweakables:
#
@@ -139,7 +140,9 @@
GMT_pstext('-F+f9,Helvetica,CornFlowerBlue+jTL -N');
printf(GMT "0.64 1.020 $LADCP{BEAM_FREQUENCY}kHz $LADCP{INSTRUMENT_TYPE} $P{ADCP_orientation}\n");
- printf(GMT "0.64 1.055 bin setup\n 0.77 1.055 : %.1fm/%1.fm/%1.fm\n",
+# printf(GMT "0.64 1.055 %s\n 0.77 1.055 : %.1fm/%1.fm/%1.fm\n",
+ printf(GMT "0.64 1.055 %s [%.1fm/%1.fm/%1.fm]\n",
+ $LADCP{BEAM_COORDINATES} ? 'beam vels' : 'Earth vels',
$LADCP{BLANKING_DISTANCE},$LADCP{TRANSMITTED_PULSE_LENGTH},$LADCP{BIN_LENGTH});
print(GMT "0.64 1.090 mean tilt\n 0.77 1.096 :\n");
print(GMT "0.64 1.130 rms a\@-pkg\@-\n 0.77 1.1315 :\n");
--- a/version.pl Wed May 25 12:14:29 2016 -0400
+++ b/version.pl Fri Aug 05 11:02:51 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 May 25 12:12:54 2016
+# dlm: Thu Jun 2 12:20:26 2016
# (c) 2015 A.M. Thurnherr
-# uE-Info: 24 20 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 27 110 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -21,8 +21,8 @@
#$VERSION = '1.1'; # Jan 4, 2016
#$VERSION = '1.2'; # May 12, 2016
-$VERSION = '1.3beta1';
+$VERSION = '1.3beta2';
$antsMinLibVersion = 6.6;
-$ADCP_tools_minVersion = 1.6; # May 19, 2016 (RDI_Coords with bin interpolation)
+$ADCP_tools_minVersion = 1.7; # May 2016 (RDI_Coords with bin interpolation & better pitch/roll rotation)