--- a/.hgignore Thu Mar 16 11:53:27 2017 -0400
+++ b/.hgignore Tue Nov 27 16:59:05 2018 -0500
@@ -1,9 +1,9 @@
#======================================================================
# . H G I G N O R E
# doc: Thu Oct 27 10:27:55 2011
-# dlm: Tue Dec 22 11:01:28 2015
+# dlm: Wed Oct 31 10:22:46 2018
# (c) 2011 A.M. Thurnherr
-# uE-Info: 16 14 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 20 16 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
syntax: glob
@@ -14,6 +14,10 @@
Plots/
Documentation/
+ADCP_tools/
+ANTSlib/
+Makefile.Archive
+LADCP_w_Software.tgz
Dec_17_2010/
Dec_24_2010/
--- a/HISTORY Thu Mar 16 11:53:27 2017 -0400
+++ b/HISTORY Tue Nov 27 16:59:05 2018 -0500
@@ -1,9 +1,9 @@
======================================================================
H I S T O R Y
doc: Mon Oct 12 16:09:24 2015
- dlm: Thu Mar 16 11:53:07 2017
+ dlm: Tue Nov 27 16:57:47 2018
(c) 2015 A.M. Thurnherr
- uE-Info: 260 19 NIL 0 0 72 3 2 4 NIL ofnI
+ uE-Info: 372 37 NIL 0 0 72 3 2 4 NIL ofnI
======================================================================
----------------------------------------------------------------------
@@ -258,3 +258,116 @@
Mar 16, 2017:
- published
+
+----------------------------------------------------------------------
+
+ Oct 3, 2017:
+ - added -q option to LADCP_w_CTD for reprocessing as 1Hz files
+
+ Oct 12, 2017:
+ - re-wrote code in [LADCP_w_ocean] to deal with earth coordinates
+ (MAJOR BUG FIX)
+
+ Oct 13, 2017:
+ - bugfix in [edit_data.pl]
+
+ Oct 17, 2017:
+ - improvements and bugfix in [LADCP_VKE]
+
+ Nov 26, 2017:
+ - significant bug fixes in [LADCP_w_ocean] related to
+ - ping-coherent residual removal
+ - bad beam
+
+ Nov 27, 2017:
+ - improved gap heuristics in [time_series.pl]
+ - added @valid_ensmeble_range [defaults.pl]
+
+ Nov 28, 2017:
+ - increased version to V1.4 [version.pl] [.hg/hgrc]
+ - worked on updating howto
+ - improvements to [LADCP_w_ocean]
+ - removed wcorr plot from [LADCP_w_postproc]
+
+ Nov 29, 2017:
+ - replaced opt_i by initial_time_lag in [defaults.pl]
+
+ Dec 9, 2017:
+ - removed common options from [LADCP_w_ocean] [LADCP_w_postproc]
+ [LADCP_wspec] [LADCP_w_CTD] [LADCP_VKE]
+
+ Dec 14, 2017:
+ - improvements to [LADCP_wspec]
+
+ Dec 17, 2017:
+ - added dependencies to [LADCP_w_ocean] [LADCP_w_CTD]
+
+ Mar 8, 2018:
+ - improvements to [LADCP_w_CTD]
+
+ Mar 9, 2018:
+ - improvements to [LADCP_w_CTD]
+
+ Mar 20, 2018:
+ - added blue background color for in-ice profiles in [plot_wprof.pl]
+ - fixed acceleration unit error in [plot_wprof.pl]
+
+ Mar 22, 2018:
+ - adapted howto
+ - massively improved time-lagging heuristic [time_lag.pl]
+ - improved [plot_time_lags.pl]
+
+ Mar 27, 2018:
+ - bugfix to new time lagging heuristic in [time_lag.pl]
+
+ Apr 24, 2018:
+ - improvements to [LADCP_w_ocean] [LADCP_VKE] [defaults.pl]
+
+ Apr 25, 2018:
+ - improvement to [LADCP_VKE]
+
+ May 1, 2018:
+ - added reflr_u filter in [LADCP_w_ocean] [edit_data.pl] [time_series.pl]
+ - added ambiguity velocity check in [LADCP_w_ocean]
+
+ May 2, 2018:
+ - bug fixes in [LADCP_w_ocean], related to
+ - reflr threshold
+ - PPI
+ - adapted [defaults.pl] to reflr_u filter
+
+ May 13, 2018:
+ - fixed bug in [LADCP_wspec]
+
+ May 16, 2018:
+ - improvement to [LADCP_wspec]
+
+ May 22, 2018:
+ - improvement to [LADCP_w_CTD]
+
+ Jul 24, 2018:
+ - improvement to [LADCP_w_CTD]
+
+ Sep 13, 2018:
+ - added '.' to library path in [version.pl]
+
+ Oct 4, 2018:
+ - improvements to [LADCP_w_CTD]
+
+ Oct 5, 2018:
+ - improvements and bugfix in [LADCP_w_CTD]
+
+ Oct 31, 2018:
+ - improvements to [LADCP_w_postproc]
+
+ Nov 1, 2018:
+ - improvements to [LADCP_w_postproc]
+
+ Nov 2, 2018:
+ - 2-beam residuals bug fix in [LADCP_w_ocean]
+
+ Nov 17, 2018:
+ - updated [HISTORY]
+ - updated prerequisites in [version.pl]
+ - updated [LADCP_w_howto.pdf]
+
--- a/LADCP_VKE Thu Mar 16 11:53:27 2017 -0400
+++ b/LADCP_VKE Tue Nov 27 16:59:05 2018 -0500
@@ -2,9 +2,9 @@
#======================================================================
# L A D C P _ V K E
# doc: Tue Oct 14 11:05:16 2014
-# dlm: Tue Mar 14 17:38:00 2017
+# dlm: Tue Jul 24 17:02:30 2018
# (c) 2012 A.M. Thurnherr
-# uE-Info: 155 45 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 399 0 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
$antsSummary = 'calculate VKE from LADCP-derived vertical-velocity profiles';
@@ -103,6 +103,11 @@
# Mar 13, 2017: - added -a)mbient <eps>
# Mar 14, 2017: - disabled -a) by default, because -a 0 is clearly bad and
# I have no evidence yet that -a something is better than -l 0
+# Oct 17, 2017: - added default 'eps' field on -k
+# - added eps.ms field to files processed without -k
+# Dec 9, 2017: - added support for $antsSuppressCommonOptions
+# Apr 24, 2018: - BUG: output was one field too wide (filled with nans) because antsBufNFields was not reset
+# Apr 25, 2018: - added -y and removed spectral bins from default output
($ANTSLIB) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
($WCALC) = ($0 =~ m{^(.*)/[^/]*$});
@@ -137,7 +142,8 @@
# Usage
#----------------------------------------------------------------------
-&antsUsage('a:bc:de:f:g:i:k:l:mno:p:q:r:s:tuw:x:z:',0,
+$antsSuppressCommonOptions = 1;
+&antsUsage('a:bc:de:f:g:i:k:l:mno:p:q:r:s:tuw:x:yz:',0,
"[poly-o)rder <n[$opt_o]> to de-mean data; -1 to disable>] [apply cosine-t)aper]",
'[-d)own/-u)pcast-only] [exclude -b)ottom window]', # LADCP_wspec options
"[-s)urface <layer depth to exclude[${opt_s}m]>",
@@ -150,6 +156,7 @@
'[o-m)it spectral correction] [spectral-tilt-correction -r)ange <max[0m]>]',
"[-e)ps-parameterization <constant[${c}s^-0.5]>",
'[include microstructure -k)e dissipation <file:field> in _VKE plot]',
+ '[-y) record spectra in output file]',
'[write output -f)iles to <directory>]',
'[write output filed with -i)ndividual spectra <basename>]',
'[output -p)lot <ps-file[#_VKE.ps]>]',
@@ -239,6 +246,7 @@
undef(@antsLayout); # shouldn't matter, because it will get overwritten
undef($antsOldHeaders); # forget those
undef(@ants_);
+ $antsBufNFields = 0;
} elsif (defined(fnrNoErr('pwrdens.0'))) {
croak("$0: -d, -u, -b, -w, -s meaningless when $0 used with spectral input\n")
@@ -398,7 +406,6 @@
my($DOF) = 0;
-
my($sumd,$sumx,$sumy) = (0,0,0); # fit kz^-2 power law
for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) {
my($i) = $f - $pg_fmin;
@@ -455,6 +462,7 @@
my(@eps_ms,@depth_ms); # output variables
if (defined($opt_k)) {
my($file,$field) = split(':',$opt_k);
+ $field = 'eps' unless defined($field);
open(ADDF,"$file") || croak("$file: $!\n"); # open file
my(@afl) = &antsFileLayout(ADDF); # read layout
my($akf,$aef);
@@ -485,8 +493,7 @@
$slpf = &antsNewField('p0fit.slope'); # power-law slope
$sslpf = &antsNewField('p0fit.slope.sig'); # power-law slope stddev
$wepsf = &antsNewField('eps.VKE'); # epsilon from VKE
-$msepsf = &antsNewField('eps.ms') # externally supplied microstructure eps
- if defined($opt_k);
+$msepsf = &antsNewField('eps.ms'); # externally supplied microstructure eps if available
my(@outLayout) = @antsNewLayout; # save for later
for ($f=0; $f<@outLayout; $f++) { # determine last spectral field in input
@@ -599,7 +606,9 @@
$sum += $eps_ms[$i]; $n++;
}
$ants_[$r][$msepsf] = $n ? $sum / $n : nan;
- }
+ } else {
+ $ants_[$r][$msepsf] = nan;
+ }
#---------------
# produce output
@@ -727,5 +736,12 @@
$antsOldHeaders = $Hbuf;
$antsHeadersPrinted = 0;
-&antsFlush(); # output record with results
+unless (defined($opt_y)) { # remove spectral bins from output
+ splice(@antsNewLayout,$pg_fmin,$nfreq);
+ for (my($r)=0; $r<@ants_; $r++) {
+ splice(@{$ants_[$r]},$pg_fmin,$nfreq);
+ }
+}
+
+&antsFlush(); # output results
&antsExit();
--- a/LADCP_w_CTD Thu Mar 16 11:53:27 2017 -0400
+++ b/LADCP_w_CTD Tue Nov 27 16:59:05 2018 -0500
@@ -2,9 +2,9 @@
#======================================================================
# L A D C P _ W _ C T D
# doc: Mon Nov 3 17:34:19 2014
-# dlm: Fri Jul 29 11:20:47 2016
+# dlm: Fri Oct 5 14:52:00 2018
# (c) 2014 A.M. Thurnherr
-# uE-Info: 69 48 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 85 93 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
$antsSummary = 'pre-process SBE 9plus CTD data for LADCP_w';
@@ -63,6 +63,27 @@
# May 26, 2016: - renamed w[.raw] to CTD_w[.raw]
# - added winch_w
# Jul 29, 2016: - change w_CTD plot label to reflect sign convention
+# Oct 3, 2017: - cosmetics
+# - added -q to suppress plots when generating 1Hz files
+# Dec 9, 2017: - added $antsSuppressCommonOptions = 1;
+# Dec 17, 2017: - added dependencies
+# Mar 8, 2018: - made error messages work with -v
+# - renamed -l)owpass option to -c)utoff
+# - added new -l)atitude
+# Mar 9, 2018: - changed -l) to location (include lon as well)
+# May 22, 2018: - added code to remove initial at surface interval
+# Jul 24, 2018: - improved CTD load error message
+# Oct 4, 2018: - minor improvement in log
+# - BUG: profiles collected with LOTS of wave heave (AAIW 2005, 004)
+# failed the contiguous pressure test at 2dbar resolution
+# => increased to 3dbar
+# Oct 5, 2018: - BUG: previous bug "fix" did not solve a real problem, because
+# problem with AAIW data is lack of pressure resolution
+# => returned to 2dbar
+# - added plotting errors
+# - improved log message
+# - BUG: initial in-air scans were not handled correctly (nscans not updated)
+
# NOTES:
# w_CTD is positive during the downcast to make the sign of the apparent
@@ -76,39 +97,49 @@
require "$ANTS/ants.pl";
require "$ANTS/fft.pl";
require "$ANTS/libstats.pl";
+require "$ANTS/libconv.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:w:',1,
+$antsSuppressCommonOptions = 1;
+&antsUsage('ai:l:orp:qs: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_CTD <cutoff[2s]>] [-w)inch-speed <granularity[10s]>]',
- '[profile -i)d <id>]',
+ '[lowpass w_CTD -c)utoff <limit[2s]>] [-w)inch-speed <granularity[10s]>]',
+ '[profile -i)d <id>] [station -l)ocation <lat/lon>]',
'[-p)lot_basenames <[%03d_w_CTD.ps],[%03d_sspd.ps]>]',
+ '[-q)uiet (no plots)]',
'<SBE CNV file>');
-&antsFloatOpt(\$opt_l,2); # default low-pass cutoff for w_CTD
-&antsCardOpt(\$opt_s,6); # default output samplint rate (Hz)
+&antsFloatOpt(\$opt_c,2); # default low-pass cutoff for w_CTD
+&antsCardOpt(\$opt_s,6); # default output sampling rate (Hz)
&antsFloatOpt(\$opt_w,10); # winch velocity granularity
&antsCardOpt(\$opt_v,$ENV{VERB}); # support VERB env variable
$CNVfile = $ARGV[0]; # open CNV file
open(F,&antsFileArg());
+&antsAddDeps($CNVfile);
&antsActivateOut(); # activate ANTS file
#----------------------------------------------------------------------
# Read Data
#----------------------------------------------------------------------
+sub _croak(@)
+{
+ print(STDERR "\n") if ($opt_v);
+ croak(@_);
+}
+
print(STDERR "Reading $CNVfile...") if ($opt_v);
chomp($rec = <F>);
-croak("$CNVfile: no data\n")
+_croak("$CNVfile: no data\n")
unless defined($rec);
if ($rec =~ /^\*/) { # SBE CNV file
@@ -116,10 +147,21 @@
($nfields,$nscans,$sampint,$badval,$ftype,$lat,$lon) = # decode SBE header
SBE_parseHeader(F,0,0); # SBE field names, no time check
- croak("$CNVfile: unexpected sampling interval $sampint\n")
+ _croak("$CNVfile: unexpected sampling interval $sampint\n")
unless (abs($sampint-1/24) < 1e-5);
- croak("$CNVfile: no latitude in header\n")
+
+ if (defined($opt_l)) { # set/override station location with -l
+ my($slat,$slon) = split('[,/]',$opt_l);
+ $lat = GMT2deg($slat);
+ $lon = GMT2deg($slon);
+ _croak("$0: cannot decode -l $opt_l\n")
+ unless numberp($lat) && numberp($lon);
+ }
+ _croak("$CNVfile: no latitude in header => -l required\n")
unless numberp($lat);
+
+ &antsAddParams('lat',$lat);
+ &antsAddParams('lon',$lon);
$pressF = fnr('prDM');
@@ -151,15 +193,15 @@
if (@ants_ != $nscans) {
if ($opt_v > 1) {
- printf(STDERR "\n\nWARNING: $CNVfile has wrong number of scans in header\n");
+ printf(STDERR "\n\nWARNING: $CNVfile has wrong number of scans in header (%d instead of %d)\n",$nscans,scalar(@ants_));
} else {
- printf(STDERR "WARNING: $CNVfile has wrong number of scans in header\n");
+ printf(STDERR "WARNING: $CNVfile has wrong number of scans in header (%d instead of %d)\n",$nscans,scalar(@ants_));
}
$nscans = @ants_;
}
} else { # simple CSV ASCII file format:
($lat,$lon,$station) = split(',',$rec);
- croak("$CNVfile: ASCII file format error (1st rec must be lat,lon[,id])\n") # header: lat,lon[,station]
+ _croak("$CNVfile: ASCII file format error (1st rec must be lat,lon[,id])\n") # header: lat,lon[,station]
unless numberp($lat) && numberp($lon) &&
$lat>=-90 && $lat<=90 &&
$lon>=-360 && $lon<=360;
@@ -171,7 +213,7 @@
for ($nscans=0; <F>; $nscans++) {
chomp;
@{$ants_[$nscans]} = split(','); # CSV format
- croak("$CNVfile: unexpected scan #$ants_[$nscans][0] (%d expected)\n",$nscans+1)
+ _croak(sprintf("$CNVfile: unexpected scan #$ants_[$nscans][0] (%d expected)\n",$nscans+1))
unless ($ants_[$nscans][0] == $nscans+1);
$ants_[$nscans][$pressF] = nan unless defined($ants_[$nscans][$pressF]); # missing values
$ants_[$nscans][$tempF] = nan unless defined($ants_[$nscans][$tempF]);
@@ -179,7 +221,8 @@
}
}
-printf(STDERR "\n\t%d scans",$nscans) if ($opt_v > 1);
+printf(STDERR "\n\t%d scans (%d seconds)",$nscans,int($nscans*$sampint))
+ if ($opt_v > 1);
printf(STDERR "\n") if ($opt_v);
#----------------------------------------------------------------------
@@ -187,9 +230,9 @@
#----------------------------------------------------------------------
$id = defined($opt_i) ? $opt_i : &antsParam('station');
-croak("$CNVfile: no station information in header => -i required\n")
+_croak("$CNVfile: no station information in header => -i required\n")
unless defined($id);
-croak("$CNVfile: non-numeric station information <$id> in header => -i required\n")
+_croak("$CNVfile: non-numeric station information <$id> in header => -i required\n")
unless numberp($id);
if (-t STDOUT) {
@@ -205,6 +248,8 @@
open(STDOUT,">$outfile") || die("$outfile: $!\n");
}
+undef($opt_p) if $opt_q; # suppress all plots on -q
+
#----------------------------------------------------------------------
# Edit Data
# - pressure outliers & spikes
@@ -227,7 +272,7 @@
print(STDERR "Editing Data...") if ($opt_v);
#----------------------------------------
- # trim initial records with
+ # trim initial scans with
# - nan pressure
# - nan conductivity
# - conductivity <= 10 mS/cm
@@ -240,36 +285,38 @@
numberp($ants_[0][$condF]) &&
(($P{'cond.unit'} eq 'mS/cm' && $ants_[0][$condF] > 10) ||
($P{'cond.unit'} eq 'S/m' && $ants_[0][$condF] > 1));
- croak("\n$CNVfile: no valid records (wrong conductivity units?)\n")
- unless (@ants_);
- printf(STDERR "\n\t%d initial at-surface records trimmed",$trimmed) if ($opt_v > 1);
+ _croak("\n$CNVfile: no valid scans (wrong conductivity units?)\n")
+ unless (@ants_);
+ $nscans -= $trimmed;
+ printf(STDERR "\n\t%d initial in-air scans trimmed",$trimmed) if ($opt_v > 1);
my($lvp) = $ants_[0][$pressF];
my($lvc) = $ants_[0][$condF];
- #------------------------------------------------
+ #-------------------------------------------------------------------------
# edit pressure outliers outside contiguous range
- # - 2dbar resolution
+ # - 2dbar resolution increased to 3dbar because of 2005 AAIW profile 004
# - histogram shifted by 100dbar to allow for negative values
- #------------------------------------------------
+ #-------------------------------------------------------------------------
+ my($press_rez) = 2;
my($outliers) = my($modeSamp) = 0; my($modeBin,$min,$max); local(@hist);
for (my($s)=0; $s<$nscans; $s++) {
next unless ($ants_[$s][$pressF]>=-100 && $ants_[$s][$pressF]<=6500);
- my($b) = ($ants_[$s][$pressF]+100) / 2;
+ my($b) = ($ants_[$s][$pressF]+100) / $press_rez;
$hist[$b]++;
next unless ($hist[$b] > $modeSamp);
$modeSamp = $hist[$b]; $modeBin = $b;
}
- printf(STDERR "\n\tvalid pressure guess: %d dbar (%d samples)",2*$modeBin-100,$modeSamp)
+ printf(STDERR "\n\tvalid pressure guess: %d dbar (%d samples)",$press_rez*$modeBin-100,$modeSamp)
if ($opt_v > 1);
($min,$max) = validRange($modeBin);
- $min = 2*$min-100; $max = 2*$max-100;
+ $min = $press_rez*$min-100; $max = $press_rez*$max-100;
for (my($s)=0; $s<$nscans; $s++) {
if ($ants_[$s][$pressF] > $max) { $outliers++; $ants_[$s][$pressF] = nan; }
if ($ants_[$s][$pressF] < $min) { $outliers++; $ants_[$s][$pressF] = nan; }
}
&antsAddParams("pressure_outliers",sprintf("%d",$outliers));
printf(STDERR "\n\tcontinuous pressure range: %d..%d dbar (%d outliers removed)",$min,$max,$outliers) if ($opt_v > 1);
- croak("$CNVfile: pressure editing removed too many 'outliers'\n")
+ _croak("$CNVfile: pressure editing removed too many 'outliers'\n")
unless ($opt_o || $outliers < 100);
#----------------------------------------------------
@@ -292,7 +339,7 @@
}
&antsAddParams("conductivity_outliers",sprintf("%d",$outliers));
printf(STDERR "\n\tcontinuous conductivity range: %.1f..%.1f S/m (%d outliers removed)",$min,$max,$outliers) if ($opt_v > 1);
- croak("$CNVfile: conductivity editing removed too many 'outliers'\n")
+ _croak("$CNVfile: conductivity editing removed too many 'outliers'\n")
unless ($opt_o || $outliers/$nscans < 0.4);
#----------------------------------------------------
@@ -324,7 +371,7 @@
&antsAddParams("temperature_outliers",sprintf("%d",$outliers));
printf(STDERR "\n\tcontinuous temperature range: %.1f..%.1f degC (%d outliers removed)",$min,$max,$outliers)
if ($opt_v > 1);
- croak("$CNVfile: temperature editing removed too many 'outliers'\n")
+ _croak("$CNVfile: temperature editing removed too many 'outliers'\n")
unless ($opt_o || $outliers/$nscans < 0.4);
#----------------------------------------
@@ -359,7 +406,7 @@
}
}
&antsAddParams("pressure_spikes_removed",sprintf("%d+%d",$ns1,$ns2));
- printf(STDERR "\n\t%d+%d pressure spikes removed",$ns1,$ns2) if ($opt_v > 1);
+ printf(STDERR "\n\t%d+%d pressure spikes removed",$ns1,$ns2) if ($opt_v>1 && $ns1+$ns2>0);
#--------------------------------------------------
# edit conductivity spikes based on large gradients
@@ -436,7 +483,7 @@
$minP = $ants_[$s][$pressF]
if numberp($ants_[$s][$pressF]) && ($ants_[$s][$pressF] < $minP);
}
-croak("$CNVfile: no valid CTD pressure data below 25dbar\n")
+_croak("$CNVfile: no valid CTD pressure data below 25dbar\n")
unless ($minP < 9e99);
if ($minP < 25) {
@@ -452,6 +499,27 @@
printf(STDERR "\n") if ($opt_v);
#----------------------------------------------------------------------
+# Removing Initial At-Surface Data
+#----------------------------------------------------------------------
+
+print(STDERR "Removing intial at-surface data...") if ($opt_v);
+
+my($trimmed) = 0;
+while ($ants_[$trimmed][$pressF] < 0.5) { $trimmed++; }
+for (my($r)=$trimmed; $r<$nscans; $r++) {
+ $ants_[$r][$elapsedF] -= $ants_[$trimmed][$elapsedF];
+}
+splice(@ants_,0,$trimmed);
+$nscans -= $trimmed;
+
+&antsAddParams('surface_data_trimmed',int($trimmed*$sampint));
+&antsAddParams('cast_duration',int($nscans*$sampint));
+
+printf(STDERR "\n\t%d seconds of data trimmed",int($trimmed*$sampint)) if ($opt_v > 1);
+
+printf(STDERR "\n") if ($opt_v);
+
+#----------------------------------------------------------------------
# Binning data
#----------------------------------------------------------------------
@@ -524,9 +592,9 @@
# - interpolate missing vertical velocities first
#----------------------------------------------------------------------
-if ($opt_l > 0) {
+if ($opt_c > 0) {
print(STDERR "Low-pass filtering vertical package velocity...") if ($opt_v);
- &antsAddParams('w_lowpass_cutoff',$opt_l);
+ &antsAddParams('w_lowpass_cutoff',$opt_c);
my($trimmed) = 0;
shift(@w),shift(@depth),shift(@elapsed),shift(@sspd),$trimmed++
@@ -563,7 +631,7 @@
$fftbuf[2*$r] = $w[$r];
$fftbuf[2*$r+1] = 0;
}
- printf(STDERR "\n\t%d zero records added",$pot-$r) if ($opt_v > 1);
+ printf(STDERR "\n\tzero-padded %d scans",$pot-$r) if ($opt_v > 1);
while ($r < $pot) { # pad with zeroes
$fftbuf[2*$r] = $fftbuf[2*$r+1] = 0;
$r++;
@@ -579,7 +647,7 @@
my($in) = 2*$n-$ip; # -ve freq fco
my($f) = $ip/2/$n*$opt_s; # frequency
$fco[$ip] = $fco[$ip+1] = $fco[$in] = $fco[$in+1] = 0
- if ($f > 1/$opt_l); # low-pass filter
+ if ($f > 1/$opt_c); # low-pass filter
}
@w_lp = &FOUR1(1,@fco); # inverse FFT
@@ -622,13 +690,16 @@
if (defined($opt_p)) {
print(STDERR "Plotting data...\n") if ($opt_v);
my(@pfmt) = split(',',$opt_p);
- croak("$0: cannot decode -p $opt_p\n")
+ _croak("$0: cannot decode -p $opt_p\n")
unless (@pfmt == 2);
my($xmin) = $elapsed[0]/60;
my($xmax) = $elapsed[$#elapsed]/60;
my($ymin) = -3; my($ymax) = 3;
my($plotsize) = '13c';
+
+ _croak(sprintf("%s: invalid region of interest (-R$xmin/$xmax/$ymin/$ymax)\n",sprintf($pfmt[0],$id)))
+ unless ($xmax > $xmin && $ymax > $ymin);
GMT_begin(sprintf($pfmt[0],$id),"-JX${plotsize}","-R$xmin/$xmax/$ymin/$ymax",'-X6 -Y4 -P');
GMT_psxy('-W1,coral');
for ($r=0; $r<@w; $r++) {
@@ -653,6 +724,8 @@
my($xmax) = round($max_sspd+3,5);
my($ymin) = 0; my($ymax) = round($depth[$atBtm]+70,100);
my($plotsize) = '13c';
+ _croak(sprintf("%s: invalid region of interest (-R$xmin/$xmax/$ymin/$ymax)\n",sprintf($pfmt[1],$id)))
+ unless ($xmax > $xmin && $ymax > $ymin);
GMT_begin(sprintf($pfmt[1],$id),"-JX${plotsize}/-${plotsize}","-R$xmin/$xmax/$ymin/$ymax",'-X6 -Y4 -P');
GMT_psbasemap('-Bg10a10f1:"Speed of Sound [m/s]":/g1000a500f100:"Depth [m]":WeSn');
GMT_psxy('-W2,coral');
Binary file LADCP_w_howto.pdf has changed
--- a/LADCP_w_ocean Thu Mar 16 11:53:27 2017 -0400
+++ b/LADCP_w_ocean Tue Nov 27 16:59:05 2018 -0500
@@ -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 Mar 6 13:46:31 2017
+# dlm: Tue Nov 27 14:04:39 2018
# (c) 2010 A.M. Thurnherr
-# uE-Info: 280 0 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 282 15 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# TODO:
@@ -277,6 +277,22 @@
# Dec 22, 2016: - moved $opt_p to [defaults.pl]
# Dec 23, 2016: - BUG: -u did not set required variables to proceed
# Mar 6, 2017: - BUG: division by zero when water-depth ~ max(CTD_depth)
+# Oct 12, 2017: - BUG: beampair w did not work for earth-coord vels; major re-write
+# of earthcoord code to unify processing
+# Nov 26, 2017: - BUG: $bad_beam did not work correctly with bin interpolation
+# - BUG: ping-coherent residual removal did not respect missing values
+# Nov 28, 2017: - added $initial_time_lag
+# - expanded semantics of -q to disable time-lagging and residual filters
+# Dec 9, 2017: - added $antsSuppressCommonOptions = 1;
+# Dec 17, 2017: - added dependencies
+# Apr 24, 2018: - added support for $water_depth_db_cmd
+# May 1, 2018: - added threshold for reference-layer horizontal speed
+# - added ambiguity velocity check
+# May 2, 2018: - BUG: ref-lr threshold did not work
+# - BUG: BT code was called for UL when -h was used
+# - replaced $PPI_seabed_editing_required by &PPI_seabed_editing_required
+# - BUG: surface PPI editing code could not be enabled; added &PPI_surface_editing_required
+# Nov 2, 2018: - BUG: for 3-beam solutions, residual{12,34} with affected beam was wrong
# HISTORY END
# CTD REQUIREMENTS
@@ -371,13 +387,13 @@
require "$WCALC/defaults.pl"; # load default/option parameters
$antsParseHeader = 0;
+$antsSuppressCommonOptions = 1;
&antsUsage('3:4a:b:c:de:g:h:i:k:lm:n:o:p:qr:s:t:uv:Vw:x:',0,
"[print software -V)ersion] [-v)erbosity <level[$opt_v]>]",
- "[-q)uick (no single-ping denoising)]",
"[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[,delta.max][$opt_r m/s]>]",
+ "[max -r)esidual <rms.max[,delta.max][$opt_r m/s]>]",
"[-h water <depth|filename>]",
"[max LADCP time-series -g)ap <length[$opt_g s]>]",
"[-i)nitial CTD time offset <guestimate> [-u)se as final]]",
@@ -387,13 +403,14 @@
"[pressure-sensor -a)cceleration-derivative correction <residual/CTD_w_tt>]",
"[-o)utput bin <resolution[$opt_o m]>] [-k) require <min[$opt_k]> samples]",
"[e-x)ecute <perl-expr>]",
+ "[-q)uick-and-dirty (no single-ping denoising, residual and time-lagging filters)]",
"<profile-id> [run-label]");
if ($opt_V) {
- printf(STDERR "+-------------------------+\n");
- printf(STDERR "| LADCP_w Software V%s |\n",$VERSION);
- printf(STDERR "|(c) 2015- A.M. Thurnherr |\n");
- printf(STDERR "+-------------------------+\n");
+ printf(STDERR "+------------------------------+\n");
+ printf(STDERR "| LADCP_w Software V%s |\n",$VERSION);
+ printf(STDERR "| (c) 2015-2017 A.M. Thurnherr |\n");
+ printf(STDERR "+------------------------------+\n");
exit(0);
}
@@ -433,9 +450,17 @@
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 .= ' -l' if defined($opt_l);
-$processing_options .= ' -q' if defined($opt_q);
+
+if (defined($opt_q)) { # quick-and-dirty
+ $processing_options .= ' -q';
+ $opt_l = 1; # disable time-lagging filter
+}
+
+if (defined($opt_i)) { # set initial time lag
+ $processing_options .= " -i $opt_i";
+ $initial_time_lag = &antsFloatOpt($opt_i);
+}
if (defined($opt_x)) { # eval cmd-line expression to override anything
$processing_options .= " -x $opt_x";
@@ -447,7 +472,7 @@
$RDI_Coords::minValidVels = 4;
}
-if ($opt_d) {
+if ($opt_d) { # disable bin mapping
$processing_options .= ' -d';
$RDI_Coords::binMapping = 'none';
}
@@ -478,7 +503,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('RDI_Coords::binMapping',$RDI_Coords::binMapping); # must be set below since binmapping depends on coords
&antsAddParams('processing_options',$processing_options);
&antsAddParams('out_basename',$out_basename);
&antsAddParams('profile_id',$PROF,'run_label',$RUN);
@@ -551,12 +576,26 @@
progress("Reading LADCP data from <$LADCP_file>...\n");
readData($LADCP_file,\%LADCP);
+&antsAddDeps($LADCP_file);
if ($LADCP{BEAM_COORDINATES}) {
progress("\t%d ensembles (beam coordinates)\n",scalar(@{$LADCP{ENSEMBLE}}));
} else {
progress("\t%d ensembles (Earth coordinates)\n",scalar(@{$LADCP{ENSEMBLE}}));
}
+if ($valid_ensemble_range[0] > 0) { # remove leading invalid records
+ my($ens) = 0;
+ while ($ens < @{$LADCP{ENSEMBLE}} && $LADCP{ENSEMBLE}[$ens]->{NUMBER}<$valid_ensemble_range[0]) { $ens++ }
+ splice(@{$LADCP{ENSEMBLE}},0,$ens);
+ progress("\t%d invalid leading ensembles removed\n",$ens)
+}
+if ($valid_ensemble_range[1] > 0) { # remove trailing invalid records
+ my($ens) = 0;
+ while ($ens < @{$LADCP{ENSEMBLE}} && $LADCP{ENSEMBLE}[$ens]->{NUMBER}<$valid_ensemble_range[1]) { $ens++ }
+ splice(@{$LADCP{ENSEMBLE}},$ens);
+ progress("\t%d invalid trailing ensembles removed\n",@{$LADCP{ENSEMBLE}}-$ens)
+}
+
error("$LADCP_file: cannot process multi-ping ensembles\n")
unless ($LADCP{PINGS_PER_ENSEMBLE} == 1);
warning(2,"$LADCP_file: wide-bandwidth setting\n")
@@ -598,11 +637,15 @@
'ADCP_frequency',$LADCP{BEAM_FREQUENCY},
'ADCP_blanking_distance',$LADCP{BLANKING_DISTANCE});
+
#------------------------------------------------------------
-# Edit beam-velocity data
-# 0) beam-vel mask on -m
-# mask file has three columns: from_ens to_ens ignore_beam
-# 1) correlation threshold
+# Edit velocity data
+# beam coords:
+# 1) beam-vel mask on -m
+# mask file has three columns: from_ens to_ens ignore_beam
+# 2) correlation threshold
+# Earth coords:
+# 1) correlation threshold
#------------------------------------------------------------
if ($LADCP{BEAM_COORDINATES}) {
@@ -663,6 +706,27 @@
progress("\tcorrelation threshold (-c %d counts): %d velocites removed (%d%% of total)\n",$opt_c,$cte,round(100*$cte/$nvv));
}
+#------------------------------------------------------------
+# Create beam coordinate velocities for Earth-velocity data
+# - velocities are replaced "in place"
+# - BT velocities are treated separately in [find_seabed.pl]
+# - this transformation will remove all 3-beam solutions
+# - disable bin mapping because Earth coords are typically bin-remapped
+#------------------------------------------------------------
+
+unless ($LADCP{BEAM_COORDINATES}) {
+ progress("Replacing earth- with beam-velocities...\n");
+ for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
+ for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
+ @{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]} =
+ &velEarthToBeam(\%LADCP,$ens,@{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]});
+ }
+ }
+ $RDI_Coords::binMapping = 'none';
+}
+
+&antsAddParams('RDI_Coords::binMapping',$RDI_Coords::binMapping); # finally, bin mapping is known
+
#-------------------------------------------------------------------
# Calculate earth velocities
# - this is done for all bins (not just valid ones), to allow
@@ -673,85 +737,59 @@
# been very useful so far)
#-------------------------------------------------------------------
-if ($LADCP{BEAM_COORDINATES}) {
- my($dummy);
- progress("Calculating earth-coordinate velocities...\n");
- if ($bad_beam) {
- progress("\tdiscarding velocities from beam $bad_beam\n");
- &antsAddParams('bad_beam_discarded',$bad_beam);
- }
- $nvw = 0;
- for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
- $LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH} =
- gimbal_pitch($LADCP{ENSEMBLE}[$ens]->{PITCH},$LADCP{ENSEMBLE}[$ens]->{ROLL});
+my($dummy);
+progress("Calculating earth-coordinate velocities...\n");
+if ($bad_beam) {
+ progress("\tdiscarding velocities from beam $bad_beam\n");
+ &antsAddParams('bad_beam_discarded',$bad_beam);
+}
+$nvw = 0;
+for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
+ $LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH} =
+ gimbal_pitch($LADCP{ENSEMBLE}[$ens]->{PITCH},$LADCP{ENSEMBLE}[$ens]->{ROLL});
- 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]->{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);
+ if ($bad_beam) {
+ for (my($bin)=0; $bin<=$LADCP{N_BINS}-1; $bin++) {
+ undef($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$bad_beam-1]);
+ undef($LADCP{ENSEMBLE}[$ens]->{BT_VELOCITY}[$bin][$bad_beam-1]);
}
+ }
- for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
- $LADCP{ENSEMBLE}[$ens]->{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]);
+ for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
+ ($LADCP{ENSEMBLE}[$ens]->{INTERP_U}[$bin],
+ $LADCP{ENSEMBLE}[$ens]->{INTERP_V}[$bin],
+ $LADCP{ENSEMBLE}[$ens]->{INTERP_W}[$bin],
+ $LADCP{ENSEMBLE}[$ens]->{INTERP_ERRVEL}[$bin]) = earthVels(\%LADCP,$ens,$bin);
+ ($LADCP{ENSEMBLE}[$ens]->{V12}[$bin],$LADCP{ENSEMBLE}[$ens]->{W12}[$bin],
+ $LADCP{ENSEMBLE}[$ens]->{V34}[$bin],$LADCP{ENSEMBLE}[$ens]->{W34}[$bin]) = BPEarthVels(\%LADCP,$ens,$bin);
+ }
- if (defined($LADCP{ENSEMBLE}[$ens]->{W}[$bin])) {
- $per_bin_nsamp[$bin]++;
- $nvw++;
- }
-
- $LADCP{ENSEMBLE}[$ens]->{U_UNEDITED}[$bin] = $LADCP{ENSEMBLE}[$ens]->{U}[$bin];
- $LADCP{ENSEMBLE}[$ens]->{V_UNEDITED}[$bin] = $LADCP{ENSEMBLE}[$ens]->{V}[$bin];
- $LADCP{ENSEMBLE}[$ens]->{W_UNEDITED}[$bin] = $LADCP{ENSEMBLE}[$ens]->{W}[$bin];
+ for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
+ $LADCP{ENSEMBLE}[$ens]->{U}[$bin] = $LADCP{ENSEMBLE}[$ens]->{INTERP_U}[$bin];
+ $LADCP{ENSEMBLE}[$ens]->{V}[$bin] = $LADCP{ENSEMBLE}[$ens]->{INTERP_V}[$bin];
+ $LADCP{ENSEMBLE}[$ens]->{W}[$bin] = $LADCP{ENSEMBLE}[$ens]->{INTERP_W}[$bin];
+ $LADCP{ENSEMBLE}[$ens]->{ERRVEL}[$bin] = $LADCP{ENSEMBLE}[$ens]->{INTERP_ERRVEL}[$bin];
+ undef($LADCP{ENSEMBLE}[$ens]->{INTERP_U}[$bin]);
+ undef($LADCP{ENSEMBLE}[$ens]->{INTERP_V}[$bin]);
+ undef($LADCP{ENSEMBLE}[$ens]->{INTERP_W}[$bin]);
+ undef($LADCP{ENSEMBLE}[$ens]->{INTERP_ERRVEL}[$bin]);
+
+ if (defined($LADCP{ENSEMBLE}[$ens]->{W}[$bin])) {
+ $per_bin_nsamp[$bin]++;
+ $nvw++;
}
+
+ $LADCP{ENSEMBLE}[$ens]->{U_UNEDITED}[$bin] = $LADCP{ENSEMBLE}[$ens]->{U}[$bin];
+ $LADCP{ENSEMBLE}[$ens]->{V_UNEDITED}[$bin] = $LADCP{ENSEMBLE}[$ens]->{V}[$bin];
+ $LADCP{ENSEMBLE}[$ens]->{W_UNEDITED}[$bin] = $LADCP{ENSEMBLE}[$ens]->{W}[$bin];
}
- progress("\t$nvw valid velocities in bins $LADCP_firstBin-$LADCP_lastBin\n");
- progress("\t3-beam solutions : $RDI_Coords::threeBeam_1 " .
- "$RDI_Coords::threeBeam_2 " .
- "$RDI_Coords::threeBeam_3 " .
- "$RDI_Coords::threeBeam_4\n")
- unless ($opt_4);
-} else { # Earth Coordinates
- 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],
- $LADCP{ENSEMBLE}[$ens]->{W}[$bin],
- $LADCP{ENSEMBLE}[$ens]->{ERRVEL}[$bin]) = @{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]};
- if (defined($LADCP{ENSEMBLE}[$ens]->{W}[$bin])) {
- $per_bin_nsamp[$bin]++;
- $nvw++;
- }
- $LADCP{ENSEMBLE}[$ens]->{U_UNEDITED}[$bin] = $LADCP{ENSEMBLE}[$ens]->{U}[$bin];
- $LADCP{ENSEMBLE}[$ens]->{V_UNEDITED}[$bin] = $LADCP{ENSEMBLE}[$ens]->{V}[$bin];
- $LADCP{ENSEMBLE}[$ens]->{W_UNEDITED}[$bin] = $LADCP{ENSEMBLE}[$ens]->{W}[$bin];
-
- $LADCP{ENSEMBLE}[$ens]->{V12}[$bin] = $LADCP{ENSEMBLE}[$ens]->{V34}[$bin] = nan;
-
- # for 3-beam solutions, w12 = w34 = w (I think)
- ($LADCP{ENSEMBLE}[$ens]->{W12}[$bin],$LADCP{ENSEMBLE}[$ens]->{W34}[$bin]) =
- velEarthToBPw($LADCP,$ens,@{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]});
- }
- }
- progress("\t$nvw valid velocities in bins $LADCP_firstBin-$LADCP_lastBin\n");
}
+progress("\t$nvw valid velocities in bins $LADCP_firstBin-$LADCP_lastBin\n");
+progress("\t3-beam solutions : $RDI_Coords::threeBeam_1 " .
+ "$RDI_Coords::threeBeam_2 " .
+ "$RDI_Coords::threeBeam_3 " .
+ "$RDI_Coords::threeBeam_4\n")
+ unless ($opt_4);
error("$LADCP_file: insufficient valid velocities\n") unless ($nvw >= $min_valid_vels);
@@ -856,7 +894,9 @@
#----------------------------------------------------------------------
# More editing
-# - this requires ${first,last}GoodEns to be known
+# - attitude threshold
+# - data in far bins (beyond reliable range)
+# - at this stage ${first,last}GoodEns are known
# - TILT field is set as a side-effect
#----------------------------------------------------------------------
@@ -896,15 +936,16 @@
progress("\tvelocities beyond bin $first_bad_bin (<%d%% valid values): %d velocites removed (%d%% of total in bins $LADCP_firstBin-$LADCP_lastBin)\n",
round(100*$per_bin_valid_frac_lim),$fprm,round(100*$fprm/$nvw));
-#--------------
+#----------------------------------------------------------------------
# Read CTD data
-#--------------
+#----------------------------------------------------------------------
progress("Reading CTD data from <$CTD_file>...\n");
open(STDIN,$CTD_file) || error("$CTD_file: $!\n");
error("$CTD_file: no data\n") unless (&antsIn());
undef($antsOldHeaders);
+&antsAddDeps($CTD_file);
&antsAddParams('lat',$P{lat}) if defined($P{lat});
&antsAddParams('lon',$P{lon}) if defined($P{lon});
@@ -978,9 +1019,9 @@
# Determine time lag
#-------------------
-if (defined($opt_i)) {
+if (defined($initial_time_lag)) {
progress("Setting initial time lag...\n");
- $CTD{TIME_LAG} = $opt_i;
+ $CTD{TIME_LAG} = $initial_time_lag;
progress("\t-i => elapsed(CTD) ~ elapsed(LADCP) + %.1fs\n",$CTD{TIME_LAG});
} else {
progress("Guestimating time lag...\n");
@@ -1191,10 +1232,10 @@
# 2) PPI
#----------------------------------------------------------------------------
-error("$0: conflicting water-depth information provided by user\n") # this can only happen if
- if defined($opt_h) && defined($water_depth); # the user is setting $water_depth
- # explicitly?!
-if (defined($opt_h)) { # deal with user-provided water-depth info
+error("$0: conflicting water-depth information provided by user\n") # only happens when $water_depth is set explicitly
+ if defined($opt_h) && defined($water_depth);
+
+if (defined($opt_h)) { # handle user-provided water-depth info
if (numberp($opt_h)) {
$water_depth = $opt_h;
} elsif (-f $opt_h) {
@@ -1207,9 +1248,6 @@
}
}
-#die("assertion failed (water_depth = $water_depth)")
-# unless (!defined($water_depth) || numberp($water_depth));
-
if (!defined($water_depth) && # find seabed in data
$LADCP{ENSEMBLE}[$LADCP_atbottom]->{XDUCER_FACING_DOWN}) {
progress("Finding seabed...\n");
@@ -1226,11 +1264,11 @@
warning(1,"using water_depth from ADCP BT data\n"); # explicitly requested
$SS_use_BT = 1;
}
- if ($SS_use_BT) { # water depth from BT data
+ if ($SS_use_BT && numberp($water_depth_BT)) { # water depth from BT data
&antsAddParams('water_depth_from','BT_data');
$water_depth = $water_depth_BT;
$sig_water_depth = $sig_water_depth_BT;
- } else { # water depth from WT data
+ } elsif (defined($water_depth)) { # water depth from WT data
&antsAddParams('water_depth_from','echo_amplitudes');
}
}
@@ -1252,6 +1290,15 @@
}
}
+if (!defined($water_depth) && defined($water_depth_db_cmd)) { # set water depth from data base
+ error("$0: lat/lon required for running $water_depth_db_cmd\n")
+ unless numbersp($P{lat},$P{lon});
+ chomp($water_depth = `$water_depth_db_cmd $P{lon} $P{lat}`);
+ error("$0: command '$water_depth_db_cmd $P{lon} $P{lat}' did not return valid water depth\n")
+ unless numberp($water_depth);
+ &antsAddParams('water_depth_from',"$water_depth_db_cmd $P{lon} $P{lat}");
+}
+
if (defined($water_depth)) { # set %PARAMs
&antsAddParams('water_depth',$water_depth,'water_depth.sig',$sig_water_depth);
} else {
@@ -1276,8 +1323,7 @@
&antsAddParams('sidelobe_editing','seabed');
}
- $PPI_editing |= eval($PPI_seabed_editing_required);
- if ($PPI_editing) {
+ if (&PPI_seabed_editing_required()) {
&antsAddParams('PPI_editing','seabed');
&antsAddParams('PPI_extend_upper_limit',$PPI_extend_upper_limit)
if numberp($PPI_extend_upper_limit);
@@ -1323,7 +1369,7 @@
} else {
&antsAddParams('sidelobe_editing','surface','vessel_draft',$vessel_draft);
}
- if ($PPI_editing) {
+ if (&PPI_surface_editing_required()) {
&antsAddParams('PPI_editing','surface');
&antsAddParams('PPI_extend_upper_limit',$PPI_extend_upper_limit)
if numberp($PPI_extend_upper_limit);
@@ -1341,9 +1387,39 @@
}
#----------------------------------------------------------------------
+# Check For Ambiguity Velocity Problems
+#----------------------------------------------------------------------
+
+progress("Checking for ambiguity velocity violations...\n");
+
+my($ambiguity_velocity) = ambiguity_velocity($LADCP{BEAM_FREQUENCY},
+ $LADCP{BEAM_ANGLE},
+ $LADCP{SPEED_OF_SOUND},
+ $LADCP{TRANSMIT_LAG_DISTANCE});
+&antsAddParams('ambiguity_velocity',$ambiguity_velocity);
+my($nbad) = 0;
+for (my($ens)=$firstGoodEns; $ens<=$lastGoodEns; $ens++) {
+ next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
+ next unless ($CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}] > $ambiguity_velocity);
+ $nbad++;
+}
+
+my($badf) = $nbad / ($lastGoodEns - $firstGoodEns + 1); # fraction of bad values
+if ($bad > 0.01) { # allow 1% violations before warning is triggered
+ warning(2,"%d ensembles (%d%% of total) have CTD_w > ambiguity velocity of %.1 m/s\n",
+ $nbad,round(100*$badf),$ambiguity_velocity);
+} elsif ($nbad > 0) {
+ info("\t%d ensembles (%d%% of total) have CTD_w > ambiguity velocity of %.1 m/s",
+ $nbad,round(100*$badf),$ambiguity_velocity);
+} else {
+ info("\tnone found\n");
+}
+
+#----------------------------------------------------------------------
# Data Editing after LADCP and CTD data have been merged
# 1) surface layer editing
-# 2) Execute user-supplied $edit_data_hook
+# 2) reference-layer horizontal velocity threshold
+# 3) Execute user-supplied $edit_data_hook
#----------------------------------------------------------------------
progress("Removing data from instrument at surface...\n");
@@ -1351,6 +1427,16 @@
$nerm = editSurfLayer($firstGoodEns,$lastGoodEns,$surface_layer_depth);
progress("\t$nerm ensembles removed\n");
+progress("Removing data collected at large horizontal package speeds...\n");
+$max_hspeed = &max_hspeed(); # defined in [defaults.h]
+&antsAddParams('max_hspeed',$max_hspeed);
+$nerm = editLargeHSpeeds($firstGoodEns,$lastGoodEns,$max_hspeed);
+my($nermf) = $nerm / ($lastGoodEns - $firstGoodEns + 1);
+info("\treference-layer horizontal speed threshold (max_hspeed = %g m/s): %d ensembles removed (%d%% of total)\n",
+ $max_hspeed,$nerm,round(100*$nermf));
+warning(2,"large fraction (%d%%) of samples exceed reference-layer horizontal speed threshold\n",round(100*$nermf))
+ if ($nermf > 0.05);
+
if (defined($post_merge_hook)) {
progress("Executing user-supplied \$post_merge_hook...\n");
&{$post_merge_hook}($firstGoodEns,$lastGoodEns);
@@ -1569,11 +1655,14 @@
for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) { # subtract from ocean velocities
next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin] -=
- $LADCP{ENSEMBLE}[$ens]->{MEDIAN_RESIDUAL_W};
+ $LADCP{ENSEMBLE}[$ens]->{MEDIAN_RESIDUAL_W}
+ if numberp($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin]);
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin] -=
- $LADCP{ENSEMBLE}[$ens]->{MEDIAN_RESIDUAL_W};
+ $LADCP{ENSEMBLE}[$ens]->{MEDIAN_RESIDUAL_W}
+ if numberp($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin]);
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin] -=
- $LADCP{ENSEMBLE}[$ens]->{MEDIAN_RESIDUAL_W};
+ $LADCP{ENSEMBLE}[$ens]->{MEDIAN_RESIDUAL_W}
+ if numberp($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin]);
}
$LADCP{ENSEMBLE}[$ens]->{REFLR_W} -= # NB: this can be nan here
@@ -1587,24 +1676,30 @@
for (my($bi)=0; $bi<=$#{$DNCAST{ENSEMBLE}}; $bi++) { # bin data
for (my($i)=0; $i<@{$DNCAST{W}[$bi]}; $i++) { # code works if MEDIAN_RESIDUAL_W is nan (possible?)
- $DNCAST{W} [$bi][$i] -= $LADCP{ENSEMBLE}[$DNCAST{ENSEMBLE}[$bi][$i]]->{MEDIAN_RESIDUAL_W};
- $DNCAST{W12}[$bi][$i] -= $LADCP{ENSEMBLE}[$DNCAST{ENSEMBLE}[$bi][$i]]->{MEDIAN_RESIDUAL_W};
- $DNCAST{W34}[$bi][$i] -= $LADCP{ENSEMBLE}[$DNCAST{ENSEMBLE}[$bi][$i]]->{MEDIAN_RESIDUAL_W};
+ $DNCAST{W} [$bi][$i] -= $LADCP{ENSEMBLE}[$DNCAST{ENSEMBLE}[$bi][$i]]->{MEDIAN_RESIDUAL_W}
+ if numberp($DNCAST{W}[$bi][$i]);
+ $DNCAST{W12}[$bi][$i] -= $LADCP{ENSEMBLE}[$DNCAST{ENSEMBLE}[$bi][$i]]->{MEDIAN_RESIDUAL_W}
+ if numberp($DNCAST{W12}[$bi][$i]);
+ $DNCAST{W34}[$bi][$i] -= $LADCP{ENSEMBLE}[$DNCAST{ENSEMBLE}[$bi][$i]]->{MEDIAN_RESIDUAL_W}
+ if numberp($DNCAST{W34}[$bi][$i]);
}
- $DNCAST{MEDIAN_W} [$bi] = median(@{$DNCAST{W}[$bi]});
- $DNCAST{MEDIAN_W12}[$bi] = median(@{$DNCAST{W12}[$bi]});
- $DNCAST{MEDIAN_W34}[$bi] = median(@{$DNCAST{W34}[$bi]});
+ $DNCAST{MEDIAN_W} [$bi] = median(@{$DNCAST{W}[$bi]}) if numberp($DNCAST{MEDIAN_W} [$bi]);
+ $DNCAST{MEDIAN_W12}[$bi] = median(@{$DNCAST{W12}[$bi]}) if numberp($DNCAST{MEDIAN_W12}[$bi]);
+ $DNCAST{MEDIAN_W34}[$bi] = median(@{$DNCAST{W34}[$bi]}) if numberp($DNCAST{MEDIAN_W34}[$bi]);
$DNCAST{MAD_W} [$bi] = mad2($DNCAST{MEDIAN_W}[$bi],@{$DNCAST{W}[$bi]});
}
for (my($bi)=0; $bi<=$#{$UPCAST{ENSEMBLE}}; $bi++) {
for (my($i)=0; $i<@{$UPCAST{W}[$bi]}; $i++) {
- $UPCAST{W} [$bi][$i] -= $LADCP{ENSEMBLE}[$UPCAST{ENSEMBLE}[$bi][$i]]->{MEDIAN_RESIDUAL_W};
- $UPCAST{W12}[$bi][$i] -= $LADCP{ENSEMBLE}[$UPCAST{ENSEMBLE}[$bi][$i]]->{MEDIAN_RESIDUAL_W};
- $UPCAST{W34}[$bi][$i] -= $LADCP{ENSEMBLE}[$UPCAST{ENSEMBLE}[$bi][$i]]->{MEDIAN_RESIDUAL_W};
+ $UPCAST{W} [$bi][$i] -= $LADCP{ENSEMBLE}[$UPCAST{ENSEMBLE}[$bi][$i]]->{MEDIAN_RESIDUAL_W}
+ if numberp($UPCAST{W}[$bi][$i]);
+ $UPCAST{W12}[$bi][$i] -= $LADCP{ENSEMBLE}[$UPCAST{ENSEMBLE}[$bi][$i]]->{MEDIAN_RESIDUAL_W}
+ if numberp($UPCAST{W12}[$bi][$i]);
+ $UPCAST{W34}[$bi][$i] -= $LADCP{ENSEMBLE}[$UPCAST{ENSEMBLE}[$bi][$i]]->{MEDIAN_RESIDUAL_W}
+ if numberp($UPCAST{W34}[$bi][$i]);
}
- $UPCAST{MEDIAN_W} [$bi] = median(@{$UPCAST{W}[$bi]});
- $UPCAST{MEDIAN_W12}[$bi] = median(@{$UPCAST{W12}[$bi]});
- $UPCAST{MEDIAN_W34}[$bi] = median(@{$UPCAST{W34}[$bi]});
+ $UPCAST{MEDIAN_W} [$bi] = median(@{$UPCAST{W}[$bi]}) if numberp($UPCAST{MEDIAN_W} [$bi]);
+ $UPCAST{MEDIAN_W12}[$bi] = median(@{$UPCAST{W12}[$bi]}) if numberp($UPCAST{MEDIAN_W12}[$bi]);
+ $UPCAST{MEDIAN_W34}[$bi] = median(@{$UPCAST{W34}[$bi]}) if numberp($UPCAST{MEDIAN_W34}[$bi]);
$UPCAST{MAD_W} [$bi] = mad2($UPCAST{MEDIAN_W}[$bi],@{$UPCAST{W}[$bi]});
}
@@ -1614,7 +1709,7 @@
# remove ensembles with large rms residuals
#----------------------------------------------------------------------
-if (defined($opt_r)) {
+if (defined($opt_r) && !$opt_q) {
progress("Applying residuals filters...\n");
progress("\trms residuals > $residuals_rms_max: ");
@@ -1711,8 +1806,10 @@
next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
my($bi) = $bindepth[$bin]/$opt_o;
next unless numberp($DNCAST{MEDIAN_W}[$bi]);
- push(@{$DNCAST{RESIDUAL12}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin]-$DNCAST{MEDIAN_W}[$bi]);
- push(@{$DNCAST{RESIDUAL34}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin]-$DNCAST{MEDIAN_W}[$bi]);
+ push(@{$DNCAST{RESIDUAL12}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin]-$DNCAST{MEDIAN_W}[$bi])
+ if numberp($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin]);
+ push(@{$DNCAST{RESIDUAL34}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin]-$DNCAST{MEDIAN_W}[$bi])
+ if numberp($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin]);
}
}
for (my($bi)=0; $bi<=$#{$DNCAST{ENSEMBLE}}; $bi++) {
@@ -1728,8 +1825,10 @@
next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
my($bi) = $bindepth[$bin]/$opt_o;
next unless numberp($UPCAST{MEDIAN_W}[$bi]);
- push(@{$UPCAST{RESIDUAL12}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin]-$UPCAST{MEDIAN_W}[$bi]);
- push(@{$UPCAST{RESIDUAL34}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin]-$UPCAST{MEDIAN_W}[$bi]);
+ push(@{$UPCAST{RESIDUAL12}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin]-$UPCAST{MEDIAN_W}[$bi])
+ if numberp($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin]);
+ push(@{$UPCAST{RESIDUAL34}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin]-$UPCAST{MEDIAN_W}[$bi])
+ if numberp($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin]);
}
}
for (my($bi)=0; $bi<=$#{$UPCAST{ENSEMBLE}}; $bi++) {
@@ -1842,7 +1941,8 @@
# Calculate BT-referenced vertical-velocity profile
#--------------------------------------------------
-if (defined($water_depth)) {
+if ($LADCP{ENSEMBLE}[$LADCP_atbottom]->{XDUCER_FACING_DOWN} && defined($water_depth)) {
+
progress("Calculating BT-referenced vertical velocities\n");
calc_BTprof($firstGoodEns,$lastGoodEns,$water_depth,$sig_water_depth);
@@ -1887,6 +1987,12 @@
$of = ">$of" unless ($of =~ /^$|^\s*\|/);
open(STDOUT,$of) || error("$of: $!\n");
undef($antsActiveHeader) unless ($ANTS_TOOLS_AVAILABLE);
+
+ sub res($$)
+ {
+ my($meas,$mean) = @_;
+ return numberp($meas) ? ($meas - $mean) : undef;
+ }
for ($ens=$firstGoodEns; $ens<$LADCP_atbottom; $ens++) { # downcast
next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
@@ -1902,9 +2008,9 @@
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin],
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin],
$LADCP{ENSEMBLE}[$ens]->{V12}[$bin],$LADCP{ENSEMBLE}[$ens]->{V34}[$bin],
- $LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin] - $DNCAST{MEDIAN_W}[$bi],
- $LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin] - $DNCAST{MEDIAN_W}[$bi],
- $LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin] - $DNCAST{MEDIAN_W}[$bi],
+ res($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin],$DNCAST{MEDIAN_W}[$bi]),
+ res($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin],$DNCAST{MEDIAN_W}[$bi]),
+ res($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin],$DNCAST{MEDIAN_W}[$bi]),
$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$CTD{W_t}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$CTD{W_tt}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
@@ -1942,9 +2048,9 @@
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin],
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin],
$LADCP{ENSEMBLE}[$ens]->{V12}[$bin],$LADCP{ENSEMBLE}[$ens]->{V34}[$bin],
- $LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin] - $UPCAST{MEDIAN_W}[$bi],
- $LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin] - $UPCAST{MEDIAN_W}[$bi],
- $LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin] - $UPCAST{MEDIAN_W}[$bi],
+ res($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin],$UPCAST{MEDIAN_W}[$bi]),
+ res($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin],$UPCAST{MEDIAN_W}[$bi]),
+ res($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin],$UPCAST{MEDIAN_W}[$bi]),
$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$CTD{W_t}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$CTD{W_tt}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
--- a/LADCP_w_postproc Thu Mar 16 11:53:27 2017 -0400
+++ b/LADCP_w_postproc Tue Nov 27 16:59:05 2018 -0500
@@ -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: Thu May 26 18:21:11 2016
+# dlm: Thu Nov 1 10:55:34 2018
# (c) 2015 A.M. Thurnherr
-# uE-Info: 101 35 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 71 77 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
$antsSummary = 'edit and re-grid LADCP vertical-velocity samples';
@@ -65,6 +65,11 @@
# Apr 14, 2016: - added profile id to warning messages
# May 24, 2016: - improved plot
# May 26, 2016: - added automatic directory creation on -d
+# Nov 28, 2017: - removed wcorr plot
+# Dec 9, 2017: - added $antsSuppressCommonOptions = 1;
+# Oct 31, 2018: - improved label (no longer explained/residual stddev)
+# Nov 1, 2018: - made layout consistent for dual- and single-head profiles
+
($ANTS) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
($WCALC) = ($0 =~ m{^(.*)/[^/]*$});
@@ -79,6 +84,7 @@
require "$ANTS/libGMT.pl";
&antsAddParams('LADCP_w_postproc::version',$VERSION);
+$antsSuppressCommonOptions = 1;
&antsUsage('b:d:i:k:l:o:p:v:w:',1,
'[profile -i)d <id>]',
'[-o)utput bin <resolution>] [-k) require <min> samples]',
@@ -150,8 +156,10 @@
#----------------------------------------------------------------------
if (-t STDOUT) {
- $opt_p = defined($opt_d) ? "$opt_d/%03d_wprof.ps,$opt_d/%03d_wcorr.ps"
- : '%03d_wprof.ps,%03d_wcorr.ps'
+# $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);
+ $opt_p = defined($opt_d) ? "$opt_d/%03d_wprof.ps" : '%03d_wprof.ps'
unless defined($opt_p);
$outfile = defined($opt_d) ? sprintf('%s/%03d.wprof',$opt_d,$id)
: sprintf('%03d.wprof',$id);
@@ -255,9 +263,8 @@
$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);
+ my($rms) = sqrt($sxx/$n);
+ return ($R,$rms);
}
sub uc_corr(@)
@@ -285,9 +292,8 @@
$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);
+ my($sig) = sqrt(($sxx+$syy)/(2*$n)); # stddev of average w signal from DL & UL
+ return ($R,$sig);
}
#----------------------------------------------------------------------
@@ -309,7 +315,7 @@
if (defined($opt_p)) {
($sumPF,$corrPF) = split(/,/,$opt_p);
croak("$0: cannot decode -p $opt_p\n")
- unless (length($corrPF)>0);
+ unless (length($sumPF)>0);
}
my($R,$R2);
@@ -474,11 +480,9 @@
? $DL_uc_median[$bi] - $UL_uc_median[$bi] : nan;
}
- ($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);
+ ($dc_R,$dc_sig) = &dc_corr(); # correlation statistics
+ ($uc_R,$uc_sig) = &uc_corr();
+ &antsAddParams('dc_R',$dc_R,'uc_R',$uc_R,'dc_w.sig',$dc_sig,'uc_w.sig',$uc_sig);
if (defined($opt_p)) { # plot 2nd-instrument profiles
GMT_psxy('-W1,coral,.');
@@ -499,11 +503,12 @@
# 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',
+ 'dc_w.diff','uc_w.diff'); # DL-UL differences
+
if ($dual_head) { # dual-head output
- @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'); # DL-UL differences
&antsAddParams('profile_id',$id,'lat',$P{lat},'lon',$P{lon}); # selected %PARAMs
&antsAddParams($dayNoP,$dn,'run_label',"$first_label & $P{run_label}");
&antsAddParams('outgrid_dz',$opt_o,'outgrid_minsamp',$opt_k);
@@ -512,10 +517,6 @@
&antsAddParams('ADCP_bin_length',$bin_length,'ADCP_pulse_length',$pulse_length);
&antsAddParams('ADCP_blanking_distance',$blanking_dist);
undef($antsOldHeaders);
-} else {
- @antsNewLayout = ('depth','hab', # single-head output
- '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)")
@@ -538,8 +539,11 @@
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);
+ if ($dual_head) {
+ push(@{$out[$bi]},$dc_diff[$bi],$uc_diff[$bi]);
+ } else {
+ push(@{$out[$bi]},nan,nan);
+ }
&antsOut(@{$out[$bi]});
}
@@ -591,25 +595,21 @@
GMT_unitcoords();
if ($dc_R < 0.3 || !numberp($dc_R)) { # correlation statistics
&antsInfo("WARNING: low dc correlation (r = %.1f) between UL and DL data in profile #$id",$dc_R);
- GMT_pstext('-F+f12,Helvetica,coral+jTL -Gred');
+ 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.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",
- 0.99,0.01,100*$dc_esig,100*$dc_rsig) if ($dual_head);
- if ($dual_head) {
+ printf(GMT "%f %f r = %.1f\n",0.62,0.01,$dc_R);
+ GMT_pstext('-F+f12,Helvetica,coral+jTR -Gwhite');
+ printf(GMT "%f %f [%.1f cm/s stddev]\n",0.99,0.01,100*$dc_sig);
if ($uc_R < 0.3 || !numberp($uc_R)) {
&antsInfo("WARNING: low uc correlation (r = %.1f) between UL and DL data in profile #$id",$uc_R);
- GMT_pstext('-F+f12,Helvetica,SeaGreen+jTL -Gred');
+ 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.61,0.05,$uc_R);
+ printf(GMT "%f %f r = %.1f\n",0.62,0.05,$uc_R);
+ GMT_pstext('-F+f12,Helvetica,SeaGreen+jTR -Gwhite');
+ printf(GMT "%f %f [%.1f cm/s stddev]\n",0.99,0.05,100*$uc_sig);
}
- GMT_pstext('-F+f12,Helvetica,SeaGreen+jTR -Gwhite');
- printf(GMT "%f %f [%.1f/%.1f cm/s @~s@~\@-e/r\@-]\n",
- 0.99,0.05,100*$uc_esig,100*$uc_rsig) if ($dual_head);
GMT_pstext('-F+f14,Helvetica,blue+jTL -N');
if (defined($outfile)) { print(GMT "0.01 -0.06 $outfile [$P{run_label}]\n"); }
@@ -621,7 +621,7 @@
GMT_end();
- if ($dual_head) { # correlation plot
+ if ($dual_head && length($corrPF)>0) { # 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]);
@@ -682,7 +682,7 @@
GMT_psbasemap('-Bf0.01a0.05:"DL Vertical Velocity [m/s]":/f0.01a0.05:"UL Vertical Velocity [m/s]":WeSn');
GMT_end();
- } # if dual_head
+ } # if dual_head && length(corrPF) > 0
}
--- a/LADCP_wspec Thu Mar 16 11:53:27 2017 -0400
+++ b/LADCP_wspec Tue Nov 27 16:59:05 2018 -0500
@@ -2,9 +2,9 @@
#======================================================================
# L A D C P _ W S P E C
# doc: Thu Jun 11 12:02:49 2015
-# dlm: Sun Mar 12 12:01:59 2017
+# dlm: Thu Nov 1 10:09:48 2018
# (c) 2012 A.M. Thurnherr
-# uE-Info: 48 17 NIL 0 0 72 10 2 4 NIL ofnI
+# uE-Info: 39 29 NIL 0 0 72 10 2 4 NIL ofnI
#======================================================================
$antsSummary = 'calculate VKE window spectra from LADCP profiles';
@@ -32,6 +32,11 @@
# - BUG: nspec was nan insted of 0
# - replaced wspec:: %PARAM-prefix with LADCP_wspec::
# Mar 12, 2017: - removed ANTSBIN (which is not public)
+# Dec 9, 2017: - added $antsSuppressCommonOptions = 1;
+# Dec 14, 2017: - added w.rms to output
+# May 13, 2018: - BUG: Removal of higher order polynomials (-o > 0) did not work
+# May 16, 2018: - modified depth.{min,max} to respect input resolution
+# Nov 1, 2018: - cosmetics
($ANTSLIB) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
($WCALC) = ($0 =~ m{^(.*)/[^/]*$});
@@ -52,13 +57,14 @@
# Usage
#----------------------------------------------------------------------
+$antsSuppressCommonOptions = 1;
&antsUsage('bc:dg:o:s:tuw:',0,
'[poly-o)rder <n[0]> to de-mean data; -1 to disable>] [suppress cosine-t)aper]',
'[-d)own/-u)pcast-only] [exclude -b)ottom window]',
'[shortwave -c)utoff <kz or lambda>]',
'[-s)urface <layer depth to exclude[150m]>',
'[-g)ap <max depth layer to fill with interpolation[40m]>]',
- '[-w)indow <power-of-two input-records[32]>]',
+ '[-w)indow <power-of-two input-records>]',
'[LADCP-profile(s)]');
&antsIntOpt(\$opt_o,0); # polynomial order to remove
@@ -132,9 +138,9 @@
$opt_g = int(($opt_g - 1) / $dz); # [m] -> [records]
-unless (defined($opt_w)) { # window size
+unless (defined($opt_w)) { # default window size: largest pwr-of-two <= 600m
for ($opt_w=32; $opt_w*$dz>600; $opt_w/=2) {}
- &antsInfo("%d-m windows ($opt_w records)",$opt_w*$dz);
+ &antsInfo("%d-m windows ($opt_w samples)",$opt_w*$dz);
}
&antsAddParams('LADCP_wspec::window_size',$opt_w,'LADCP_wspec::output_depth_resolution',$dz*$opt_w);
@@ -146,7 +152,7 @@
$resolution_bandwidth *= 2*$PI;
&antsAddParams('LADCP_wspec::resolution_bandwidth',$resolution_bandwidth);
-push(@antsNewLayout,'widx','depth','depth.min','depth.max','hab','nspec','w.nsamp.avg');
+push(@antsNewLayout,'widx','depth','depth.min','depth.max','hab','w.rms','nspec','w.nsamp.avg');
for (my($i)=0; $i<$opt_w/2+1; $i++) {
my($kz) = 2*$PI*$i/$zrange;
last if ($kz > $kzlim);
@@ -233,17 +239,34 @@
$opt_b = 1;
}
+ #--------------------------------------------------
+ # calculate rms w in window
+ # - also determines if there are missing y values
+ #--------------------------------------------------
+
+ my($dc_gap) = $opt_u; # exclude dc with -d, uc with -u
+ my($uc_gap) = $opt_d;
+ my($sumsq,$n) = (0,0);
+ for (my($r)=$fromR; $r<$fromR+$opt_w; $r++) {
+ if (numberp($ants_[$r][$dcwfnr])) {
+ $sumsq += $ants_[$r][$dcwfnr]**2;
+ $n++;
+ } else {
+ $dc_gap = 1;
+ }
+ if (numberp($ants_[$r][$ucwfnr])) {
+ $sumsq += $ants_[$r][$ucwfnr]**2;
+ $n++;
+ } else {
+ $uc_gap = 1;
+ }
+ }
+ my($wrms) = ($n > 0) ? sqrt($sumsq/$n) : nan;
+
#-----------------------------------
# output nan on non-numeric y values
#-----------------------------------
- my($dc_gap) = $opt_u; # exclude dc with -d, uc with -u
- my($uc_gap) = $opt_d;
- for (my($r)=$fromR; $r<$fromR+$opt_w; $r++) {
- $dc_gap |= !numberp($ants_[$r][$dcwfnr]);
- $uc_gap |= !numberp($ants_[$r][$ucwfnr]);
- }
-
if ($dc_gap && $uc_gap) {
push(@out,$ants_[$fromR+$opt_w/2][$zfnr]);
if ($ants_[0][$zfnr] > $ants_[1][$zfnr]) {
@@ -255,6 +278,7 @@
}
push(@out,defined($habfnr) ? # hab
avgF($habfnr,$fromR) : nan);
+ push(@out,$wrms); # rms w
push(@out,0); # nspec
push(@out,nan); # w.nsamp.avg
for ($i=0; $i<=$opt_w/2; $i++) { # power
@@ -285,7 +309,7 @@
&lfit($zfnr,$dcwfnr,-1,\@A,\@iA,\@covar,\&modelEvaluate,$fromR,$fromR+$opt_w-1);
&modelCleanup();
for (my($i)=0; $i<$opt_w; $i++) {
- modelEvaluate($i,$zfnr,\@afunc);
+ modelEvaluate($fromR+$i,$zfnr,\@afunc);
for ($calc=0,my($p)=1; $p<=$modelNFit; $p++) {
$calc += $A[$p] * $afunc[$p];
}
@@ -298,7 +322,7 @@
&lfit($zfnr,$ucwfnr,-1,\@A,\@iA,\@covar,\&modelEvaluate,$fromR,$fromR+$opt_w-1);
&modelCleanup();
for (my($i)=0; $i<$opt_w; $i++) {
- modelEvaluate($i,$zfnr,\@afunc);
+ modelEvaluate($fromR+$i,$zfnr,\@afunc);
for ($calc=0,my($p)=1; $p<=$modelNFit; $p++) {
$calc += $A[$p] * $afunc[$p];
}
@@ -341,16 +365,18 @@
@uc_pwr = &pgram_onesided($opt_w,@uc_coeff)
unless ($uc_gap);
- push(@out,$ants_[$fromR+$opt_w/2][$zfnr]); # middle z
+# push(@out,$ants_[$fromR+$opt_w/2][$zfnr]); # middle z
+ push(@out,avgF($zfnr,$fromR)); # average z
if ($ants_[0][$zfnr] > $ants_[1][$zfnr]) { # input descending
- push(@out,$ants_[$fromR+$opt_w-1][$zfnr]); # z.min
- push(@out,$ants_[$fromR][$zfnr]); # z.max
+ push(@out,$ants_[$fromR+$opt_w-1][$zfnr]-$dz/2); # z.min
+ push(@out,$ants_[$fromR][$zfnr]+$dz/2); # z.max
} else { # input ascending
- push(@out,$ants_[$fromR][$zfnr]); # z.min
- push(@out,$ants_[$fromR+$opt_w-1][$zfnr]); # z.max
+ push(@out,$ants_[$fromR][$zfnr]-$dz/2); # z.min
+ push(@out,$ants_[$fromR+$opt_w-1][$zfnr]+$dz/2); # z.max
}
push(@out,defined($habfnr) ? # hab
avgF($habfnr,$fromR) : nan);
+ push(@out,$wrms); # w.rms
my($nspec) = !$dc_gap + !$uc_gap; # nspec
push(@out,$nspec);
my($nsamp_sum) = my($nsn) = 0; # w.nsamp.avg
--- a/Makefile Thu Mar 16 11:53:27 2017 -0400
+++ b/Makefile Tue Nov 27 16:59:05 2018 -0500
@@ -1,13 +1,25 @@
#======================================================================
# M A K E F I L E
# doc: Mon Oct 17 13:29:27 2011
-# dlm: Tue Apr 21 20:55:33 2015
+# dlm: Wed Oct 31 10:19:00 2018
# (c) 2011 A.M. Thurnherr
-# uE-Info: 19 0 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 22 33 NIL 0 0 72 0 2 4 NIL ofnI
+#======================================================================
+
+# GO_SHIP archive target
+
+PROGS = LADCP_w_CTD LADCP_w_ocean LADCP_w_postproc LADCP_wspec LADCP_VKE
+LIBS = *.pl
+ANTSLIB = ANTSlib/.[ln]* ANTSlib/*
+A_TOOLS = ADCP_tools/RDI*pl
+
+LADCP_w_Software.tgz: ${PROGS} ${LIBS} ${ANTSLIB} ${A_TOOLS}
+ tar cvfz $@ $^
+
#======================================================================
MAKE_DIR = /Data/Makefiles
-include ${MAKE_DIR}/Makefile.GMT
+include ${MAKE_DIR}/Makefile.GMT5
w.cpt:
mkCPT -oc polar -- -0.07 -0.05 -0.04 -0.03 -0.02 -0.01 0.01 0.02 0.03 0.04 0.05 0.07 > $@
--- a/bottom_tracking.pl Thu Mar 16 11:53:27 2017 -0400
+++ b/bottom_tracking.pl Tue Nov 27 16:59:05 2018 -0500
@@ -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 May 24 16:34:43 2016
+# dlm: Tue May 1 21:48:54 2018
# (c) 2010 A.M. Thurnherr
-# uE-Info: 116 5 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 21 14 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -15,8 +15,10 @@
# Mar 4, 2014: - removed old unused code
# Jan 26, 2016: - added %PARAMs
# May 24, 2016: - calc_binDepths() -> binDepths()
+# May 1, 2018: - log-file cosmetics
-# This code is essentially identical to the one used in LADCPproc. Differences:
+# This code is derived from the one used in LADCPproc, with the following
+# differences:
# 1) velocity editing is simpler: no wake editing, no PPI editing, no shear
# editing, no w outlier
# 2) median/mad calculated instead of mean/stddev
@@ -129,7 +131,7 @@
progress("\t$nBTfound BT ensembles found\n");
progress("\t$nBTdepthFlag flagged bad because of wrong bottom depth\n");
progress("\t$nBTvalidVelFlag flagged bad because of no valid velocities\n");
- progress("\t$nBTwFlag flagged bad because of incorrect vertical velocity\n");
+ progress("\t$nBTwFlag flagged bad because of inconsistent vertical velocity\n");
for (my($gi)=0; $gi<@BTw; $gi++) { # calc grid medians & mads
$BT{N_SAMP}[$gi] = @{$BTw[$gi]};
--- a/defaults.pl Thu Mar 16 11:53:27 2017 -0400
+++ b/defaults.pl Tue Nov 27 16:59:05 2018 -0500
@@ -1,9 +1,9 @@
#======================================================================
# D E F A U L T S . P L
# doc: Tue Oct 11 17:11:21 2011
-# dlm: Sun Mar 12 12:53:44 2017
+# dlm: Wed May 2 14:11:50 2018
# (c) 2011 A.M. Thurnherr
-# uE-Info: 459 36 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 308 44 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -82,6 +82,11 @@
# this file is read before the options are parsed
# Aug 5, 2016: - updated header
# Dec 22, 2016: - added $opt_p
+# Nov 27, 2017: - added @valid_ensemble_range
+# Nov 29, 2017: - replaced opt_i by initial_time_lag
+# Apr 24, 2018: - added $water_depth_db_cmd
+# May 2, 2018: - added max_hspeed
+# - replaced $PPI_seabed_editing_required by &PPI_seabed_editing_required
#======================================================================
# Output Log Files
@@ -113,6 +118,12 @@
# Input Data
#======================================================================
+# The two elements in the @valid_ensemble_range array limit the minimum
+# and maximum ensemble numbers considered during processing. This is
+# useful primarily for files with lots of on-deck data.
+
+# @valid_ensemble_range = (3000,10000)
+
# Set $opt_4 to 1 (or use the -4 option) to suppress 3-beam LADCP
# solutions. This only has an effect for beam-coordinate data.
@@ -247,12 +258,32 @@
$surface_layer_depth = 25;
+# Water depth is important for precious ping interference editing
+# (see below) and for setting the height-above bottom field.
+# - by default, water depth for dowwnward-facing ADCPs is determined
+# from the seabed echo return
+# - when water depth is set explicitly either via the $water_depth or
+# the $opt_h variable no search for the seabed is done
+# - the variable $water_depth_db_cmd can be set to the name of an
+# external command, which is used to get nominal water depth
+# from a data base if there is no other water depth information
+# (from echo return or supplied by user). The command will be
+# called with longitude and latitude as the only arguments and
+# is expected to return the water depth in meters.
+
+#$water_depth = 2048; # uncomment to set water depth to 2048m
+#$opt_h = 2048; # uncomment to set water depth to 2048m
+#$water_depth_db_cmd = 'waterdepth'; # uncomment to use 'waterdepth' command to get water depth
+
+
# Previous Ping Interference editing as described in [edit_data.pl]
-# - enabled by default for WH150 data
-# - 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;
+# - enabled by default seabed editing of WH150 data but nothing else
+# - PPI_seabed_editing_required is a function that is called
+# once the data are loaded for the downlooker only; if it
+# returns true, seabed PPI editing is enabled
+# - PPI_surface_editing_required is a function that is called
+# once the data are loaded for the uplooker only; if it
+# returns true, sea surface PPI editing is enabled
# - 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
@@ -263,10 +294,18 @@
# set by the shortest acoustic path between the ADCP and the
# seabed.
-$PPI_seabed_editing_required = '($LADCP{BEAM_FREQUENCY} < 300)';
+sub PPI_seabed_editing_required()
+{
+# return 1; # uncomment to enable unconditional PPI editing
+ return ($LADCP{BEAM_FREQUENCY} < 300); # low-frequency instruments require PPI editing
+}
-#$PPI_editing = 1; # uncomment to enable PPI always
-#$PPI_extend_upper_limit = 1.03; # see comments above
+sub PPI_surface_editing_required()
+{
+ return 0; # no sea surface PPI editing by default
+}
+
+#$PPI_extend_upper_limit = 1.03; # see comments above
# The following variables control the "non-obvious" sidelobe editing for
@@ -285,14 +324,35 @@
$vessel_draft = 6; # in meters
+# The following function, which is called after the LADCP data have been
+# read, must return the maximum horizontal reference-layer
+# speed that is allowed. The following values are based on 2018 GO-SHIP
+# S4P profile #106 where the CTD rosette was dragged quickly during
+# the latter part of the upcast. Of course, it is possible that the
+# differnces between the UL and DL data could be due to tilt-
+# sensor differences, rather than due to instrument type.
+
+sub max_hspeed()
+{
+ if (abs($LADCP{BEAM_FREQUENCY}-300) <= 25) { # 300kHz Workhorse
+ $max_hspeed = 0.55; # meters/second
+ } elsif (abs($LADCP{BEAM_FREQUENCY}-150) <= 25) { # 150kHz Workhorse
+ $max_hspeed = 0.35; # meters/second
+ } else {
+ warning(2,"unknown horizontal speed limit for this instrument frequency ($LADCP{BEAM_FREQUENCY} kHz)\n");
+ $max_hspeed = 9e99;
+ }
+}
+
#======================================================================
# Time Lagging
#======================================================================
-# The -i option allows defining an initial guess for the time lag between
-# the LADCP and the CTD data.
+# The following variable allows specifying an initial guess for the time
+# lag between the LADCP and the CTD data. The -i option overrides any value
+# set in the [ProcessingParams] file.
-# $opt_i = 567;
+# $initial_time_lag = 567;
# The following variables define the bins used to calculate the reference-
--- a/edit_data.pl Thu Mar 16 11:53:27 2017 -0400
+++ b/edit_data.pl Tue Nov 27 16:59:05 2018 -0500
@@ -1,9 +1,9 @@
#======================================================================
-# E D I T _ D A T A . P L
+# / D A T A / S R C / O C E A N O G R A P H Y / L A D C P _ V E R T I C A L _ V E L O C I T Y / E D I T _ D A T A . P L
# doc: Sat May 22 21:35:55 2010
-# dlm: Mon Jun 6 21:13:28 2016
+# dlm: Tue Nov 27 11:07:33 2018
# (c) 2010 A.M. Thurnherr
-# uE-Info: 543 0 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 46 71 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -41,6 +41,9 @@
# - verified that removed velocities are counted correctly
# Jun 3, 2016: - BUG: applyTiltCorrection() did not use gimbal_pitch
# Jun 6, 2016: - removed applyTiltCorrection()
+# Oct 13, 2017: - BUG: editPPI() only allowed for nominal transducer frequencies
+# May 1, 2018: - added editLargeHSpeeds()
+# Nov 17, 2018: - BUG: spurious letter "z" had crept in at some stage
# NOTES:
# - all bins must be edited (not just the ones between $LADCP_firstBin
@@ -53,6 +56,7 @@
#======================================================================
# correctAttitude($ens,$pitch_bias,$roll_bias,$heading_bias)
+# - attitude bias correction
# - this is called before gimbal_pitch is calculated
#======================================================================
@@ -371,11 +375,11 @@
my($nerm) = 0; # of ensembles affected
unless (defined($bha)) {
- if ($LADCP{BEAM_FREQUENCY} == 1200) { $bha = 2.4; }
- elsif ($LADCP{BEAM_FREQUENCY} == 600) { $bha = 2.5; }
- elsif ($LADCP{BEAM_FREQUENCY} == 300) { $bha = 3.7; }
- elsif ($LADCP{BEAM_FREQUENCY} == 150) { $bha = 6.7; }
- elsif ($LADCP{BEAM_FREQUENCY} == 75) { $bha = 8.4; }
+ if (abs($LADCP{BEAM_FREQUENCY}-1200)/1200 <= 0.1) { $bha = 2.4; }
+ elsif (abs($LADCP{BEAM_FREQUENCY}-600) / 600 <= 0.1) { $bha = 2.5; }
+ elsif (abs($LADCP{BEAM_FREQUENCY}-300) / 300 <= 0.1) { $bha = 3.7; }
+ elsif (abs($LADCP{BEAM_FREQUENCY}-150) / 150 <= 0.1) { $bha = 6.7; }
+ elsif (abs($LADCP{BEAM_FREQUENCY}-75) / 75 <= 0.1) { $bha = 8.4; }
else { error("$0: unexpected transducer frequency $LADCP{BEAM_FREQUENCY}\n"); }
}
@@ -542,5 +546,26 @@
}
#======================================================================
+# $nerm = editLargeHSpeeds($fe,$te,$max_hspeed)
+# - filter based on 2018 GO-SHIP LADCP profile #106, where UL
+# velocities become bad when ship starts dragging rosette
+# - only valid bin range is edited/counted
+#======================================================================
+
+sub editLargeHSpeeds($$$)
+{
+ my($fe,$te,$limit) = @_;
+ my($nerm) = 0;
+ for (my($ens)=$fe; $ens<=$te; $ens++) {
+ next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
+ next unless (vel_speed($LADCP{ENSEMBLE}[$ens]->{REFLR_U},
+ $LADCP{ENSEMBLE}[$ens]->{REFLR_U}) > $limit);
+ undef($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
+ $nerm++;
+ }
+ return $nerm;
+}
+
+#======================================================================
1;
--- a/plot_time_lags.pl Thu Mar 16 11:53:27 2017 -0400
+++ b/plot_time_lags.pl Tue Nov 27 16:59:05 2018 -0500
@@ -1,9 +1,9 @@
#======================================================================
# P L O T _ T I M E _ L A G S . P L
# doc: Tue Jul 28 13:21:09 2015
-# dlm: Tue Mar 7 12:04:21 2017
+# dlm: Thu Mar 22 10:56:17 2018
# (c) 2015 A.M. Thurnherr
-# uE-Info: 35 0 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 40 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -12,7 +12,8 @@
# Mar 16, 2016: - adapted to gmt5
# May 18, 2016: - added version
# May 24, 2016: - fixed for partial-depth casts
-# mar 7, 2017: - added time lines for -p
+# Mar 7, 2017: - added time lines for -p
+# Mar 22, 2018: - removed plotting of yellow runs on -l
require "$ANTS/libGMT.pl";
@@ -35,11 +36,13 @@
printf(GMT "%f $ymin\n%f $ymax\n>\n",$x,$x);
}
- GMT_psxy('-W8,yellow'); # offset runs
- for (my($i)=0; $i<@bmo_buf; $i++) {
- printf(GMT ">\n%f %f\n%f %f\n",
- $fg_buf[$i]/60-0.5,$bmo_buf[$i],
- $lg_buf[$i]/60+0.5,$bmo_buf[$i]);
+ unless ($opt_l) {
+ GMT_psxy('-W8,yellow'); # indicate valid runs
+ for (my($i)=0; $i<@bmo_buf; $i++) {
+ printf(GMT ">\n%f %f\n%f %f\n",
+ $fg_buf[$i]/60-0.5,$bmo_buf[$i],
+ $lg_buf[$i]/60+0.5,$bmo_buf[$i]);
+ }
}
GMT_psxy('-Sc0.1 -Gcoral'); # individual offsets
--- a/plot_wprof.pl Thu Mar 16 11:53:27 2017 -0400
+++ b/plot_wprof.pl Tue Nov 27 16:59:05 2018 -0500
@@ -1,9 +1,9 @@
#======================================================================
# P L O T _ W P R O F . P L
# doc: Sun Jul 26 11:08:50 2015
-# dlm: Thu May 26 11:29:32 2016
+# dlm: Tue Mar 20 15:26:21 2018
# (c) 2015 A.M. Thurnherr
-# uE-Info: 20 0 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 22 53 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -18,6 +18,8 @@
# - fixed for partial-depth profiles
# - suppress plotting of nsamp == 0
# May 26, 2016: - added instrument coord system to plot labels
+# Mar 20, 2018: - BUG: units of vertical package acceleration were wrong
+# - added blue background for likely in-ice package accelerations
# Tweakables:
#
@@ -164,19 +166,23 @@
}
printf(GMT "0.91 1.090 %.1f\\260\n",$P{uc_mean_tilt});
- if ($P{dc_rms_accel_pkg} < 0.7) {
+ if ($P{dc_rms_accel_pkg} < 0.1) {
+ GMT_pstext('-F+f9,Helvetica,coral+jTL -Gblue -N');
+ } elsif ($P{dc_rms_accel_pkg} < 0.7) {
GMT_pstext('-F+f9,Helvetica,coral+jTL -N');
} else {
GMT_pstext('-F+f9,Helvetica,coral+jTL -Gyellow -N');
}
- printf(GMT "0.78 1.125 %.1fm\@+2\@+/s\n",$P{dc_rms_accel_pkg});
+ printf(GMT "0.78 1.125 %.1fm/s\@+2\@+\n",$P{dc_rms_accel_pkg});
- if ($P{uc_rms_accel_pkg} < 0.7) {
+ if ($P{uc_rms_accel_pkg} < 0.1) {
+ GMT_pstext('-F+f9,Helvetica,SeaGreen+jTL -Gblue -N');
+ } elsif ($P{uc_rms_accel_pkg} < 0.7) {
GMT_pstext('-F+f9,Helvetica,SeaGreen+jTL -N');
} else {
GMT_pstext('-F+f9,Helvetica,SeaGreen+jTL -Gyellow -N');
}
- printf(GMT "0.89 1.125 %.1fm\@+2\@+/s\n",$P{uc_rms_accel_pkg});
+ printf(GMT "0.89 1.125 %.1fm/s\@+2\@+\n",$P{uc_rms_accel_pkg});
my($depth_tics) = ($plot_wprof_ymax-$plot_prof_ymin < 1000 ) ? 'f10a100' : 'f100a500'; # AXES
setR1();
--- a/time_lag.pl Thu Mar 16 11:53:27 2017 -0400
+++ b/time_lag.pl Tue Nov 27 16:59:05 2018 -0500
@@ -1,9 +1,9 @@
#======================================================================
# T I M E _ L A G . P L
# doc: Fri Dec 17 21:59:07 2010
-# dlm: Tue Mar 7 09:03:20 2017
+# dlm: Tue Oct 16 16:26:06 2018
# (c) 2010 A.M. Thurnherr
-# uE-Info: 73 63 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 80 38 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -71,6 +71,13 @@
# Mar 7, 2016: - BUG: editing did not work correctly in all cases
# Mar 6, 2017: - BUG: assertion in mad_w failed with 2017 P18 DL#206
# Mar 9, 2017: - tightened timelag editing (good_diff: 2->1)
+# Mar 22, 2018: - re-wrote heuristics to remove lags with large mads
+# - BUG: bestLag with 1 valid sample returned 0 mad
+# - BUG: timelag editing did not work correctly when there was not a sufficiently long valid lag
+# Mar 27, 2018: - BUG: re-written heuristic could fail when there was a valid but unpopular lag with very low mad.
+# Solution: remove very unpopular lags first
+# Oct 4, 2018: - added timelagging debug code
+# Oct 16, 2018: - removed debug code
# DIFFICULT STATIONS:
# NBP0901#131 this requires the search-radius doubling heuristic
@@ -102,22 +109,21 @@
$CTD_mean_w += $CTD{W}[$s+$so];
$nsamp++;
}
- return 9e99 unless ($nsamp);
+ return 9e99 unless ($nsamp > 1);
$LADCP_mean_w /= $nsamp;
$CTD_mean_w /= $nsamp;
- my($sad) = 0; # now, calculate mad
+ my($sad) = $nsamp = 0; # now, calculate mad
for (my($e)=$fe; $e<=$le; $e++) {
- next unless numberp($LADCP{ENSEMBLE}[$e]->{REFLR_W});
my($s) = int(($LADCP{ENSEMBLE}[$e]->{ELAPSED} + $CTD{TIME_LAG} - $CTD{ELAPSED}[0]) / $CTD{DT} + 0.5);
next unless ($s>=0 && $s<=$#{$CTD{ELAPSED}});
next unless numberp($LADCP{ENSEMBLE}[$e]->{REFLR_W});
my($dw) = $LADCP{ENSEMBLE}[$e]->{REFLR_W}-$LADCP_mean_w - ($CTD{W}[$s+$so]-$CTD_mean_w);
# print(STDERR "dw = $dw ($LADCP{ENSEMBLE}[$e]->{REFLR_W}-$LADCP_mean_w - ($CTD{W}[$s+$so]-$CTD_mean_w)\n");
next unless (abs($dw) <= $w_max_lim);
- $sad += abs($dw);
+ $sad += abs($dw); $nsamp++;
}
- return $sad/$nsamp;
+ return $nsamp ? $sad/$nsamp : 9e99;
}
@@ -198,18 +204,34 @@
$n_valid_windows++;
$nBest{$so}++; $madBest{$so} += $mad;
}
+ my($maxN) = 0;
foreach my $i (keys(%nBest)) {
+ $maxN = $nBest{$i} if ($nBest{$i} > $maxN);
$madBest{$i} /= $nBest{$i};
}
- my($med_mad) = median(values(%madBest)); # remove lags with large mads
- my($mad_mad) = mad2($med_mad,values(%madBest));
+ foreach my $lag (keys(%nBest)) { # remove unpopular lags
+ next if ($nBest{$lag} >= $maxN/10);
+ $n_valid_windows -= $nBest{$lag};
+ $nBest{$lag} = 0; $madBest{$lag} = 9e99;
+ }
+
+## my($med_mad) = median(values(%madBest)); # remove lags with large mads
+## my($mad_mad) = mad2($med_mad,values(%madBest));
+## foreach my $lag (keys(%nBest)) {
+## next if ($madBest{$lag} <= $med_mad+$mad_mad);
+## $n_valid_windows -= $nBest{$lag};
+## $nBest{$lag} = 0;
+## }
+
+ my($min_mad) = min(values(%madBest)); # remove lags with large mads
foreach my $lag (keys(%nBest)) {
- next if ($madBest{$lag} <= $med_mad+$mad_mad);
+ next if ($madBest{$lag} <= 3*$min_mad);
$n_valid_windows -= $nBest{$lag};
$nBest{$lag} = 0;
}
+
my(@best_lag); # find 3 most popular lags
foreach my $lag (keys(%nBest)) {
$best_lag[0] = $lag if ($nBest{$lag} > $nBest{$best_lag[0]});
@@ -236,7 +258,6 @@
$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]} # require quorum
>= $opt_3*$n_valid_windows) {
if (max(@best_lag)-min(@best_lag) > $TL_max_allowed_three_lag_spread) {
@@ -278,6 +299,16 @@
#----------------------------------------------------
# Here, either $failed is set, or we have a valid lag.
+ #----------------------------------------------------
+
+# if ($failed) {
+# for (my($wi)=0; $wi<$n_windows; $wi++) {
+# print(STDERR "$wi $so[$wi] $mad[$wi]\n");
+# }
+# }
+
+ #----------------------------------------------------
+ # 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
@@ -354,6 +385,7 @@
($lg[$i] == $n_windows-1) ? $LADCP{ENSEMBLE}[$lastGoodEns]->{ELAPSED} : $elapsed[$lg[$i]]);
}
# print(STDERR "elim = @elim\n");
+ $failed = 1 unless (@elim);
$nerm = $failed
? editBadTimeLagging($first_ens,$last_ens,-1)
: editBadTimeLagging($first_ens,$last_ens,@elim);
--- a/time_series.pl Thu Mar 16 11:53:27 2017 -0400
+++ b/time_series.pl Tue Nov 27 16:59:05 2018 -0500
@@ -1,9 +1,9 @@
#======================================================================
# T I M E _ S E R I E S . P L
# doc: Sun May 23 16:40:53 2010
-# dlm: Wed Apr 17 17:05:16 2013
+# dlm: Wed May 2 11:23:48 2018
# (c) 2010 A.M. Thurnherr
-# uE-Info: 20 63 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 24 57 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -18,12 +18,20 @@
# Oct 12, 2011: - re-worked ref_lr_w()
# - stopped depth integration across gaps >= 5s
# Apr 17, 2013: - improved gap message (added ensemble range)
+# Nov 27, 2017: - BUG: gap heuristic could not deal with P06#001
+# - BUG: gap heuristic could not deal with P06#025
+# May 1, 2018: - added reflr u and v calculations
+# - BUG: reflr u and v calcs did not work
# NOTES:
# - resulting DEPTH field based on integrated w without any sound speed correction
# - single-ping ensembles assumed, i.e. no percent-good tests applied
# - specified bin numbers are 1-relative
+#----------------------------------------------------------------------
+# Reference-Layer Velocities
+#----------------------------------------------------------------------
+
sub ref_lr_w($$$$) # calc ref-layer vert vels
{
my($dta,$ens,$rl_b0,$rl_b1) = @_;
@@ -39,6 +47,22 @@
$dta->{ENSEMBLE}[$ens]->{REFLR_W_NSAMP} = @w;
}
+sub ref_lr_uv($$$$) # calc ref-layer horiz vels
+{
+ my($dta,$ens,$rl_b0,$rl_b1) = @_;
+ my(@u,@v);
+
+ for (my($bin)=$rl_b0-1; $bin<=$rl_b1-1; $bin++) {
+ next unless defined($dta->{ENSEMBLE}[$ens]->{U}[$bin]);
+ die unless numbersp($dta->{ENSEMBLE}[$ens]->{U}[$bin],$dta->{ENSEMBLE}[$ens]->{V}[$bin]);
+ push(@u,$dta->{ENSEMBLE}[$ens]->{U}[$bin]);
+ push(@v,$dta->{ENSEMBLE}[$ens]->{V}[$bin]);
+ }
+ return unless (@u);
+ $dta->{ENSEMBLE}[$ens]->{REFLR_U} = avg(@u); $dta->{ENSEMBLE}[$ens]->{REFLR_V} = avg(@v);
+ $dta->{ENSEMBLE}[$ens]->{REFLR_UV_NSAMP} = @u;
+}
+
#======================================================================
# ($firstgood,$lastgood,$atbottom,$w_gap_time) =
# calcLADCPts($dta,$skip_ens,$lr_b0,$lr_b1,$min_corr,$max_e,$max_gap);
@@ -76,19 +100,22 @@
$dta->{ENSEMBLE}[$lastgood]->{UNIX_TIME}; # ... last good ens
if ($dt > $max_gap) {
- if ($max_depth>50 && $depth<0.1*$max_depth) {
- warning(1,"long gap (%ds) near end of profile --- terminated at ensemble #$dta->{ENSEMBLE}[$e]->{NUMBER}\n",$dt);
- last;
- }
- if ($depth < 10) {
- warning(1,"long gap (%ds) near beginning of profile --- restarted at ensemble #$dta->{ENSEMBLE}[$e]->{NUMBER}\n",$dt);
- $firstgood = $lastgood = $e;
- undef($atbottom); undef($max_depth);
- $depth = 0;
- $dta->{ENSEMBLE}[$e]->{ELAPSED} = 0;
- $dta->{ENSEMBLE}[$e]->{DEPTH} = $depth;
- $w_gap_time = 0;
- next;
+ if (($max_depth>50 && abs($depth)<0.1*$max_depth) && # looks like a profile
+ (@{$dta->{ENSEMBLE}}-$e < 0.25*@{$dta->{ENSEMBLE}})) { # in the final quartile of the data
+ warning(1,"long gap (%ds) after likely profile (0->%d->%dm) --- finishing at ens#$dta->{ENSEMBLE}[$e]->{NUMBER}\n",
+ $dt,$max_depth,$depth);
+ last;
+ } elsif ((abs($depth) < 10) || # shallow gap at the beginning
+ ($depth == $max_depth)) { # biased in-air data
+ warning(1,"long surface gap (%ds) --- restarting at ens#$dta->{ENSEMBLE}[$e]->{NUMBER}\n",$dt);
+ warning(1,"[depth = $depth, max_depth = $max_depth]\n");
+ $firstgood = $lastgood = $e;
+ undef($atbottom); undef($max_depth);
+ $depth = 0;
+ $dta->{ENSEMBLE}[$e]->{ELAPSED} = 0;
+ $dta->{ENSEMBLE}[$e]->{DEPTH} = $depth;
+ $w_gap_time = 0;
+ next;
}
if ($dta->{ENSEMBLE}[$e]->{ELAPSED} < 200) {
warning(1,"long gap (%ds) at ensembles #$dta->{ENSEMBLE}[$lastgood]->{NUMBER}-$dta->{ENSEMBLE}[$e]->{NUMBER}, %ds into the profile\n",
@@ -107,7 +134,13 @@
$lastgood = $e;
}
+ for (my($e)=$firstgood; $e<=$lastgood; $e++) { # calculate u and v
+ ref_lr_uv($dta,$e,$rl_b0,$rl_b1);
+ }
+
return ($firstgood,$lastgood,$atbottom,$w_gap_time);
}
+#----------------------------------------------------------------------
+
1;
--- a/version.pl Thu Mar 16 11:53:27 2017 -0400
+++ b/version.pl Tue Nov 27 16:59:05 2018 -0500
@@ -1,9 +1,9 @@
#======================================================================
# V E R S I O N . P L
# doc: Tue Oct 13 10:40:57 2015
-# dlm: Sun Mar 12 12:11:05 2017
+# dlm: Tue Nov 27 13:15:48 2018
# (c) 2015 A.M. Thurnherr
-# uE-Info: 26 15 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 35 29 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -18,13 +18,22 @@
# May 12, 2016: - V1.2
# May 19, 2016: - updated ADCP tools to V1.6
# Aug 5, 2017: - updated ANTS lib to V6.7
-# Mar 12, 2017: - updated ANTS lit to V6.8
+# Mar 12, 2017: - updated ANTS lib to V6.8
+# Mar 15, 2017: - V1.3
+# Nov 28, 2017: - V1.4 (perl-tools 2.0; ANTSlib 6.9)
+# Sep 13, 2018: - added '.' to library path to allow do without full pathname
+# - added 1; to the end
+# Nov 27, 2018: - updated ANTS lib to V7.1
+# - updated ADCP tools to V2.2
#$VERSION = '1.1'; # Jan 4, 2016
#$VERSION = '1.2'; # May 12, 2016
-
-$VERSION = '1.3';
+#$VERSION = '1.3'; # Mar 15, 2017
+$VERSION = '1.4'; # Nov 28, 2017
-$antsMinLibVersion = 6.8; # Feb 2017: - .lsfit.poly, .nminterp.linear added
-$ADCP_tools_minVersion = 1.9; # Feb 2017: - round() namespace clash resolved
+$antsMinLibVersion = 7.1;
+$ADCP_tools_minVersion = 2.2;
+use lib '.';
+
+1;