--- a/HISTORY Mon Jan 04 11:19:09 2016 +0000
+++ b/HISTORY Tue Feb 02 14:55:24 2016 +0000
@@ -1,9 +1,9 @@
======================================================================
H I S T O R Y
doc: Mon Oct 12 16:09:24 2015
- dlm: Mon Jan 4 10:57:43 2016
+ dlm: Sun Jan 24 13:47:15 2016
(c) 2015 A.M. Thurnherr
- uE-Info: 30 59 NIL 0 0 72 3 2 4 NIL ofnI
+ uE-Info: 39 40 NIL 0 0 72 3 2 4 NIL ofnI
======================================================================
V1.0:
@@ -28,3 +28,12 @@
Jan 4, 2016:
- decreased default vertical resolution to 20m [defaults.pl] [LADCP_w_ocean]
- removed beta from version [version.pl] [.hg/hgrc]
+
+V1.2:
+ Jan 5, 2016: V1.2beta
+ - updated [version.pl]
+ - added support for [ADCP_tools_lib.pl]: [version.pl] [LADCP_w_ocean]
+ Jan 24, 2016:
+ - added QC for mean residuals
+ - added QC for dual-head wprofs
+ - added automatic editing of bad time lagged data
--- a/LADCP_VKE Mon Jan 04 11:19:09 2016 +0000
+++ b/LADCP_VKE Tue Feb 02 14:55:24 2016 +0000
@@ -2,9 +2,9 @@
#======================================================================
# L A D C P _ V K E
# doc: Tue Oct 14 11:05:16 2014
-# dlm: Tue Dec 29 11:04:23 2015
+# dlm: Mon Jan 25 09:31:12 2016
# (c) 2012 A.M. Thurnherr
-# uE-Info: 43 58 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 44 49 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
$antsSummary = 'calculate VKE from LADCP-derived vertical-velocity profiles';
@@ -40,7 +40,8 @@
# Nov 30, 2015: - BUG: -a was not allowed in usage
# - added -f
# Dec 27, 2015: - reduced minlim for eps to 1e-13 W/kg
-# Dec 29, 2015: - added 3rd consistency check (p0 limit)
+# Dec 29, 2015: - added 3rd consistency check (p0 limit)
+# Jan 25, 2016: - added software version %PARAM
($ANTSLIB) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
($WCALC) = ($0 =~ m{^(.*)/[^/]*$});
@@ -51,6 +52,7 @@
require "$ANTSLIB/ants.pl";
require "$ANTSLIB/libLADCP.pl";
require "$ANTSLIB/libGMT.pl";
+&antsAddParams('LADCP_VKE',"Version $VERSION");
use FileHandle;
use IPC::Open2;
--- a/LADCP_w_CTD Mon Jan 04 11:19:09 2016 +0000
+++ b/LADCP_w_CTD Tue Feb 02 14:55:24 2016 +0000
@@ -2,9 +2,9 @@
#======================================================================
# L A D C P _ W _ C T D
# doc: Mon Nov 3 17:34:19 2014
-# dlm: Tue Oct 13 10:51:14 2015
+# dlm: Mon Jan 25 09:32:32 2016
# (c) 2014 A.M. Thurnherr
-# uE-Info: 50 0 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 51 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
$antsSummary = 'pre-process SBE 9plus CTD data for LADCP_w';
@@ -47,6 +47,7 @@
# Oct 12, 2015: - require ANTSlibs V6.2 for release
# Oct 13, 2015: - added code to deal with CNV files with wrong number of scans in header
# Oct 13, 2015: - adapted to [version.pl]
+# Jan 25, 2016: - added software version %PARAM
($ANTS) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
($WCALC) = ($0 =~ m{^(.*)/[^/]*$});
@@ -58,6 +59,7 @@
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,
--- a/LADCP_w_ocean Mon Jan 04 11:19:09 2016 +0000
+++ b/LADCP_w_ocean Tue Feb 02 14:55:24 2016 +0000
@@ -2,9 +2,9 @@
#======================================================================
# L A D C P _ W _ O C E A N
# doc: Fri Dec 17 18:11:13 2010
-# dlm: Mon Jan 4 10:54:26 2016
+# dlm: Wed Jan 27 20:55:29 2016
# (c) 2010 A.M. Thurnherr
-# uE-Info: 196 64 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 450 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# TODO:
@@ -194,6 +194,15 @@
# 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
# HISTORY END
# CTD REQUIREMENTS
@@ -270,8 +279,8 @@
require "$WCALC/bottom_tracking.pl";
require "$ANTS/ants.pl";
require "$ANTS/libstats.pl";
-require "$ADCP_TOOLS/RDI_PD0_IO.pl";
-require "$ADCP_TOOLS/RDI_Coords.pl";
+require "$ADCP_TOOLS/ADCP_tools_lib.pl";
+&antsAddParams('LADCP_w_ocean',"Version $VERSION");
use IO::Handle;
@@ -290,7 +299,6 @@
'[-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]>]',
- '[offset CTD -d)epth <by[0m]>]',
'[-i)nitial CTD time offset <guestimate> [-u)se as final]]',
'[calculate -n) <lags,lags[10,100]>] [lag -w)indow <sz,sz[240s,20s]>] [lag-p)iece <CTD_elapsed_min|+[,...]>]',
'[require top-3) lags to account for <frac[0.6]> of all]',
@@ -328,23 +336,41 @@
#-----------------------------
require "$WCALC/defaults.pl"; # load default/option parameters
do "$processing_param_file"; # load processing parameters
-eval($opt_x) if defined($opt_x); # eval cmd-line expression to override anything
+
+$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);
-$RDI_Coords::minValidVels = 4 if ($opt_4); # decode cmd-line arguments
-($LADCP_firstBin,$LADCP_lastBin) = split(',',$opt_b);
+if (defined($opt_x)) { # eval cmd-line expression to override anything
+ $processing_options .= " -x $opt_x";
+ eval($opt_x);
+}
+
+if ($opt_4) { # disallow 3-beam solutions
+ $processing_options .= " -4";
+ $RDI_Coords::minValidVels = 4;
+}
+
+($LADCP_firstBin,$LADCP_lastBin) = split(',',$opt_b); # select valid bins
croak("$0: cannot decode -b $opt_b\n")
unless (numberp($LADCP_firstBin) &&
($LADCP_lastBin eq '*' || numberp($LADCP_lastBin)));
+$processing_options .= " -b $opt_b";
+
+$outGrid_firstBin = $LADCP_firstBin if ($outGrid_firstBin eq '*'); # default bins to use in gridded output
+$outGrid_lastBin = $LADCP_lastBin if ($outGrid_lastBin eq '*'); # NB: can still be '*'!!!
+
@number_of_timelag_windows = split(',',$opt_n);
croak("$0: cannot decode -n $opt_n\n")
unless numberp($number_of_timelag_windows[0]) && numberp($number_of_timelag_windows[1]);
+$processing_options .= " -n $opt_n";
@length_of_timelag_windows = split(',',$opt_w);
croak("$0: cannot decode -w $opt_w\n")
unless numberp($length_of_timelag_windows[0]) && numberp($length_of_timelag_windows[1]);
-
+$processing_options .= " -w $opt_w";
croak("$0: \$out_basename undefined\n") # all plotting routines use this
unless defined($out_basename);
+&antsAddParams('processing_options',$processing_options);
&antsAddParams('out_basename',$out_basename);
&antsAddParams('profile_id',$PROF,'run_label',$RUN);
@@ -419,15 +445,17 @@
unless ($refLr_lastBin>=1 && $refLr_lastBin<=$LADCP{N_BINS});
error("$0: first reference-layer bin > last reference-layer bin\n")
unless ($refLr_firstBin <= $refLr_lastBin);
+&antsAddParams('refLr_firstBin',$refLr_firstBin,'refLr_lastBin',$refLr_lastBin);
-$LADCP_lastBin = $LADCP{N_BINS}
- if ($LADCP_lastBin eq '*');
+$LADCP_lastBin = $LADCP{N_BINS} if ($LADCP_lastBin eq '*');
+$outGrid_lastBin = $LADCP{N_BINS} if ($outGrid_lastBin eq '*');
error("$0: first valid LADCP bin outside valid range\n")
unless ($LADCP_firstBin>=1 && $LADCP_firstBin<=$LADCP{N_BINS});
error("$0: last valid LADCP bin outside valid range\n")
unless ($LADCP_lastBin>=1 && $LADCP_lastBin<=$LADCP{N_BINS});
error("$0: first valid LADCP bin > last valid LADCP bin\n")
unless ($LADCP_firstBin <= $LADCP_lastBin);
+&antsAddParams('LADCP_firstBin',$LADCP_firstBin,'LADCP_lastBin',$LADCP_lastBin);
warning(0,"first reference-layer bin < first valid LADCP bin\n")
unless ($refLr_firstBin >= $LADCP_firstBin);
@@ -478,11 +506,18 @@
}
$nvv = $cte = 0;
+
+ if ($pitch_bias || $roll_bias || $heading_bias) {
+ &antsAddParams('pitch_bias',$pitch_bias,
+ 'roll_bias',$roll_bias,
+ 'heading_bias',$heading_bias);
+ for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
+ correctAttitude($ens,$pitch_bias,$roll_bias,$heading_bias);
+ }
+ }
for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
- correctAttitude($ens,$pitch_bias,$roll_bias,$heading_bias);
$nvv += countValidBeamVels($ens);
$cte += editCorr($ens,$opt_c);
-# $pte += editTilt($ens,$opt_t);
}
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));
@@ -495,7 +530,6 @@
for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
$nvv += countValidBeamVels($ens);
$cte += editCorr_Earthcoords($ens,$opt_c);
-# $pte += editTilt($ens,$opt_t);
}
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));
@@ -515,8 +549,10 @@
if ($LADCP{BEAM_COORDINATES}) {
my($dummy);
progress("Calculating earth-coordinate velocities...\n");
- progress("\tdiscarding velocities from beam $bad_beam\n")
- if ($bad_beam);
+ if ($bad_beam) {
+ progress("\tdiscarding velocities from beam $bad_beam\n");
+ &antsAddParams('bad_beam_discarded',$bad_beam);
+ }
$nvw = 0;
for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
@@ -586,11 +622,14 @@
progress("Editing earth-coordinate velocity data...\n");
+&antsAddParams('per_ens_outliers_mad_limit',$per_ens_outliers_mad_limit)
+&antsAddParams('farthest_valid_bins_truncated',$truncate_farthest_valid_bins)
+ if ($truncate_farthest_valid_bins);
$evrm = $trrm = $worm = 0;
for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
$evrm += editErrVel($ens,$opt_e);
$worm += editWOutliers($ens,$per_ens_outliers_mad_limit);
- $trrm += editTruncRange($ens,$truncate_farthest_valid_bins) if ($truncate_farthest_valid_bins > 0);
+ $trrm += editTruncRange($ens,$truncate_farthest_valid_bins) if ($truncate_farthest_valid_bins);
}
progress("\terror-velocity threshold (-e %.2f m/s): %d velocites removed (%d%% of total in bins $LADCP_firstBin-$LADCP_lastBin)\n",
$opt_e,$evrm,round(100*$evrm/$nvw));
@@ -676,6 +715,7 @@
#----------------------------------------------------------------------
progress("Editing additional earth-coordinate velocity data...\n");
+&antsAddParams('per_bin_valid_frac_lim',$per_bin_valid_frac_lim);
my($first_bad_bin);
for (my($bin)=$LADCP_firstBin-1; $bin<$LADCP_lastBin-1; $bin++) {
@@ -738,14 +778,6 @@
error("$0: CTD start depth must be numeric\n")
unless numberp($CTD{DEPTH}[0]);
-if (($CTD{DEPTH}[0] < -$opt_d) && !defined($opt_d)) {
- $opt_d = -1 * round($CTD{DEPTH}[0]);
-}
-if ($opt_d > 0) {
- progress("\tadding $opt_d dbar offset to pressure data\n");
- for (my($i)=0; $i<@{$CTD{DEPTH}}; $i++) { $CTD{DEPTH}[$i] += $opt_d; }
- $CTD_maxdepth += $opt_d;
-}
progress("\tstart depth = %.1fm\n",$CTD{DEPTH}[0]);
progress("\tmax depth = %dm (# $CTD_atbottom)\n",$CTD_maxdepth);
@@ -880,6 +912,8 @@
progress("Merging CTD with LADCP data...\n");
+&antsAddParams('w_max_lim',$w_max_lim);
+
my($cli) = 0; # current-lag index
my($lag) = $CTD_time_lag[$cli]; # current lag
@@ -946,10 +980,11 @@
&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(1,"%.2f cm/s (large) reference-layer w_ocean away from boundaries\n",100*sqrt($sumWsqI/$nWsqI));
- } elsif (sqrt($sumWsqI/$nWsqI) > 0.1) {
- warning(2,"%.2f cm/s (very large) reference-layer w_ocean away from boundaries\n",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) {
+ warning(2,"%.2f cm/s (implausibly large) reference-layer w_ocean away from boundaries\n",100*sqrt($sumWsqI/$nWsqI));
}
} elsif ($nWsq > 0) {
&antsAddParams('rms_w_reflr_err',sqrt($sumWsq/$nWsq),'rms_w_reflr_err_interior',nan);
@@ -965,8 +1000,10 @@
progress("Calculating volume-scattering coefficients (Sv)...\n");
calc_backscatter_profs($firstGoodEns,$lastGoodEns);
-progress("\tapplying empirical Sv correction\n");
-correct_backscatter($firstGoodEns,$lastGoodEns);
+if (@nSv) {
+ progress("\tapplying empirical Sv correction\n");
+ correct_backscatter($firstGoodEns,$lastGoodEns);
+}
#----------------------------------------------------------------------------
# Find Seabed & Edit data
@@ -999,7 +1036,7 @@
find_seabed(\%LADCP,$LADCP_atbottom,$LADCP{BEAM_COORDINATES});
if (defined($water_depth) && defined($water_depth_BT)) {
my($dd) = abs($water_depth_BT - $water_depth);
- warning(2,sprintf("Large instrument vs. backscatter-derived water-depth difference (%.1fm)\n",$dd))
+ 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)) {
@@ -1007,9 +1044,12 @@
$SS_use_BT = 1;
}
if ($SS_use_BT) {
+ &antsAddParams('water_depth_from','BT_data');
$water_depth = $water_depth_BT;
$sig_water_depth = $sig_water_depth_BT;
- }
+ } else {
+ &antsAddParams('water_depth_from','echo_amplitudes');
+ }
}
if (defined($water_depth)) {
@@ -1027,7 +1067,6 @@
&antsAddParams('water_depth','unknown','water_depth.sig','nan');
}
-
if ($LADCP{ENSEMBLE}[$LADCP_atbottom]->{XDUCER_FACING_DOWN}) { # DOWNLOOKER
&antsAddParams('ADCP_orientation','downlooker');
@@ -1040,9 +1079,15 @@
progress("Editing data to remove sidelobe interference from sea surface...\n");
($nvrm,$nerm) = editSideLobes($firstGoodEns,$lastGoodEns,undef);
progress("\t$nvrm velocities from $nerm ensembles removed\n");
+ &antsAddParams('sidelobe_editing','surface+seabed','vessel_draft',$vessel_draft);
+ } else {
+ &antsAddParams('sidelobe_editing','seabed');
}
if ($PPI_editing) {
+ &antsAddParams('PPI_editing','seabed');
+ &antsAddParams('PPI_extend_upper_limit',$PPI_extend_upper_limit)
+ if numberp($PPI_extend_upper_limit);
progress("Editing data to remove PPI from seabed...\n");
progress("\tConstructing depth-average soundspeed profile...\n");
die("assertion failed") unless defined($water_depth);
@@ -1075,9 +1120,14 @@
} else {
warning(2,"unknown water depth --- cannot edit UL data for sidelobe interference from seabed\n");
}
- }
-
+ &antsAddParams('sidelobe_editing','surface+seabed','vessel_draft',$vessel_draft);
+ } else {
+ &antsAddParams('sidelobe_editing','surface','vessel_draft',$vessel_draft);
+ }
if ($PPI_editing) {
+ &antsAddParams('PPI_editing','surface');
+ &antsAddParams('PPI_extend_upper_limit',$PPI_extend_upper_limit)
+ if numberp($PPI_extend_upper_limit);
progress("Editing data to remove PPI from sea surface...\n");
progress("\tConstructing depth-average soundspeed profile...\n");
$DASSprof[0] = my($sum) = 0;
@@ -1098,6 +1148,7 @@
#----------------------------------------------------------------------
progress("Removing data from instrument at surface...\n");
+&antsAddParams('surface_layer_depth',$surface_layer_depth);
$nerm = editSurfLayer($firstGoodEns,$lastGoodEns,$surface_layer_depth);
progress("\t$nerm ensembles removed\n");
@@ -1347,9 +1398,9 @@
} # unless ($opt_q);
-#-------------------------
-# Output Per-Bin Resiudals
-#-------------------------
+#--------------------------------------------
+# Calculate and Output Bin-Averaged Resiudals
+#--------------------------------------------
if (@out_BR) {
progress("Binning residuals...");
@@ -1375,14 +1426,49 @@
- $UPCAST{MEDIAN_W}[$bindepth[$bin]/$opt_o]);
}
}
- local(@dc_avg_bres,@uc_avg_bres);
- for (my($bin)=0; $bin<max(scalar(@dc_bres),scalar(@uc_bres)); $bin++) { # calc means/stddevs
- $dc_avg_bres[$bin] = avg(@{$dc_bres[$bin]});
- $uc_avg_bres[$bin] = avg(@{$uc_bres[$bin]});
- $dc_sig_bres[$bin] = stddev2($dc_avg_bres[$bin],@{$dc_bres[$bin]}),scalar(@{$dc_bres[$bin]});
- $uc_sig_bres[$bin] = stddev2($uc_avg_bres[$bin],@{$uc_bres[$bin]}),scalar(@{$uc_bres[$bin]});
+ local(@dc_avg_bres,@uc_avg_bres,@dc_avg_bres_nsamp,@uc_avg_bres_nsamp); # calc means/stddevs
+ for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
+ $dc_bres_nsamp[$bin] = @{$dc_bres[$bin]};
+ if ($dc_bres_nsamp[$bin] > 0) {
+ $dc_avg_bres[$bin] = avg(@{$dc_bres[$bin]});
+ $dc_sig_bres[$bin] = stddev2($dc_avg_bres[$bin],@{$dc_bres[$bin]});
+ } else {
+ $dc_avg_bres[$bin] = nan;
+ $dc_sig_bres[$bin] = nan;
+ }
+ $uc_bres_nsamp[$bin] = @{$uc_bres[$bin]};
+ if ($uc_bres_nsamp[$bin] > 0) {
+ $uc_avg_bres[$bin] = avg(@{$uc_bres[$bin]});
+ $uc_sig_bres[$bin] = stddev2($uc_avg_bres[$bin],@{$uc_bres[$bin]});
+ } else {
+ $uc_avg_bres[$bin] = nan;
+ $uc_sig_bres[$bin] = nan;
+ }
}
+ local($uc_bres_max_nsamp,$dc_bres_max_nsamp) = (0,0); # calc rms in block of well-determined bins
+ for (my($bin)=$LADCP_firstBin; $bin<=$LADCP_lastBin-1; $bin++) { # SKIP 1ST BIN!!!
+ next if ($bin+1<$outGrid_firstBin || $bin+1>$outGrid_lastBin); # skip bins not included in gridded output
+ $dc_bres_max_nsamp = $dc_bres_nsamp[$bin] # nsamp in best sampled bin
+ if ($dc_bres_nsamp[$bin] > $dc_bres_max_nsamp);
+ $uc_bres_max_nsamp = $uc_bres_nsamp[$bin]
+ if ($uc_bres_nsamp[$bin] > $uc_bres_max_nsamp);
+ }
+ my($dc_sumsq,$uc_sumsq,$dc_n,$uc_n) = (0,0,0,0); # calc rms residual
+ for (my($bin)=$LADCP_firstBin; $bin<=$LADCP_lastBin-1; $bin++) { # SKIP 1ST BIN
+ next if ($bin+1<$outGrid_firstBin || $bin+1>$outGrid_lastBin); # skip bins not included in gridded output
+ if ($dc_bres_nsamp[$bin] >= $dc_bres_max_nsamp/3) { # skip bins with < 1/3 max(nsamp)
+ $dc_sumsq += $dc_avg_bres[$bin]**2;
+ $dc_n++;
+ }
+ if ($uc_bres_nsamp[$bin] >= $uc_bres_max_nsamp/3) {
+ $uc_sumsq += $uc_avg_bres[$bin]**2;
+ $uc_n++;
+ }
+ }
+ local($dc_bres_rms) = sqrt($dc_sumsq/$dc_n);
+ local($uc_bres_rms) = sqrt($uc_sumsq/$uc_n);
+
progress("\nWriting binned residuals to ");
my($saveParams) = $antsCurParams;
@@ -1402,8 +1488,8 @@
open(STDOUT,$of) || error("$of: $!\n");
undef($antsActiveHeader) unless ($ANTS_TOOLS_AVAILABLE);
for (my($bin)=0; $bin<max(scalar(@dc_bres),scalar(@uc_bres)); $bin++) {
- &antsOut($bin+1,$dc_avg_bres[$bin],$dc_sig_bres[$bin],scalar(@{$dc_bres[$bin]}),
- $uc_avg_bres[$bin],$uc_sig_bres[$bin],scalar(@{$uc_bres[$bin]}));
+ &antsOut($bin+1,$dc_avg_bres[$bin],$dc_sig_bres[$bin],$dc_bres_nsamp[$bin],
+ $uc_avg_bres[$bin],$uc_sig_bres[$bin],$uc_bres_nsamp[$bin]);
}
&antsOut('EOF'); open(STDOUT,">&2");
}
--- a/LADCP_w_postproc Mon Jan 04 11:19:09 2016 +0000
+++ b/LADCP_w_postproc Tue Feb 02 14:55:24 2016 +0000
@@ -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 Oct 13 10:51:31 2015
+# dlm: Sat Jan 30 14:50:53 2016
# (c) 2015 A.M. Thurnherr
-# uE-Info: 50 0 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 640 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
$antsSummary = 'edit and re-grid LADCP vertical-velocity samples';
@@ -47,6 +47,14 @@
# - added run-label(s) to plot
# - require ANTSlibs V6.2 for release
# Oct 13, 2015: - adapted to [version.pl]
+# Jan 20, 2016: - added dc_w.diff, uc_w.diff for dual-head profiles
+# Jan 21, 2016: - added %_wcorr.ps plot output
+# - added correlation-based QC to wprof plot
+# Jan 22, 2016: - added output -d)ir option
+# Jan 24, 2016: - BUG: uc_/dc_corr returned single nan on failure
+# Jan 25, 2016: - added software version %PARAM
+# Jan 26, 2016: - added -v
+# Jan 30, 2016: - added -w
($ANTS) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
($WCALC) = ($0 =~ m{^(.*)/[^/]*$});
@@ -59,23 +67,62 @@
require "$ANTS/ants.pl";
require "$ANTS/libstats.pl";
require "$ANTS/libGMT.pl";
+&antsAddParams('LADCP_w_postproc',"Version $VERSION");
-&antsUsage('b:i:k:o:p:',1,
+&antsUsage('b:d:i:k:l:o:p:v:w:',1,
'[profile -i)d <id>]',
'[-o)utput bin <resolution>] [-k) require <min> samples]',
+ '[-v)alid bins <DL first>,<DL last>[,<UL first>,<UL_last>]',
+ '[-w) <DL_dc_field>,<DL_uc_field>[,<UL_dc_field>,<UL_uc_field>]',
+ '[-s)urface-layer <limit[300m]>]',
+ '[ouptput -d)ir <name>]',
'[-p)lot <[%03d_wprof.eps]> [-b)t <wprof>]]',
- '<.wsamp file> [.wsamp file]');
+ '<DL.wsamp file> [UL.wsamp file] (or only <UL.wsamp file>)');
+
+$dual_head = (@ARGV==1); # single or dual head
+
+$id = defined($opt_i) ? $opt_i : &antsParam('profile_id'); # ensure profile id
+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");
+}
+
+&antsCardOpt(\$opt_s,300); # surface layer depth limit
-$dual_head = (@ARGV==1);
+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)));
+ $fvBin = &antsRequireParam('LADCP_firstBin') if ($fvBin eq '*'); # corresponding UL values set below
+ $lvBin = &antsRequireParam('LADCP_lastBin') if ($lvBin eq '*');
+ &antsAddParams('DL_first_valid_bin',$fvBin,
+ 'DL_last_valid_bin', $lvBin);
+} else {
+ &antsAddParams('DL_first_valid_bin',&antsRequireParam('outgrid_firstbin'),
+ 'DL_last_valid_bin',&antsRequireParam('outgrid_lastbin'));
+}
-if (defined($opt_o)) {
+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,
+ 'UL_dc_w_field',$Udwf,'UL_uc_w_field',$Uuwf);
+} else {
+ ($Ddwf,$Duwf,$Udwf,$Uuwf) = ('w','w','w','w');
+}
+
+if (defined($opt_o)) { # output grid resolution
$opt_o_override = 1;
&antsCardOpt($opt_o);
} else {
$opt_o = &antsRequireParam('outgrid_dz');
}
-if (defined($opt_k)) {
+if (defined($opt_k)) { # minimum number of required samples
$opt_k_override = 1;
&antsCardOpt($opt_k);
} else {
@@ -84,20 +131,19 @@
}
#----------------------------------------------------------------------
-# Redirect STDOUT to %.wprof & create %_wprof.ps if STDOUT is a tty
+# Redirect STDOUT to %.wprof & create plots if STDOUT is a tty
#----------------------------------------------------------------------
-$id = defined($opt_i) ? $opt_i : &antsParam('profile_id');
-croak("$0: no profile_id in first file => -i required\n")
- unless defined($id);
-
if (-t STDOUT) {
- $opt_p = '%03d_wprof.ps' unless defined($opt_p);
- $outfile = sprintf('%03d.wprof',$id);
+ $opt_p = defined($opt_d) ? "$opt_d/%03d_wprof.ps,$opt_d/%03d_wcorr.ps"
+ : '%03d_wprof.ps,%03d_wcorr.ps'
+ unless defined($opt_p);
+ $outfile = defined($opt_d) ? sprintf('%s/%03d.wprof',$opt_d,$id)
+ : sprintf('%03d.wprof',$id);
open(STDOUT,">$outfile") || die("$outfile: $!\n");
}
-croak("$0: -b only makes sense when a plot is produced\n")
+croak("$0: -b only makes sense when plots are produced\n")
if $opt_b && !defined($opt_p);
#----------------------------------------------------------------------
@@ -174,10 +220,72 @@
}
#----------------------------------------------------------------------
-# Edit Data
+# Main Program
#----------------------------------------------------------------------
-$wF = &fnr('w');
+sub dc_corr()
+{
+ my($n) = 0;
+ my($ax,$ay) = (0,0);
+
+ for (my($bi)=int($opt_s/$opt_o); $bi<@DL_dc_median; $bi++) {
+ next unless numberp($DL_dc_median[$bi]) && numberp($UL_dc_median[$bi]);
+ $n++;
+ $ax += $DL_dc_median[$bi];
+ $ay += $UL_dc_median[$bi];
+ }
+ return (nan,nan,nan) unless ($n > 2);
+ $ax /= $n;
+ $ay /= $n;
+
+ my($syy,$sxy,$sxx) = (0,0,0);
+ for (my($bi)=int($opt_s/$opt_o); $bi<@DL_dc_median; $bi++) {
+ next unless numberp($DL_dc_median[$bi]) && numberp($UL_dc_median[$bi]);
+ my($xt) = $DL_dc_median[$bi] - $ax;
+ my($yt) = $UL_dc_median[$bi] - $ay;
+ $sxx += $xt * $xt;
+ $syy += $yt * $yt;
+ $sxy += $xt * $yt;
+ }
+ my($R) = $sxy/(sqrt($sxx * $syy) + 1e-16);
+ my($esig) = sqrt(($R**2)*$sxx/$n);
+ my($rsig) = sqrt((1-$R**2)*$sxx/$n);
+ return ($R,$esig,$rsig);
+}
+
+sub uc_corr(@)
+{
+ my($n) = 0;
+ my($ax,$ay) = (0,0);
+
+ for (my($bi)=int($opt_s/$opt_o); $bi<@DL_dc_median; $bi++) {
+ next unless numberp($DL_uc_median[$bi]) && numberp($UL_uc_median[$bi]);
+ $n++;
+ $ax += $DL_uc_median[$bi];
+ $ay += $UL_uc_median[$bi];
+ }
+ return (nan,nan,nan) unless ($n > 2);
+ $ax /= $n;
+ $ay /= $n;
+
+ my($syy,$sxy,$sxx) = (0,0,0);
+ for (my($bi)=int($opt_s/$opt_o); $bi<@DL_dc_median; $bi++) {
+ next unless numberp($DL_uc_median[$bi]) && numberp($UL_uc_median[$bi]);
+ my($xt) = $DL_uc_median[$bi] - $ax;
+ my($yt) = $UL_uc_median[$bi] - $ay;
+ $sxx += $xt * $xt;
+ $syy += $yt * $yt;
+ $sxy += $xt * $yt;
+ }
+ my($R) = $sxy/(sqrt($sxx * $syy) + 1e-16);
+ my($esig) = sqrt(($R**2)*$sxx/$n);
+ my($rsig) = sqrt((1-$R**2)*$sxx/$n);
+ return ($R,$esig,$rsig);
+}
+
+#----------------------------------------------------------------------
+
+$dcwF = &fnr($Ddwf); $ucwF = &fnr($Duwf);
$eF = &fnr('elapsed');
$dF = &fnr('depth');
$bF = &fnr('bin');
@@ -191,8 +299,14 @@
croak("$0: cannot determine day number\n")
unless defined($dayNoP);
+if (defined($opt_p)) {
+ ($sumPF,$corrPF) = split(/,/,$opt_p);
+ croak("$0: cannot decode -p $opt_p\n")
+ unless (length($corrPF)>0);
+}
+
my($R,$R2);
-if ($opt_p) { # begin plot
+if (defined($opt_p)) { # begin summary plot
$xmin = -0.1; $x2min = -200;
$xmax = 0.35; $x2max = 200;
$ymin = 0;
@@ -201,7 +315,7 @@
$plotsize = 13;
$R = "-R$xmin/$xmax/$ymin/$ymax";
$R2 = "-R$x2min/$x2max/$ymin/$ymax";
- GMT_begin(sprintf($opt_p,$id),"-JX$plotsize/-$plotsize",$R,'-P -X6 -Y4');
+ GMT_begin(sprintf($sumPF,$id),"-JX$plotsize/-$plotsize",$R,'-P -X6 -Y4');
GMT_psxy('-W1');
print(GMT "0 $ymin\n0 $ymax");
GMT_psxy('-L -G200');
@@ -237,19 +351,36 @@
$min_depth = 9e99; # sentinels
$max_depth = -9e99;
-$curF = ''; # current input file (sentinel)
+$curF = $P{PATHNAME}; # current input file (sentinel)
$filt = 0;
-for (my($curF),$r=0; &antsIn(); $r++) {
- if ($P{PATHNAME} ne $curF) { # new file
+for ($r=0; &antsIn(); $r++) {
+ if ($P{PATHNAME} ne $curF) { # 2nd file (UL data)
$curF = $P{PATHNAME};
- &antsInfo("WARNING: inconsistent %%outgrid_dz in input files")
+ $dcwF = &fnr($Udwf); $ucwF = &fnr($Uuwf);
+
+ &antsInfo("WARNING: inconsistent %%outgrid_dz in input files") # consistency checks
if (defined($P{outgrid_dz}) && $P{outgrid_dz}!=$opt_o &&!$opt_o_override);
&antsInfo("WARNING: inconsistent %%outgrid_minsamp in input files")
if defined($P{outgrid_minsamp}) && !$opt_k_override &&
(( $dual_head && $P{outgrid_minsamp}*2!=$opt_k) ||
(!$dual_head && $P{outgrid_minsamp}!=$opt_k));
+ if (defined($opt_v)) {
+ $fvBin = &antsRequireParam('LADCP_firstBin') if ($UL_fvBin eq '*'); # valid bin ranges
+ $lvBin = &antsRequireParam('LADCP_lastBin') if ($UL_lvBin eq '*');
+ &antsAddParams('UL_first_valid_bin',$fvBin,
+ 'UL_last_valid_bin', $lvBin);
+ } else {
+ &antsAddParams('UL_first_valid_bin',&antsRequireParam('outgrid_firstbin'),
+ 'UL_last_valid_bin',&antsRequireParam('outgrid_lastbin'));
+ }
+
+ for (my($bi)=0; $bi<=$#dcw1; $bi++) { # calc DL median profile
+ $DL_dc_median[$bi] = median(@{$dcw1[$bi]});
+ $DL_uc_median[$bi] = median(@{$ucw1[$bi]});
+ }
+
$bin_length = sprintf('%d & %d',round($bin_length),($P{ADCP_bin_length}))
unless (round($bin_length) == round($P{ADCP_bin_length}));
$pulse_length = sprintf('%d & %d',round($pulse_length),round($P{ADCP_pulse_length}))
@@ -259,48 +390,73 @@
$PROF = $STN = $ID = $id; $RUN = antsRequireParam('run_label');
undef(@rngMin); undef(@rngMax); undef(@bins);
- unless ($return = do "./EditParams") { # man perlfunc
+ unless ($return = do "./EditParams") { # man perlfunc
croak("./EditParams: $@\n") if ($@);
-# croak("do ./EditParams: $!\n") unless defined($return); # happens when EditParams does not exist
-# croak("couldn't run ./EditParams\n") unless $return; # happens when EditParams does not exist
}
- if ($opt_p && @dcw1) { # 2nd file in dual-head profile => plot 1st
+ if (defined($opt_p)) { # 2nd file in dual-head profile => plot 1st
GMT_psxy('-W2,255/127/80,-');
for (my($bi)=0; $bi<=$#dcw1; $bi++) {
- printf(GMT "%f %f\n",median(@{$dcw1[$bi]}),($bi+0.5)*$opt_o);
+ printf(GMT "%f %f\n",$DL_dc_median[$bi],($bi+0.5)*$opt_o);
}
GMT_psxy('-W2,46/139/87,-');
for (my($bi)=0; $bi<=$#ucw1; $bi++) {
- printf(GMT "%f %f\n",median(@{$ucw1[$bi]}),($bi+0.5)*$opt_o);
+ printf(GMT "%f %f\n",$DL_uc_median[$bi],($bi+0.5)*$opt_o);
}
undef(@dcw1); undef(@ucw1);
}
+ } # of 2nd file started
+
+ if (defined($opt_v)) { # explicit ranges of validity given
+ next if ($ants_[0][$bF]<$fvBin-1 || # => apply them
+ $ants_[0][$bF]>$lvBin-1);
+ } else { # no range of valid bins given
+ next if ($ants_[0][$bF]<$P{outgrid_firstbin}-1 || # => use values from [LADCP_w_ocean]
+ $ants_[0][$bF]>$P{outgrid_lastbin}-1);
}
- next if ($ants_[0][$bF]<$P{outgrid_firstbin}-1 || $ants_[0][$bF]>$P{outgrid_lastbin}-1);
- $filt++,next if &isBad();
- $min_depth = $ants_[0][$dF] if ($ants_[0][$dF] < $min_depth);
+
+ $filt++,next if &isBad(); # additional editing
+
+ $min_depth = $ants_[0][$dF] if ($ants_[0][$dF] < $min_depth); # update depth limits
$max_depth = $ants_[0][$dF] if ($ants_[0][$dF] > $max_depth);
+
my($bi) = $ants_[0][$dF]/$opt_o;
- if ($ants_[0][$dcF]) { # downcast
- push(@{$dcw[$bi]},$ants_[0][$wF]); # vertical velocity
- push(@{$dcw1[$bi]},$ants_[0][$wF]) if ($dual_head);
- push(@{$dce[$bi]},$ants_[0][$eF]); # elapsed time
- } else { # upcast
- push(@{$ucw[$bi]},$ants_[0][$wF]);
- push(@{$ucw1[$bi]},$ants_[0][$wF]) if ($dual_head);
+ if ($ants_[0][$dcF]) { # downcast
+ push(@{$dcw[$bi]},$ants_[0][$dcwF]); # vertical velocity
+ push(@{$dcw1[$bi]},$ants_[0][$dcwF]) if ($dual_head); # single-instrument w
+ push(@{$dce[$bi]},$ants_[0][$eF]); # elapsed time
+ } else { # upcast
+ push(@{$ucw[$bi]},$ants_[0][$ucwF]);
+ push(@{$ucw1[$bi]},$ants_[0][$ucwF]) if ($dual_head);
push(@{$uce[$bi]},$ants_[0][$eF]);
}
-}
+} # file-read loop
+
+if ($dual_head) {
+ for (my($bi)=0; $bi<=$#dcw1; $bi++) { # calc UL median & difference profiles
+ $UL_dc_median[$bi] = median(@{$dcw1[$bi]});
+ $UL_uc_median[$bi] = median(@{$ucw1[$bi]});
+ $dc_diff[$bi] = numberp($DL_dc_median[$bi]) && numberp($UL_dc_median[$bi])
+ ? $DL_dc_median[$bi] - $UL_dc_median[$bi] : nan;
+ $uc_diff[$bi] = numberp($DL_uc_median[$bi]) && numberp($UL_uc_median[$bi])
+ ? $DL_uc_median[$bi] - $UL_uc_median[$bi] : nan;
+ }
-if ($opt_p && @dcw1) { # 2nd file in dual-head profile in buffer => plot it
- GMT_psxy('-W2,255/127/80,.');
- for (my($bi)=0; $bi<=$#dcw1; $bi++) {
- printf(GMT "%f %f\n",median(@{$dcw1[$bi]}),($bi+0.5)*$opt_o);
- }
- GMT_psxy('-W2,46/139/87,.');
- for (my($bi)=0; $bi<=$#ucw1; $bi++) {
- printf(GMT "%f %f\n",median(@{$ucw1[$bi]}),($bi+0.5)*$opt_o);
+ ($dc_R,$dc_esig,$dc_rsig) = &dc_corr(); # correlation statistics
+ ($uc_R,$uc_esig,$uc_rsig) = &uc_corr();
+ &antsAddParams('dc_R',$dc_R,'uc_R',$uc_R,
+ 'dc_explained_stddev',$dc_esig,'uc_explained_stddev',$uc_esig,
+ 'dc_residual_stddev',$dc_rsig,'uc_residual_stddev',$uc_rsig);
+
+ if (defined($opt_p)) { # plot 2nd-instrument profiles
+ GMT_psxy('-W2,255/127/80,.');
+ for (my($bi)=0; $bi<=$#dcw1; $bi++) {
+ printf(GMT "%f %f\n",$UL_dc_median[$bi],($bi+0.5)*$opt_o);
+ }
+ GMT_psxy('-W2,46/139/87,.');
+ for (my($bi)=0; $bi<=$#ucw1; $bi++) {
+ printf(GMT "%f %f\n",$UL_uc_median[$bi],($bi+0.5)*$opt_o);
+ }
}
}
@@ -308,23 +464,24 @@
if ($filt > 0);
#----------------------------------------------------------------------
-# Average and Output Profile, Add to Plot
+# Average and Output Profiles, Add to Summary Plot
#----------------------------------------------------------------------
-@antsNewLayout = ('depth','hab',
- 'dc_elapsed','dc_w','dc_w.mad','dc_w.nsamp',
- 'uc_elapsed','uc_w','uc_w.mad','uc_w.nsamp');
-
if ($dual_head) { # selected %PARAMs only
+ @antsNewLayout = ('depth','hab',
+ 'dc_elapsed','dc_w','dc_w.mad','dc_w.nsamp',
+ 'uc_elapsed','uc_w','uc_w.mad','uc_w.nsamp',
+ 'dc_w.diff','uc_w.diff');
&antsAddParams('profile_id',$id,'lat',$P{lat},'lon',$P{lon});
&antsAddParams($dayNoP,$dn,'run_label',"$first_label & $P{run_label}");
&antsAddParams('outgrid_dz',$opt_o,'outgrid_minsamp',$opt_k);
- &antsAddParams('outgrid_firstbin',$P{outgrid_firstbin},'outgrid_lastbin',$P{outgrid_lastbin});
- &antsAddParams('ADCP_bin_length',round($bin_length),'ADCP_pulse_length',round($pulse_length))
- &antsAddParams('ADCP_blanking_distance',round($blanking_dist));
&antsAddParams('min_depth',round($min_depth),'max_depth',round($max_depth));
&antsAddParams('water_depth',$P{water_depth},'water_depth.sig',$P{water_depth.sig});
undef($antsOldHeaders);
+} else {
+ @antsNewLayout = ('depth','hab',
+ 'dc_elapsed','dc_w','dc_w.mad','dc_w.nsamp',
+ 'uc_elapsed','uc_w','uc_w.mad','uc_w.nsamp');
}
&antsInfo("WARNING: unknown water depth (no height-above-bottom)")
@@ -347,10 +504,35 @@
avg(@{$uce[$bi]}),
(($ucns[$bi]>=$opt_k)?$ucwm[$bi]:nan),(($ucns[$bi]>=$opt_k)?$ucwmad[$bi]:nan),
scalar(@{$ucw[$bi]}));
+ push(@{$out[$bi]},$dc_diff[$bi],$uc_diff[$bi])
+ if ($dual_head);
&antsOut(@{$out[$bi]});
}
-if ($opt_p) { # complete plot
+if (defined($opt_p)) { # complete summary plot
+ GMT_psxy('-W2/100/100/255');
+ print(GMT "-0.1 $opt_s\n0.07 $opt_s\n");
+ if ($dc_R < 0.3 || !numberp($dc_R)) {
+ &antsInfo("WARNING: low dc correlation (r = %.1f) between UL and DL data",$dc_R);
+ GMT_pstext('-Gwhite -Wred');
+ } elsif ($dc_R < 0.5) { GMT_pstext('-Gblack -Wyellow'); }
+ else { GMT_pstext('-Gblack -Wgreen'); }
+ printf(GMT "%f %f 12 0 0 BL %.1f\n",-0.07,0.9*$ymax,$dc_R);
+ GMT_pstext('-G255/127/80');
+ printf(GMT "%f %f 12 0 0 BL dc\n",-0.095,0.9*$ymax);
+ printf(GMT "%f %f 12 0 0 BL [%.1f/%.1f cm/s @~s@~\@-e/r\@-]\n",
+ 0.02,0.9*$ymax,100*$dc_esig,100*$dc_rsig);
+ if ($uc_R < 0.3 || !numberp($uc_R)) {
+ &antsInfo("WARNING: low uc correlation (r = %.1f) between UL and DL data",$uc_R);
+ GMT_pstext('-Gwhite -Wred');
+ } elsif ($uc_R < 0.5) { GMT_pstext('-Gblack -Wyellow'); }
+ else { GMT_pstext('-Gblack -Wgreen'); }
+ printf(GMT "%f %f 12 0 0 BL %.1f\n",-0.07,0.95*$ymax,$uc_R);
+ GMT_pstext('-G46/139/87');
+ printf(GMT "%f %f 12 0 0 BL uc\n",-0.095,0.95*$ymax);
+ printf(GMT "%f %f 12 0 0 BL [%.1f/%.1f cm/s @~s@~\@-e/r\@-]\n",
+ 0.02,0.95*$ymax,100*$uc_esig,100*$uc_rsig);
+
GMT_setR($R);
GMT_psxy('-W4,255/127/80'); # median profiles
@@ -397,6 +579,70 @@
GMT_pstext();
print(GMT '0.62 0.98 12 0 0 MR m.a.d.');
GMT_end();
+
+ if ($dual_head) { # correlation plot
+ my($mwm) = 0.05; # max(|w|) for axes
+ for (my($bi)=0; $bi<@DL_dc_median; $bi++) {
+ next unless numberp($DL_dc_median[$bi]) && numberp($UL_dc_median[$bi]);
+ $mwm = abs($DL_dc_median[$bi]) if abs($DL_dc_median[$bi]) > $mwm;
+ $mwm = abs($UL_dc_median[$bi]) if abs($UL_dc_median[$bi]) > $mwm;
+ }
+ for (my($bi)=0; $bi<@DL_uc_median; $bi++) {
+ next unless numberp($DL_uc_median[$bi]) && numberp($UL_uc_median[$bi]);
+ $mwm = abs($DL_uc_median[$bi]) if abs($DL_uc_median[$bi]) > $mwm;
+ $mwm = abs($UL_uc_median[$bi]) if abs($UL_uc_median[$bi]) > $mwm;
+ }
+ $mwm = int(100*$mwm+0.9999) / 100;
+ $R = "-R-$mwm/$mwm/-$mwm/$mwm";
+
+ GMT_begin(sprintf($corrPF,$id),"-JX$plotsize/$plotsize",$R,'-P -X6 -Y4');
+ GMT_psxy('-Ggrey80 -L');
+ printf(GMT "%g %g\n%g %g\n%g %g\n%g %g\n%g %g\n%g %g\n",
+ -$mwm, -$mwm+0.01,
+ -$mwm, -$mwm,
+ -$mwm+0.01, -$mwm,
+ $mwm, $mwm-0.01,
+ $mwm, $mwm,
+ $mwm-0.01, $mwm);
+
+ GMT_psxy('-W4,grey50');
+ print(GMT "-$mwm -$mwm\n$mwm $mwm\n");
+ GMT_psxy('-Sc0.12c -G255/127/80 -W1,blue,-');
+ for (my($bi)=0; $bi<@DL_dc_median; $bi++) {
+ next unless numberp($DL_dc_median[$bi]) && numberp($UL_dc_median[$bi]);
+ my($depth) = ($bi+0.5)*$opt_o;
+ last if ($depth > $opt_s);
+ print(GMT "$DL_dc_median[$bi] $UL_dc_median[$bi]\n");
+ }
+ GMT_psxy('-Sc0.12c -G255/127/80');
+ for (my($bi)=0; $bi<@DL_dc_median; $bi++) {
+ next unless numberp($DL_dc_median[$bi]) && numberp($UL_dc_median[$bi]);
+ my($depth) = ($bi+0.5)*$opt_o;
+ next unless ($depth > $opt_s);
+ print(GMT "$DL_dc_median[$bi] $UL_dc_median[$bi]\n");
+ }
+ GMT_psxy('-Sc0.12c -G46/139/87 -W1,blue,-');
+ for (my($bi)=0; $bi<@DL_uc_median; $bi++) {
+ next unless numberp($DL_uc_median[$bi]) && numberp($UL_uc_median[$bi]);
+ my($depth) = ($bi+0.5)*$opt_o;
+ last if ($depth > $opt_s);
+ print(GMT "$DL_uc_median[$bi] $UL_uc_median[$bi]\n");
+ }
+ GMT_psxy('-Sc0.12c -G46/139/87');
+ for (my($bi)=0; $bi<@DL_uc_median; $bi++) {
+ next unless numberp($DL_uc_median[$bi]) && numberp($UL_uc_median[$bi]);
+ my($depth) = ($bi+0.5)*$opt_o;
+ next unless ($depth > $opt_s);
+ print(GMT "$DL_uc_median[$bi] $UL_uc_median[$bi]\n");
+ }
+ GMT_pstext('-Gblue -N');
+ if (defined($outfile)) { printf(GMT "%f %f 14 0 0 TL $outfile [$P{run_label}]\n",-$mwm,1.1*$mwm); }
+ else { printf(GMT "%f %f 14 0 0 TL %03d\n [$P{run_label}]",$id,-$mwm,1.1*$wmw); }
+ GMT_psbasemap('-Bf0.01a0.05:"DL Vertical Velocity [m/s]":/f0.01a0.05:"UL Vertical Velocity [m/s]":WeSn');
+ GMT_end();
+
+ } # if dual_head
+
}
&antsExit(0);
--- a/LADCP_wspec Mon Jan 04 11:19:09 2016 +0000
+++ b/LADCP_wspec Tue Feb 02 14:55:24 2016 +0000
@@ -2,9 +2,9 @@
#======================================================================
# L A D C P _ W S P E C
# doc: Thu Jun 11 12:02:49 2015
-# dlm: Tue Oct 13 10:50:38 2015
+# dlm: Mon Jan 25 09:33:54 2016
# (c) 2012 A.M. Thurnherr
-# uE-Info: 21 43 NIL 0 0 72 10 2 4 NIL ofnI
+# uE-Info: 39 0 NIL 0 0 72 10 2 4 NIL ofnI
#======================================================================
$antsSummary = 'calculate VKE window spectra from LADCP profiles';
@@ -19,6 +19,7 @@
# - 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
($ANTSBIN) = (`which list` =~ m{^(.*)/[^/]*$});
($ANTSLIB) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
@@ -34,6 +35,7 @@
require "$ANTSLIB/nrutil.pl";
require "$ANTSBIN/.lsfit.poly";
require "$ANTSBIN/.nminterp.linear";
+&antsAddParams('LADCP_wspec',"Version $VERSION");
#----------------------------------------------------------------------
# Usage
--- a/acoustic_backscatter.pl Mon Jan 04 11:19:09 2016 +0000
+++ b/acoustic_backscatter.pl Tue Feb 02 14:55:24 2016 +0000
@@ -1,9 +1,9 @@
#======================================================================
# A C O U S T I C _ B A C K S C A T T E R . P L
# doc: Wed Oct 20 13:02:27 2010
-# dlm: Thu Jun 18 13:04:26 2015
+# dlm: Tue Jan 26 19:24:42 2016
# (c) 2010 A.M. Thurnherr
-# uE-Info: 27 82 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 171 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -25,6 +25,7 @@
# 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
#----------------------------------------------------------------------
# Volume Scattering Coefficient, following Deines (IEEE 1999)
@@ -96,9 +97,6 @@
my($range_to_bin) = abs($bd[$bin] - $LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH})
/ cos(rad($LADCP{ENSEMBLE}[$ens]->{TILT}))
/ $cosBeamAngle;
-# next
-# if ($range_to_bin < $SS_min_allowed_range ||
-# $range_to_bin > $SS_max_allowed_range);
my($temp) = defined($CTD_temp)
? $CTD{TEMP}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]
: $LADCP{ENSEMBLE}[$ens]->{TEMPERATURE};
@@ -144,6 +142,8 @@
my($bin) = $LADCP_firstBin-1;
my(@refSvProf,@refSvSamp,$depth,$i);
+ &antsAddParams('Sv_ref_bin',$Sv_ref_bin);
+
RETRY:
for ($depth=0; $depth<@nSv; $depth++) { # create reference profile
next unless ($nSv[$depth][$Sv_ref_bin-1] > 0);
@@ -168,11 +168,6 @@
}
info("\tusing bin %d as reference\n",$Sv_ref_bin);
-# for ($i=0; $i<@refSvProf; $i++) {
-# print(STDERR "$i $refSvProf[$i] $refSvSamp[$i]\n");
-# }
-# die;
-
my(@dSvProf); # create profiles for all bins
for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
my(@dSvSamp);
@@ -230,6 +225,9 @@
my($search_below) = int($_[0]); # grid index to begin search
my(@wdepth,@Sv_rng); # list of water_depth indices
+ &antsAddParams('SS_min_signal',$SS_min_signal,'SS_min_samp',$SS_min_samp,
+ 'SS_max_allowed_depth_range',$SS_max_allowed_depth_range);
+
for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) { # find backscatter min/max below $search_below in each bin
my($minSv,$maxSv,$depthmaxSv,$lastvalid) = (1e99,-1e99,-1,-1);
for (my($depth)=$search_below; $depth<@nSv; $depth++) {
--- a/bottom_tracking.pl Mon Jan 04 11:19:09 2016 +0000
+++ b/bottom_tracking.pl Tue Feb 02 14:55:24 2016 +0000
@@ -1,9 +1,9 @@
#======================================================================
# B O T T O M _ T R A C K I N G . P L
# doc: Wed Oct 20 21:05:37 2010
-# dlm: Tue Mar 4 13:54:15 2014
+# dlm: Tue Jan 26 15:25:14 2016
# (c) 2010 A.M. Thurnherr
-# uE-Info: 15 43 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 16 33 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -13,6 +13,7 @@
# Oct 24, 2011: - disabled not-very-useful %BT-params
# Apr 22, 2013: - replace output_bin_size by opt_o
# Mar 4, 2014: - removed old unused code
+# Jan 26, 2016: - added %PARAMs
# This code is essentially identical to the one used in LADCPproc. Differences:
# 1) velocity editing is simpler: no wake editing, no PPI editing, no shear
@@ -115,6 +116,10 @@
{
my($LADCP_start,$LADCP_end,$wd,$sig_wd) = @_;
+ &antsAddParams('BT_max_range',$BT_max_range,
+ 'BT_max_bin_range_diff',$BT_max_bin_range_diff,
+ 'BT_max_w_error',$BT_max_w_error);
+
for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
next unless ($wd-$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH} < $BT_max_range);
binBTprof($ens,$wd,$sig_wd);
--- a/defaults.pl Mon Jan 04 11:19:09 2016 +0000
+++ b/defaults.pl Tue Feb 02 14:55:24 2016 +0000
@@ -1,9 +1,9 @@
#======================================================================
# D E F A U L T S . P L
# doc: Tue Oct 11 17:11:21 2011
-# dlm: Mon Jan 4 10:53:34 2016
+# dlm: Thu Jan 28 16:55:23 2016
# (c) 2011 A.M. Thurnherr
-# uE-Info: 65 64 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 466 0 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -63,39 +63,41 @@
# Sep 26, 2015: - added sidelobe editing params
# Oct 13, 2015: - addded support for $ENV{VERB}
# Jan 4, 2016: - decreased default vertical resolution to 20m
+# Jan 22, 2016: - changed outGrid_firstBin default to 1
+# Jan 26, 2016: - removed -d
+# - changed outGrid_firstBin default to '*', also lastBin
+# Jan 27, 2016: - added documentation
#======================================================================
# Data Input
#======================================================================
-# CTD depth adjustment
-# - set with -d (-a up to 2013/05/16)
-# - value is added to CTD pressure
-# - use when CTD has -ve pressures
-
-&antsFloatOpt(\$opt_d,0);
-
-
-# set $opt_4 to 1 to suppress 3-beam LADCP solutions
+# Set $opt_4 to 1 (or use the -4 option) to suppress 3-beam LADCP
+# solutions
#$opt_4 = 1;
-# correct attiude sensors
+# The following variables allow bias-correcting the attiude
+# sensors.
# NB: heading is not used for vertical-velocity processing!
-$pitch_bias = $roll_bias = $heading_bias = 0;
+$pitch_bias = 0;
+$roll_bias = 0;
+$heading_bias = 0;
-# minimum valid velocities to require in a LADCP file
-# throws an error if not enough valid vels in a file
+# The following variable defines the minimum valid velocities
+# required in a LADCP file. If there are fewer data, an
+# error is produced
$min_valid_vels = 50;
-# bins to use in w calculations
-# - set with -b
-# - defaults to 2-last
+# The -b option defines the range of bins to use in w calculations.
+# The '*' indicates the last bin in the ADCP file. For data
+# collected with non-zero blanking distance, -b '1,*' should
+# likely be used.
$opt_b = '2,*' unless defined($opt_b);
@@ -114,50 +116,52 @@
$opt_v = 1 unless numberp($opt_v);
-# output base name
+# The "base name" of all output files (usually 0-padded 3-digits)
$out_basename = sprintf('%03d',$PROF);
-# output subdirectories
+# Output subdirectories
error("$RUN: no such directory\n") unless (-d $RUN);
$data_subdir = $plot_subdir = $log_subdir = $RUN;
-# min w samples required for each vertical-velocity bin
-# - value recorded in %outgrid_minsamp and used by [LADCP_w_regrid]
+# 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);
-# output grid resolution in meters
-# - value recorded in %outgrid_dz and used by [LADCP_w_regrid]
+# The -o option sets the output grid resolution in meters. The following
+# sets the default value.
&antsFloatOpt(\$opt_o,20);
-# the following variables limit the bins used to grid w_oean
+# 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 = 0;
-$outGrid_lastBin = 999;
+$outGrid_firstBin = '*'; # use $LADCP_firstBin (-b)
+$outGrid_lastBin = '*'; # use $LADCP_lastBin (-b)
#======================================================================
# Data Editing
#======================================================================
-# min correlation
+# The following sets the default correlation limit; measurements with
+# correlations below this limit are discarded.
&antsFloatOpt(\$opt_c,70);
-# max tilt (pitch/roll)
+# 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
@@ -177,56 +181,64 @@
&antsFloatOpt(\$opt_t,12);
-# max err vel
+# The following sets the default error velocity limit; measurements
+# with error velocities below this limit are discarded.
&antsFloatOpt(\$opt_e,0.1);
-# truncate farthest valid velocities
+# The following variable allows editing the velocities farthest
+# from the transducer. It defines how many velocities are to be
+# removed from each ensemble.
$truncate_farthest_valid_bins = 0;
-# discard velocities from chosen beam (1-4)
+# The following variable allows editing all data from a given
+# beam. Set to 1-4 to enable.
$bad_beam = 0;
-# max LADCP gap length in seconds
+# The following sets the maximum gap length in the w time series that
+# is simply ignored.
&antsFloatOpt(\$opt_g,60);
-# max allowed vertical velocity in m/s
+# The following variable sets the max allowed vertical ocean velocity
+# in m/s. Measurements with |w| this limit are discarded.
$w_max_lim = 1;
-# in each ensemble, vertical velocities differing more than this
-# parameter times mean absolute deviation from median, are considered
-# outliers and removed
+# In each ensemble, vertical velocities differing more than this
+# parameter times the mean absolute deviation from the median, are
+# considered outliers and removed.
$per_ens_outliers_mad_limit = 2;
-# data from bins with less valid velocities than the following parameter
-# are considered bad and removed
+# Data from bins with less valid velocities than the following parameter
+# are considered bad and removed. It is not clear whether this really
+# makes sense, but this editing is likely safe because it only affects
+# ensebles with the largest ranges.
$per_bin_valid_frac_lim = 0.15;
-# ensembles when instrument is shallower than
-# $surface_layer_depth in meters are removed.
-# possible contamination: ship's hull, thrusters, bubble clouds
-# Inspired by 2011_IWISE station 8
+# All ensembles recorded when the CTD is shallower than
+# the following parameter (depth in meters) are discarded.
+# Possible contamination: ship's hull, thrusters, bubble clouds
+# Inspired by 2011_IWISE station 8.
$surface_layer_depth = 25;
-# PPI editing as described in [edit_data.pl]
+# Previous Ping Interference editing as described in [edit_data.pl]
# - enabled by default for WH150 data
-# - variable sets an expression, to be evaluated once the data
-# have been loaded
+# - the variable defines a string with a perl expression, which is
+# evaluated once the data are loaded
# - 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
@@ -239,23 +251,21 @@
$PPI_editing_required = '($LADCP{BEAM_FREQUENCY} < 300)';
-# arbitrarily increase calculated max dist from seabed by 3%
-
-#$PPI_extend_upper_limit = 1.03;
+#$PPI_extend_upper_limit = 1.03; # see comments above
# The following variables control the "non-obvious" sidelobe editing for
# contamination from the seabed for the UL and from the sea surface for the
# DL. Tests with DoMORE-2 data (WH150 DL, WH300 UL) strongly suggest that
# it is not necessary to edit DL data for surface contamination. However,
-# at least for that instrument combination, UL contamination from the
+# at least for that instrument combination, UL (WH300) contamination from the
# seabed should clearly be removed.
$sidelobe_editing_DL_surface = 0;
$sidelobe_editing_UL_seabed = 1;
# The following variable sets the depth for sidelobe contamination
-# from the surface.
+# from the surface.
$vessel_draft = 6; # in meters
@@ -264,34 +274,42 @@
# Time Lagging
#======================================================================
-# externally supplied lag initial guess
+# The -i option allows defining an initial guess for the time lag between
+# the LADCP and the CTD data.
# $opt_i = 567;
-# reference layer bins for w for time matching
+# The following variables define the bins used to calculate the reference-
+# layer velocities used for time lagging.
($refLr_firstBin,$refLr_lastBin) = (2,6);
-# number of time lags during each of 2 lagging steps
+# The -n option defines the number of windows used to calculate
+# the optimal time lag. There's one value for each time-lagging step.
$opt_n = '10,100' unless defined($opt_n);
-# time lag search window widths for each of 2 lagging steps
-# - full width in seconds
+# The -w option defines the width of the window (in seconds) used
+# to calculate the optimal time lag. There's one value for each
+# time-lagging step.
$opt_w = '240,20' unless defined($opt_w);
-# if top 3 lags have spread greater than $TL_max_allowed_three_lag_spread
+# The following parameters control whether the top three time lags
+# are accepted or not.
+# If the top 3 lags have spread greater than $TL_max_allowed_three_lag_spread
# (in CTD scans) they must account for at least $TL_required_timelag_top_three_fraction
-# or there is an error
+# or an error is generated.
+# Notes:
# - $TL_max_allowed_three_lag_spread default was initially set to 2 but found to be
# violated quite often during 2011_IWISE
# - large spread may indicate dropped CTD scans
-# - the optimum value of $TL_max_allowed_three_lag_spread may be cast-duration dependent
+# - the optimum value of $TL_max_allowed_three_lag_spread may be
+# cast-duration dependent
$TL_max_allowed_three_lag_spread = 3;
&antsFloatOpt(\$opt_3,0.6);
@@ -319,7 +337,7 @@
# The following variable defines the minimum Sv signal in a bin (max - min)
-# required for reliable seabed detection. A limiting value of 40dB is
+# required for reliable seabed detection FROM ECHO AMPLITUDES. A limit of 40dB is
# indicated based on GoM#13, where the seabed is only visible in the last
# bin (#25). 30dB is chosen as the default to allow for variability.
# This value may need to be changed for data not collected with WH300
@@ -331,9 +349,9 @@
$SS_min_signal = 30;
-# Require at minimum nubmer of valid samples for seabed detection. Each
-# sample is a bin with a clear seabed maximum. With a proper setting
-# of $SS_min_signal, the algorithm is stable even with only a single
+# Require at minimum nubmer of valid samples for seabed detection FROM ECHO
+# AMPLITUDES. Each sample is a bin with a clear seabed maximum. With a proper
+# setting of $SS_min_signal, the algorithm is stable even with only a single
# sample (GoM#13). However, a default of 3 required samples is chosen
# to make seabed detection less sensitive to $SS_min_signal.
# This parameter is only used when $SS_use_BT == 0.
@@ -342,15 +360,16 @@
# The following numbers define the valid range of height-above bottom
-# for seabed detection. For data collected with WH300 instruments and
-# 8m bins, the maximum range needs to be greater than 250m (based on
-# GoM#13).
+# for seabed detection FROM ECHO AMPLITUDE. For data collected with WH300
+# instruments and 8m bins, the maximum range needs to be greater than 250m
+# (based on # GoM#13).
$SS_min_allowed_range = 0;
$SS_max_allowed_range = 350;
-# Number of ensembles around bottom to search. Only used with $SS_use_BT == 1.
+# Number of ensembles around bottom to search sabed IN BT DATA.
+# Only used with $SS_use_BT == 1.
$SS_search_window_halfwidth = 200;
@@ -436,11 +455,13 @@
# *_correlation.ps correlation time-depth plot
#----------------------------------------------------------------------
-@out_wsamp = ("plot_residuals($plot_subdir/${out_basename}_residuals.ps)",
- "plot_backscatter($plot_subdir/${out_basename}_backscatter.ps)",
-# "plot_correlation($plot_subdir/${out_basename}_correlation.ps)",
- "plot_wsamp($plot_subdir/${out_basename}_wsamp.ps)",
- "$data_subdir/$out_basename.wsamp");
+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)");
#----------------------------------------------------------------------
--- a/edit_data.pl Mon Jan 04 11:19:09 2016 +0000
+++ b/edit_data.pl Tue Feb 02 14:55:24 2016 +0000
@@ -1,9 +1,9 @@
#======================================================================
# E D I T _ D A T A . P L
# doc: Sat May 22 21:35:55 2010
-# dlm: Sat Sep 26 12:58:46 2015
+# dlm: Sun Jan 24 16:22:23 2016
# (c) 2010 A.M. Thurnherr
-# uE-Info: 34 56 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 433 31 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -32,6 +32,7 @@
# May 21, 2014: - got it to work correctly
# - croak -> error
# Sep 26, 2015: - added $vessel_draft to editSideLobes
+# Jan 23, 2016: - added &editBadTimeLagging()
# NOTES:
# - editCorr_Earthcoords() is overly conservative and removed most
@@ -416,4 +417,44 @@
return $nerm;
}
+
+#===============================================================================
+# $nerm = editBadTimeLagging($fromEns,$toEns,$good_from_elapsed1,$good_to_elapsed1,...)
+#===============================================================================
+
+sub editBadTimeLagging($$@)
+{
+ my($fe,$te,@elim) = @_;
+
+ my($nerm) = 0; # of ensembles removed
+ my($i) = 0;
+
+ if ($elim[0] < 0) { # entire piece is bad
+ for (my($e)=$fe; $e<=$te; $e++) {
+ undef($LADCP{ENSEMBLE}[$e]->{REFLR_W});
+ $nerm++;
+ }
+ } elsif (defined($elim[1])) { # limits in elim
+ my($e);
+ for ($e=$fe; @elim; shift(@elim),shift(@elim)) {
+# print(STDERR "deleting to $elim[0]\n");
+ while ($LADCP{ENSEMBLE}[$e]->{ELAPSED} < $elim[0]) {
+ undef($LADCP{ENSEMBLE}[$e]->{REFLR_W});
+ $nerm++;
+ $e++;
+ }
+# print(STDERR "keeping to $elim[1]\n");
+ while ($LADCP{ENSEMBLE}[$e]->{ELAPSED} < $elim[1]) { $e++; }
+ }
+# print(STDERR "deleting to $LADCP{ENSEMBLE}[$te]->{ELAPSED}\n");
+ while ($e <= $te) {
+ undef($LADCP{ENSEMBLE}[$e]->{REFLR_W});
+ $nerm++;
+ $e++;
+ }
+ }
+ return $nerm;
+}
+
+
1;
--- a/find_seabed.pl Mon Jan 04 11:19:09 2016 +0000
+++ b/find_seabed.pl Tue Feb 02 14:55:24 2016 +0000
@@ -1,9 +1,9 @@
#======================================================================
# F I N D _ S E A B E D . P L
# doc: Sun May 23 20:26:11 2010
-# dlm: Wed Oct 19 14:25:27 2011
+# dlm: Tue Jan 26 15:22:26 2016
# (c) 2010 A.M. Thurnherr
-# uE-Info: 15 0 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 16 33 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -13,6 +13,7 @@
# - increased z_offset from 10km to 15km
# Oct 19, 2011: - added $SS_max_allowed_range
# - renamed $SS_min_allowed_hab to *_range
+# Jan 26, 2016: - added %PARAMs
# NOTES:
# 1) BT range is corrected for sound speed at the transducer. This is not
@@ -36,6 +37,11 @@
my($i,$dd,$sd,$nd);
my(@guesses);
+ &antsAddParams('SS_min_allowed_range',$SS_min_allowed_range,
+ 'SS_max_allowed_range',$SS_max_allowed_range,
+ 'SS_search_window_halfwidth',$SS_search_window_halfwidth,
+ 'SS_max_allowed_depth_range',$SS_max_allowed_depth_range);
+
return undef unless ($be-$SS_search_window_halfwidth >= 0 &&
$be+$SS_search_window_halfwidth <= $#{$d->{ENSEMBLE}});
for ($i=$be-$SS_search_window_halfwidth; $i<=$be+$SS_search_window_halfwidth; $i++) {
--- a/plot_backscatter.pl Mon Jan 04 11:19:09 2016 +0000
+++ b/plot_backscatter.pl Tue Feb 02 14:55:24 2016 +0000
@@ -1,13 +1,14 @@
#======================================================================
# P L O T _ B A C K S C A T T E R . P L
# doc: Tue Jul 28 13:21:09 2015
-# dlm: Thu Jul 30 09:53:42 2015
+# dlm: Tue Jan 26 20:46:11 2016
# (c) 2015 A.M. Thurnherr
-# uE-Info: 57 59 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 19 31 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# HISTORY:
# Jul 28, 2015: - created from [LWplot_Sv]
+# Jan 26, 2016: - added return on no data to plot
require "$ANTS/libGMT.pl";
@@ -15,6 +16,8 @@
{
my($pfn) = @_;
+ return unless ($P{max_depth});
+
my($xmin) = $P{min_ens}-0.5;
my($xmax) = $P{max_ens}+0.5;
my($ymin) = 0;
--- a/plot_correlation.pl Mon Jan 04 11:19:09 2016 +0000
+++ b/plot_correlation.pl Tue Feb 02 14:55:24 2016 +0000
@@ -1,13 +1,14 @@
#======================================================================
# P L O T _ C O R R E L A T I O N . P L
# doc: Tue Jul 28 13:21:09 2015
-# dlm: Thu Jul 30 09:53:49 2015
+# dlm: Tue Jan 26 20:46:18 2016
# (c) 2015 A.M. Thurnherr
-# uE-Info: 57 59 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 19 31 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# HISTORY:
# Jul 28, 2015: - created from [LWplot_corr]
+# Jan 26, 2016: - added return on no data to plot
require "$ANTS/libGMT.pl";
@@ -15,6 +16,8 @@
{
my($pfn) = @_;
+ return unless ($P{max_depth});
+
my($xmin) = $P{min_ens}-0.5;
my($xmax) = $P{max_ens}+0.5;
my($ymin) = 0;
--- a/plot_mean_residuals.pl Mon Jan 04 11:19:09 2016 +0000
+++ b/plot_mean_residuals.pl Tue Feb 02 14:55:24 2016 +0000
@@ -1,9 +1,9 @@
#======================================================================
# P L O T _ M E A N _ R E S I D U A L S . P L
# doc: Tue Jul 28 13:21:09 2015
-# dlm: Thu Jul 30 12:38:12 2015
+# dlm: Tue Jan 26 20:13:56 2016
# (c) 2015 A.M. Thurnherr
-# uE-Info: 39 36 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 16 33 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -11,6 +11,9 @@
# Jul 29, 2015: - finished
# Jul 30, 2015: - added bin_tics
# - added outGrid_* support
+# Jan 22, 2015: - many changes
+# - added quality assessment label
+# Jan 25, 2016: - added return on no data
require "$ANTS/libGMT.pl";
@@ -18,6 +21,8 @@
{
my($pfn) = @_;
+ return unless ($P{BR_max_bin});
+
my($xmin) = -0.05;
my($xmax) = 0.05;
my($ymin) = 0.5;
@@ -26,7 +31,7 @@
my($R) = "-R$xmin/$xmax/$ymin/$ymax";
GMT_begin($pfn,'-JX10/-10',$R,'-P');
- if ($outGrid_firstBin>$LADCP_firstBin || $outGrid_lastBin<$LADCP_lastBin) {
+ if ($outGrid_firstBin>$LADCP_firstBin || $outGrid_lastBin<$LADCP_lastBin) { # mark used bins
GMT_psxy('-G200 -M -L');
printf(GMT ">\n%f %f\n%f %f\n%f %f\n%f %f\n",
$xmin,$LADCP_firstBin-0.5,
@@ -42,37 +47,56 @@
if ($outGrid_lastBin<$LADCP_lastBin);
}
- GMT_psxy('-W1');
+ GMT_psxy('-W1'); # plot zero line
printf(GMT "0 $ymin\n0 $ymax\n");
GMT_psxy('-Mn -W4/coral');
- for (my($bin)=0; $bin<scalar(@dc_bres); $bin++) {
+ for (my($bin)=$LADCP_firstBin; $bin<@dc_bres; $bin++) { # SKIP FIRST BIN
+ next if ($bin+1<$outGrid_firstBin || $bin+1>$outGrid_lastBin);
+ next unless ($dc_bres_nsamp[$bin] >= $dc_bres_max_nsamp/3);
printf(GMT "%f %d\n",$dc_avg_bres[$bin],$bin+1);
}
GMT_psxy('-Mn -Ex0.2/4/coral');
- for (my($bin)=0; $bin<scalar(@dc_bres); $bin++) {
+ for (my($bin)=$LADCP_firstBin-1; $bin<@dc_bres; $bin++) {
printf(GMT "%f %d %f\n",
$dc_avg_bres[$bin],
$bin+1,
- (scalar(@{$dc_bres[$bin]}) > 1) ?
- $dc_sig_bres[$bin]/sqrt(scalar(@{$dc_bres[$bin]})-1) : 0);
+ ($dc_bres_nsamp[$bin] > 1) ?
+ $dc_sig_bres[$bin]/sqrt($dc_bres_nsamp[$bin]-1) : 0);
}
GMT_psxy('-Mn -W4/SeaGreen');
- for (my($bin)=0; $bin<scalar(@uc_bres); $bin++) {
+ for (my($bin)=$LADCP_firstBin; $bin<@uc_bres; $bin++) { # SKIP FIRST BIN
+ next if ($bin+1<$outGrid_firstBin || $bin+1>$outGrid_lastBin);
+ next unless ($uc_bres_nsamp[$bin] >= $uc_bres_max_nsamp/3);
printf(GMT "%f %d\n",$uc_avg_bres[$bin],$bin+1);
}
GMT_psxy('-Mn -Ex0.2/4/SeaGreen');
- for (my($bin)=0; $bin<scalar(@uc_bres); $bin++) {
+ for (my($bin)=$LADCP_firstBin-1; $bin<@uc_bres; $bin++) {
printf(GMT "%f %d %f\n",
$uc_avg_bres[$bin],
$bin+1,
- (scalar(@{$uc_bres[$bin]}) > 1) ?
- $uc_sig_bres[$bin]/sqrt(scalar(@{$uc_bres[$bin]})-1) : 0);
+ ($uc_bres_nsamp[$bin] > 1) ?
+ $uc_sig_bres[$bin]/sqrt($uc_bres_nsamp[$bin]-1) : 0);
}
GMT_unitcoords(); # LABELS
- GMT_pstext(-Gblue);
- print(GMT "0.02 0.98 12 0 0 BL $P{out_basename} $P{run_label}\n");
+
+ GMT_pstext('-Gblue -N'); # profile id
+ print(GMT "0.0 -0.03 14 0 0 BL $P{out_basename} $P{run_label}\n");
+
+ GMT_pstext('-Gcoral'); # rms residuals
+ print(GMT "0.01 0.93 12 0 0 BL dc\n");
+ if ($dc_bres_rms >= 0.005) { GMT_pstext('-Gwhite -Wred'); }
+ elsif ($dc_bres_rms >= 0.001) { GMT_pstext('-Gblack -Wyellow'); }
+ else { GMT_pstext('-Gblack -Wgreen'); }
+ printf(GMT "0.10 0.93 12 0 0 BL %.1f mm/s rms\n",1000*$dc_bres_rms);
+
+ GMT_pstext('-GSeaGreen');
+ print(GMT "0.01 0.98 12 0 0 BL uc\n");
+ if ($uc_bres_rms >= 0.005) { GMT_pstext('-Gwhite -Wred'); }
+ elsif ($uc_bres_rms >= 0.001) { GMT_pstext('-Gblack -Wyellow'); }
+ else { GMT_pstext('-Gblack -Wgreen'); }
+ printf(GMT "0.10 0.98 12 0 0 BL %.1f mm/s rms\n",1000*$uc_bres_rms);
my($bin_tics) = ($ymax <= 20) ? 'f1a1' : 'f1a2';
GMT_setR($R); # FINISH PLOT
--- a/plot_residuals.pl Mon Jan 04 11:19:09 2016 +0000
+++ b/plot_residuals.pl Tue Feb 02 14:55:24 2016 +0000
@@ -1,15 +1,16 @@
#======================================================================
# P L O T _ R E S I D U A L S . P L
# doc: Tue Jul 28 13:21:09 2015
-# dlm: Thu Jul 30 09:52:39 2015
+# dlm: Tue Jan 26 20:45:58 2016
# (c) 2015 A.M. Thurnherr
-# uE-Info: 73 59 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 21 31 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# HISTORY:
# Jul 28, 2015: - created from [LWplot_residuals]
# Jul 30, 2015: - made it respect outGrid_ selection
# - modified $ens_tics
+# Jan 26, 2016: - added return on no data to plot
require "$ANTS/libGMT.pl";
@@ -17,6 +18,8 @@
{
my($pfn) = @_;
+ return unless ($P{max_depth});
+
my($xmin) = $P{min_ens}-0.5;
my($xmax) = $P{max_ens}+0.5;
my($ymin) = 0;
--- a/plot_time_lags.pl Mon Jan 04 11:19:09 2016 +0000
+++ b/plot_time_lags.pl Tue Feb 02 14:55:24 2016 +0000
@@ -1,13 +1,14 @@
#======================================================================
# P L O T _ T I M E _ L A G S . P L
# doc: Tue Jul 28 13:21:09 2015
-# dlm: Wed Jul 29 14:47:57 2015
+# dlm: Tue Jan 26 20:14:53 2016
# (c) 2015 A.M. Thurnherr
-# uE-Info: 39 30 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 19 38 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
# Jul 29, 2015: - created from [LWplot_TL]
+# Jan 26, 2016: - added return on no data to plot
require "$ANTS/libGMT.pl";
@@ -15,6 +16,8 @@
{
my($pfn) = @_;
+ return unless ($P{'elapsed.min'});
+
my($xmin) = $P{'elapsed.min'}/60;
my($xmax) = $P{'elapsed.max'}/60;
my($ymin) = -24;
@@ -34,13 +37,11 @@
printf(GMT "%f %f\n",$elapsed_buf[$wi]/60,$so_buf[$wi]);
}
- my($fel) = $P{min_elapsed}; # from-elapsed limit
GMT_psxy('-W4/grey20 -M');
for (my($i)=0; $i<@bmo_buf; $i++) {
printf(GMT ">\n%f %f\n%f %f\n",
- $fel/60, $bmo_buf[$i],
- $te_buf[$i]/60+1,$bmo_buf[$i]);
- $fel = $te_buf[$i];
+ $fg_buf[$i]/60-0.5,$bmo_buf[$i],
+ $lg_buf[$i]/60+0.5,$bmo_buf[$i]);
}
GMT_unitcoords(); # LABELS
--- a/plot_wsamp.pl Mon Jan 04 11:19:09 2016 +0000
+++ b/plot_wsamp.pl Tue Feb 02 14:55:24 2016 +0000
@@ -1,9 +1,9 @@
#======================================================================
# P L O T _ W S A M P . P L
# doc: Tue Jul 28 13:21:09 2015
-# dlm: Mon Oct 12 13:37:27 2015
+# dlm: Tue Jan 26 20:46:47 2016
# (c) 2015 A.M. Thurnherr
-# uE-Info: 13 53 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 22 31 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -11,6 +11,7 @@
# Jul 30, 2015: - added support for outGrid_*
# Sep 21, 2015: - BUG: function was still called plot_w()
# Oct 12, 2015: - move main label outside plot area
+# Jan 26, 2016: - added return on no data to plot
require "$ANTS/libGMT.pl";
@@ -18,6 +19,8 @@
{
my($pfn) = @_;
+ return unless ($P{max_depth});
+
my($xmin) = $P{min_ens}-0.5;
my($xmax) = $P{max_ens}+0.5;
my($ymin) = 0;
--- a/time_lag.pl Mon Jan 04 11:19:09 2016 +0000
+++ b/time_lag.pl Tue Feb 02 14:55:24 2016 +0000
@@ -1,9 +1,9 @@
#======================================================================
# T I M E _ L A G . P L
# doc: Fri Dec 17 21:59:07 2010
-# dlm: Wed Jul 29 14:43:08 2015
+# dlm: Tue Jan 26 15:04:54 2016
# (c) 2010 A.M. Thurnherr
-# uE-Info: 272 59 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 68 33 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -60,6 +60,12 @@
# May 15, 2015: - fiddled with assertions
# Jun 19, 2015: - disabled L2 warning on partial-depth time-lagging failures
# Jul 29, 2015: - support for new plotting system
+# Jan 22, 2016: - started adding support for timelag-filtering
+# Jan 23, 2016: - continued
+# Jan 24, 2016: - made it work
+# - BUG: time-lag plot was not produced when final lag piece had problems
+# Jan 25, 2016: - search-radius-doubling heuristic had typo
+# - added %PARAMs
# DIFFICULT STATIONS:
# NBP0901#131 this requires the search-radius doubling heuristic
@@ -127,17 +133,21 @@
{ # STATIC SCOPE
-local(@elapsed_buf,@so_buf,@mad_buf,@bmo_buf,@te_buf,$elapsed_min_buf); # available to plot routines
+local(@elapsed_buf,@so_buf,@mad_buf,@bmo_buf,@fg_buf,$lg_buf,$elapsed_min_buf); # available to plot routines
sub calc_lag($$$$$)
{
my($n_windows,$w_size,$scan_increment,$first_ens,$last_ens) = @_;
my($search_radius) = $scan_increment==1 ? 3 : $w_size;
+ &antsAddParams('TL_max_allowed_three_lag_spread',$TL_max_allowed_three_lag_spread);
+
my($ctmsg);
if ($first_ens==$firstGoodEns && $last_ens==$lastGoodEns) { $ctmsg = "full-cast"; }
else { $ctmsg = "partial-cast"; }
+ my($last_lag_piece) = ($last_ens == $lastGoodEns); # none is following
+
RETRY:
my($failed) = undef;
progress("Calculating $n_windows $ctmsg time lags from ${w_size}s-long windows at %dHz resolution...\n",
@@ -160,11 +170,10 @@
$first_ens = $approx_joint_profile_start_ens
if ($first_ens < $approx_joint_profile_start_ens);
- my($last_lag_piece) = ($last_ens == $lastGoodEns); # none is following
$last_ens = $approx_joint_profile_end_ens
if ($last_ens > $approx_joint_profile_end_ens);
- for (my($wi)=0; $wi<$n_windows; $wi++) {
+ for (my($wi)=0; $wi<$n_windows; $wi++) { # use bestLag() in each window
my($fe) = $first_ens + int(($last_ens-$first_ens-$window_ens)*$wi/($n_windows-1)+0.5);
die("assertion failed\n\tfe = $fe, first_ens = $first_ens, last_ens = $last_ens, window_ens = $window_ens, firstGoodEns = $firstGoodEns, lastGoodEns = $lastGoodEns")
unless ($fe>=$firstGoodEns && $fe+$window_ens<=$lastGoodEns);
@@ -188,7 +197,7 @@
$nBest{$lag} = 0;
}
- my(@best_lag);
+ my(@best_lag); # find 3 most popular lags
foreach my $lag (keys(%nBest)) {
$best_lag[0] = $lag if ($nBest{$lag} > $nBest{$best_lag[0]});
}
@@ -200,21 +209,23 @@
next if ($lag == $best_lag[0] || $lag == $best_lag[1]);
$best_lag[2] = $lag if ($nBest{$lag} > $nBest{$best_lag[2]});
}
- if ($nBest{$best_lag[2]}) {
+ if ($nBest{$best_lag[2]}) { # there are at least 3 lags
progress("\t3 most popular offsets: %d (%d%% %.1fcm/s mad), %d (%d%% %.1fcm/s mad), %d (%d%% %.1fcm/s mad)\n",
$best_lag[0],int(($nBest{$best_lag[0]}/$n_valid_windows)*100+0.5),100*$madBest{$best_lag[0]},
$best_lag[1],int(($nBest{$best_lag[1]}/$n_valid_windows)*100+0.5),100*$madBest{$best_lag[1]},
$best_lag[2],int(($nBest{$best_lag[2]}/$n_valid_windows)*100+0.5),100*$madBest{$best_lag[2]});
- } elsif ($nBest{$best_lag[1]}) {
+ } elsif ($nBest{$best_lag[1]}) { # there are only 2 lags
progress("\toffsets: %d (%d%% %.1fcm/s mad), %d (%d%% %.1fcm/s mad)\n",
$best_lag[0],int(($nBest{$best_lag[0]}/$n_valid_windows)*100+0.5),100*$madBest{$best_lag[0]},
$best_lag[1],int(($nBest{$best_lag[1]}/$n_valid_windows)*100+0.5),100*$madBest{$best_lag[1]});
- } else {
+ } else { # there is only 1 lag
progress("\toffset: %d (%d%% %.1fcm/s mad)\n",
$best_lag[0],int(($nBest{$best_lag[0]}/$n_valid_windows)*100+0.5),100*$madBest{$best_lag[0]});
}
- unless ($nBest{$best_lag[0]}+$nBest{$best_lag[1]}+$nBest{$best_lag[2]} >= $opt_3*$n_valid_windows) {
+
+ unless ($nBest{$best_lag[0]}+$nBest{$best_lag[1]}+$nBest{$best_lag[2]} # require quorum
+ >= $opt_3*$n_valid_windows) {
if (max(@best_lag)-min(@best_lag) > $TL_max_allowed_three_lag_spread) {
warning(2,"$0: cannot determine a valid $ctmsg lag; top 3 tags account for %d%% of total (use -3 to relax criterion)\n",
int(100*($nBest{$best_lag[0]}+$nBest{$best_lag[1]}+$nBest{$best_lag[2]})/$n_valid_windows+0.5))
@@ -226,16 +237,15 @@
}
}
- my($bmo) = $best_lag[0];
-
- if ($bmo > 0.9*$search_radius/2/$CTD{DT}) {
+ my($bmo) = $best_lag[0]; # best mean offset
+ if ($bmo > 0.9*$search_radius/2/$CTD{DT}) { # cannot be near edge of window
if ($search_radius == $w_size) {
warning(0,"lag too close to edge of search --- trying again after shifting the initial offset\n");
$CTD{TIME_LAG} += $search_radius/2;
} else {
warning(0,"lag too close to edge of search --- trying again after doubling the search radius\n");
$search_radius *= 2;
- $search_radius =- $w_size if ($search_radius > $w_size);
+ $search_radius = $w_size if ($search_radius > $w_size);
}
undef(%nBest); undef(%madBest); undef(@best_lag);
goto RETRY;
@@ -247,28 +257,123 @@
} else {
warning(0,"lag too close to edge of search --- trying again after doubling the search radius\n");
$search_radius *= 2;
- $search_radius =- $w_size if ($search_radius > $w_size);
+ $search_radius = $w_size if ($search_radius > $w_size);
}
undef(%nBest); undef(%madBest); undef(@best_lag);
goto RETRY;
+ }
+
+ #----------------------------------------------------
+ # Here, either $failed is set, or we have a valid lag.
+ # If we have a valid lag, a continuous range of good
+ # lags is determined using a finite-state machine.
+ # state == 0 no good run found yet
+ # state == 1 good run found, $fg is set
+ # A good run is at least $min_runlength long,
+ # and every $scan_runlength-long sequence contains at
+ # least $min_good scan offsets that
+ # agree with the median offset within +/- $good_diff.
+ #----------------------------------------------------
+
+ my(@fg,@lg);
+ my($min_runlength) = 7; my($scan_runlength) = 7; my($min_good) = 4; my($good_diff) = 2;
+ unless ($failed || $scan_increment>1) {
+ my($state) = 0;
+ for (my($i)=0; 1; $i++) {
+# printf(STDERR "$i: state = $state\n");
+ if ($state == 0) {
+ last if ($i >= @elapsed-$scan_runlength);
+ my($ngood) = 0;
+ for (my($j)=0; $j<$scan_runlength; $j++) {
+ $ngood += (abs($bmo-$so[$i+$j]) <= $good_diff);
+ }
+# printf(STDERR "$i: ngood = $ngood\n");
+ if ($ngood >= $min_good) { # we want at least 3 out of 5
+ $state = 1;
+ if ($i == 0) { # run begins at start
+ push(@fg,0);
+ } else { # run begins at first matching offset
+ my($fg);
+ for (my($j)=0; $j<$scan_runlength; $j++) {
+ $fg = $i+$scan_runlength-1-$j
+ if (abs($bmo-$so[$i+$scan_runlength-1-$j]) <= $good_diff);
+ }
+ push(@fg,$fg);
+ }
+ }
+ } elsif ($state == 1) { # growing run
+ die("assertion failed (i = $i)")
+ if ($i > @elapsed-$scan_runlength);
+ if ($i == @elapsed-$scan_runlength) { # run extends to end
+ push(@lg,$#elapsed);
+ last;
+ }
+ my($ngood) = 0;
+ for (my($j)=0; $j<$scan_runlength; $j++) {
+ $ngood += (abs($bmo-$so[$i+$j]) <= $good_diff);
+ }
+# printf(STDERR "$i: ngood = $ngood\n");
+ if ($ngood < $min_good) { # run ended
+ my($lg);
+ for (my($j)=0; $j<$scan_runlength; $j++) {
+ $lg = $i if (abs($bmo-$so[$i+$j]) <= 1);
+ }
+ push(@lg,$lg);
+ $state = 0;
+ }
+ } # if state == 1
+ } # for i
+# printf(STDERR "%d runs found\n",scalar(@lg));
+ } # unless $failed || scan_increment > 1
+
+ #--------------------------------------------------
+ # Filter LADCP data for measurements during times
+ # of uncertain time lags
+ #--------------------------------------------------
+
+ if ($scan_increment == 1) {
+ progress("\tEditing data with unknown time-lags...\n");
+ my(@elim);
+# print(STDERR "fg = @fg; lg = @lg\n");
+ for (my($i)=0; $i<@fg; $i++) {
+ next if ($lg[$i]-$fg[$i] < $min_runlength);
+ push(@elim,($fg[$i] == 0) ? $elapsed[$first_ens] : $elapsed[$fg[$i]],
+ ($lg[$i] == $n_windows-1) ? $elapsed[$last_ens] : $elapsed[$lg[$i]]);
+ }
+# print(STDERR "elim = @elim\n");
+ $nerm = $failed
+ ? editBadTimeLagging($first_ens,$last_ens,-1)
+ : editBadTimeLagging($first_ens,$last_ens,@elim);
+ progress("\t\t$nerm ensembles removed (%d%% of total), leaving %d run\n",round(100*$nerm/($last_ens-$first_ens+1)),scalar(@elim)/2);
}
+ #------------------------------------------------------
+ # Produce plot on fine-grained time-lagging
+ # - accumulate data into plot buffer
+ # - on last lag piece (usually upcast), plot is drawn
+ #------------------------------------------------------
+
if (@out_TL && $scan_increment==1) {
push(@elapsed_buf,@elapsed); # buffer elapsed data in static scope
push(@so_buf,@so); # scan offset
push(@mad_buf,@mad); # mean absolute deviation
- push(@bmo_buf,$bmo); # best median offset
- push(@te_buf,$elapsed[$#elapsed]); # to elapsed (from elapsed to elapsed, capisc?)
+
+ for (my($i)=0; $i<@fg; $i++) {
+ next if ($lg[$i]-$fg[$i] < $min_runlength);
+ push(@bmo_buf,$bmo); # best median offset (copy for each run)
+ push(@fg_buf,$elapsed[$fg[$i]]); # first good so in lag-piece
+ push(@lg_buf,$elapsed[$lg[$i]]); # last good so in lag-piece
+ }
$elapsed_min_buf = $elapsed[0] # min of valid elapsed (for plotting)
unless defined($elapsed_min_buf);
- if ($last_lag_piece) { # output all data
- progress("\tWriting time-lagging time series to ");
+ if ($last_lag_piece) {
+ progress("\tWriting time-lagging time series to "); # output all data
my($saveParams) = $antsCurParams;
@antsNewLayout = ('elapsed.LADCP','scan_offset','mad','downcast');
&antsAddParams('best_scan_offsets',"@bmo_buf");
- &antsAddParams('to_elapsed_limits',"@te_buf");
+# &antsAddParams('to_elapsed_limits',"@te_buf");
&antsAddParams('elapsed.min',$elapsed_min_buf);
&antsAddParams('elapsed.max',$elapsed_buf[$#elapsed_buf]);
&antsAddParams('elapsed.bot',$LADCP{ENSEMBLE}[$LADCP_atbottom]->{ELAPSED});
--- a/version.pl Mon Jan 04 11:19:09 2016 +0000
+++ b/version.pl Tue Feb 02 14:55:24 2016 +0000
@@ -1,13 +1,19 @@
#======================================================================
# V E R S I O N . P L
# doc: Tue Oct 13 10:40:57 2015
-# dlm: Mon Jan 4 10:57:14 2016
+# dlm: Mon Jan 25 09:28:36 2016
# (c) 2015 A.M. Thurnherr
-# uE-Info: 9 0 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 12 0 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
-#$VERSION = '1.1beta1'; # Oct 12, 2015
-$VERSION = '1.1'; # Jan 4, 2016
+# HISTORY:
+# Oct 13, 2015: - created
+# Jan 4, 2016: - added $ADCP_tools_minVersion
-$antsMinLibVersion = 6.2;
+#$VERSION = '1.1beta1'; # Oct 12, 2015
+#$VERSION = '1.1'; # Jan 4, 2016
+$VERSION = '1.2beta'; # Jan 5, 2016
+$antsMinLibVersion = 6.2;
+$ADCP_tools_minVersion = 1.4;
+