Version 1.2 finished. Out for testing to Jay Hooper.
--- a/HISTORY Tue Feb 02 14:55:24 2016 +0000
+++ b/HISTORY Tue Mar 08 14:09:57 2016 +0000
@@ -1,9 +1,9 @@
======================================================================
H I S T O R Y
doc: Mon Oct 12 16:09:24 2015
- dlm: Sun Jan 24 13:47:15 2016
+ dlm: Tue Mar 8 13:57:23 2016
(c) 2015 A.M. Thurnherr
- uE-Info: 39 40 NIL 0 0 72 3 2 4 NIL ofnI
+ uE-Info: 41 47 NIL 0 0 72 3 2 4 NIL ofnI
======================================================================
V1.0:
@@ -37,3 +37,5 @@
- added QC for mean residuals
- added QC for dual-head wprofs
- added automatic editing of bad time lagged data
+ Jan 25 -- Mar 8:
+ - many bug fixes and small improvements
--- a/LADCP_VKE Tue Feb 02 14:55:24 2016 +0000
+++ b/LADCP_VKE Tue Mar 08 14:09:57 2016 +0000
@@ -2,9 +2,9 @@
#======================================================================
# L A D C P _ V K E
# doc: Tue Oct 14 11:05:16 2014
-# dlm: Mon Jan 25 09:31:12 2016
+# dlm: Tue Mar 8 13:53:35 2016
# (c) 2012 A.M. Thurnherr
-# uE-Info: 44 49 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 195 24 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
$antsSummary = 'calculate VKE from LADCP-derived vertical-velocity profiles';
@@ -42,6 +42,12 @@
# Dec 27, 2015: - reduced minlim for eps to 1e-13 W/kg
# Dec 29, 2015: - added 3rd consistency check (p0 limit)
# Jan 25, 2016: - added software version %PARAM
+# Mar 7, 2016: - removed unused spectral normalization code
+# - added latitude constraint for calculation of eps.w
+# Mar 8, 2016: - renamed fs_pwr to pwr.fs for consistency
+# - renamed eps.w to eps.VKE for consistency
+# - changed default output filenames for -d and -u
+# - removed ./ from figure label
($ANTSLIB) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
($WCALC) = ($0 =~ m{^(.*)/[^/]*$});
@@ -127,12 +133,12 @@
if (defined($opt_c)) { # shortwave cutoff supplied
$lmin = ($opt_c < 1) ? 2*$PI/$opt_c : $opt_c;
- &antsAddParams('LADCP_VKE::shortwave_cutoff',2*$PI/$lmin); # ensure eps.w is calculated below
+ &antsAddParams('LADCP_VKE::shortwave_cutoff',2*$PI/$lmin); # ensure eps.VKE is calculated below
} elsif (defined(antsParam('shortwave_cutoff'))) { # cutoff already applied
$lmin = 2*$PI/antsParam('shortwave_cutoff');
} else { # use 100m default cutoff
$lmin = 100;
- &antsAddParams('LADCP_VKE::shortwave_cutoff',2*$PI/$lmin); # ensure eps.w is calculated below
+ &antsAddParams('LADCP_VKE::shortwave_cutoff',2*$PI/$lmin); # ensure eps.VKE is calculated below
}
$lmax = 9e99; # no longwave cutoff implemented yet
@@ -175,9 +181,18 @@
croak("$opt_f: not a directory\n") unless (-d $opt_f);
my($id) = &antsRequireParam('profile_id');
- $opt_p = sprintf('%s/%03d_VKE.ps',$opt_f,$id)
- unless defined($opt_p);
- $outfile = sprintf('%s/%03d.VKE',$opt_f,$id);
+
+ if ($opt_d) {
+ $opt_p = sprintf('%s/%03d_dc_VKE.ps',$opt_f,$id) unless defined($opt_p);
+ $outfile = sprintf('%s/%03d_dc.VKE',$opt_f,$id);
+ } elsif ($opt_u) {
+ $opt_p = sprintf('%s/%03d_uc_VKE.ps',$opt_f,$id) unless defined($opt_p);
+ $outfile = sprintf('%s/%03d_uc.VKE',$opt_f,$id);
+ } else {
+ $opt_p = sprintf('%s/%03d_VKE.ps',$opt_f,$id) unless defined($opt_p);
+ $outfile = sprintf('%s/%03d.VKE',$opt_f,$id);
+ }
+ $outfile =~ s@^./@@;
open(STDOUT,">$outfile") || die("$outfile: $!\n");
} elsif (defined($opt_f)) {
croak("-f can only be used without STDOUT redirection\n");
@@ -199,16 +214,6 @@
}
-sub normalize_spectral_power($) # scale pwrdens by p0
-{
- my($r) = @_;
-
- for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) {
- $ants_[$r][$f] /= $ants_[$r][$p0f];
- }
-}
-
-
sub fit_universal_w_spec($) # vertical velocity => p0
{
my($r) = @_;
@@ -254,13 +259,13 @@
# Process File
#----------------------------------------------------------------------
-$fspwrf = &antsNewField('fs_pwr'); # derived fields
+$fspwrf = &antsNewField('pwr.fs'); # derived fields
$p0f = &antsNewField('p0'); # VKE parameterization
$rmsf = &antsNewField('p0fit.rms');
$sampf = &antsNewField('p0fit.nsamp');
$rf = &antsNewField('p0fit.r');
-$wepsf = &antsNewField('eps.w');
+$wepsf = &antsNewField('eps.VKE');
my(@outLayout) = @antsNewLayout; # save for later
for ($f=0; $f<@outLayout; $f++) { # determine last spectral field in input
@@ -288,14 +293,18 @@
my($min_depth) = 9e99;
my($max_depth) = -9e99;
-for (my($r)=0; $r<@ants_; $r++) { # loop over all windows
+my($latM) = abs(&antsRequireParam('lat'));
+&antsInfo("WARNING: low latitude-profile no epsilon estimated")
+ unless ($latM > 5);
+
+for (my($r)=0; $r<@ants_; $r++) { # loop over all windows
#--------------------------
# apply spectral correction
#--------------------------
unless ($opt_m) {
- for (my($i)=0; $i<$nfreq; $i++) { # loop over wavenumbers
+ for (my($i)=0; $i<$nfreq; $i++) { # loop over wavenumbers
$ants_[$r][$i+$pg_fmin] *=
T_w(antsParam("k.$i"),antsParam('ADCP_bin_length'),
antsParam('ADCP_pulse_length'),antsParam('input_depth_resolution'),
@@ -309,9 +318,9 @@
integrate_fs_power($r); # calculate total finescale power
fit_universal_w_spec($r); # fit kz^-2 spectrum
-# normalize_spectral_power($r); # normalize by p0 (even w/o shortwave cutoff)
- # consistency checks:
- if ($ants_[$r][$rmsf] <= 0.4 && # 1) rms from kz fit < 0.4 (as in paper)
+
+ if ($latM > 5 && # Tests: 0) not equatorial
+ $ants_[$r][$rmsf] <= 0.4 && # 1) rms from kz fit < 0.4 (as in paper)
$ants_[$r][$rf] < 0 && # 2) raw spectra are red spectrum
$ants_[$r][$p0f] > 1e-7) { # 3) exclude low-p0 tail apparent at eps<2e-10 W/kg
$ants_[$r][$wepsf] = ($ants_[$r][$p0f] / $opt_e)**2;
--- a/LADCP_w_CTD Tue Feb 02 14:55:24 2016 +0000
+++ b/LADCP_w_CTD Tue Mar 08 14:09:57 2016 +0000
@@ -2,9 +2,9 @@
#======================================================================
# L A D C P _ W _ C T D
# doc: Mon Nov 3 17:34:19 2014
-# dlm: Mon Jan 25 09:32:32 2016
+# dlm: Fri Feb 19 08:01:38 2016
# (c) 2014 A.M. Thurnherr
-# uE-Info: 51 0 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 241 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
$antsSummary = 'pre-process SBE 9plus CTD data for LADCP_w';
@@ -48,6 +48,10 @@
# 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
+# Feb 16, 2016: - BUG: %sampling_interval was erroenously called %sampling_frequency
+# Feb 19, 2016: - improved on-surface data detection (conductivity <= 10mS/cm)
+# - BUG: temperatures < 0 were not allowed
+# - temperature histogram resolution increased from 1degC to 0.5degC for no good reason
($ANTS) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
($WCALC) = ($0 =~ m{^(.*)/[^/]*$});
@@ -169,7 +173,7 @@
my($guess_bin) = @_;
my($min_bin,$max_bin);
- croak("assertion failed") unless ($hist[$guess_bin]);
+ die("assertion failed") unless ($hist[$guess_bin]);
for ($max_bin=$guess_bin; $hist[$max_bin]; $max_bin++) { }
for ($min_bin=$guess_bin; $hist[$min_bin]; $min_bin--) { }
return ($min_bin,$max_bin+1);
@@ -179,12 +183,18 @@
print(STDERR "Editing Data...") if ($opt_v);
#----------------------------------------
- # trim initial records with nan pressures
+ # trim initial records with
+ # - nan pressure
+ # - nan conductivity
+ # - conductivity <= 10 mS/cm
#----------------------------------------
my($trimmed) = 0; # trim leading nan pressures
shift(@ants_),$trimmed++
- until numberp($ants_[0][$pressF]) && numberp($ants_[0][$condF]);
- printf(STDERR "\n\t%d initial records with nan pressure and/or conductivity trimmed",$trimmed) if ($opt_v > 1);
+ until numberp($ants_[0][$pressF]) &&
+ numberp($ants_[0][$condF]) &&
+ (($P{cond.unit} eq 'mS/cm' && $ants_[0][$condF] > 10) ||
+ ($P{cond.unit} eq 'S/m' && $ants_[0][$condF] > 1));
+ printf(STDERR "\n\t%d initial at-surface records trimmed",$trimmed) if ($opt_v > 1);
my($lvp) = $ants_[0][$pressF];
my($lvc) = $ants_[0][$condF];
@@ -227,7 +237,7 @@
$modeSamp = $hist[$b]; $modeBin = $b;
}
($min,$max) = validRange($modeBin);
- $min /= $condHistRes; $max /= $condHistRes;
+ $min /= $condHistRes; $max /= $condHistRes;
for (my($s)=0; $s<$nscans; $s++) {
if ($ants_[$s][$condF] > $max) { $outliers++; $ants_[$s][$condF] = nan; }
if ($ants_[$s][$condF] < $min) { $outliers++; $ants_[$s][$condF] = nan; }
@@ -241,20 +251,24 @@
# edit temperature outliers outside contiguous range
# - Stan's NBP0901 profiles require resolution of 1deg
# - otherwise 0.2deg seems to be fine
+ # - however, on Feb 19, 2016 it was found that the
+ # resolution had been left at 1degC without any
+ # apparent adverse effects
#----------------------------------------------------
$outliers = $modeSamp = 0;
undef(@hist);
for (my($s)=0; $s<$nscans; $s++) {
- next unless ($ants_[$s][$tempF] > 0);
- my($b) = $ants_[$s][$tempF]*1; # 10th of a degree histogram resolution
+ next unless ($ants_[$s][$tempF] >= -10);
+ my($b) = ($ants_[$s][$tempF] + 10) * 0.5;
$hist[$b]++;
next unless ($hist[$b] > $modeSamp);
$modeSamp = $hist[$b]; $modeBin = $b;
}
- printf(STDERR "\n\ttemperature mode: %.1f degC (%d samples)",$modeBin/10,$modeSamp)
- if ($opt_v > 1);
+# printf(STDERR "\n\ttemperature mode: %.1f degC (%d samples)",$modeBin/0.5-10,$modeSamp)
+# if ($opt_v > 1);
($min,$max) = validRange($modeBin);
- $min /= 1; $max /= 1;
+ $min = ($min / 0.5) - 10;
+ $max = ($max / 0.5) - 10;
for (my($s)=0; $s<$nscans; $s++) {
if ($ants_[$s][$tempF] > $max) { $outliers++; $ants_[$s][$tempF] = nan; }
if ($ants_[$s][$tempF] < $min) { $outliers++; $ants_[$s][$tempF] = nan; }
@@ -395,8 +409,8 @@
my($sps) = round(1 / $sampint / $opt_s);
print(STDERR "Creating ${opt_s}Hz time series ($sps samples per bin)...") if ($opt_v);
-&antsAddParams('sampling_frequency',1/$opt_s);
-&antsAddParams('sampling_rate',$opt_s);
+&antsAddParams('sampling_interval',1/$opt_s);
+&antsAddParams('sampling_frequency',$opt_s);
my(@press,@temp,@cond);
my($sp,$np,$st,$nt,$sc,$nc);
--- a/LADCP_w_ocean Tue Feb 02 14:55:24 2016 +0000
+++ b/LADCP_w_ocean Tue Mar 08 14:09:57 2016 +0000
@@ -2,18 +2,18 @@
#======================================================================
# L A D C P _ W _ O C E A N
# doc: Fri Dec 17 18:11:13 2010
-# dlm: Wed Jan 27 20:55:29 2016
+# dlm: Tue Mar 8 14:07:23 2016
# (c) 2010 A.M. Thurnherr
-# uE-Info: 450 0 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 215 54 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# TODO:
+# ! use instrument tilt in sidelobe editing
+# ! worry about water-depth differences (disabled warning)
# - make upcast-flag valid for yoyo casts
# - make diagnostic output 3-beam field work for Earth coordinates
-# - remove water-depth from BT code, which is not really used and a bit of an outlier
+# ? remove water-depth from BT code, which is not really used and a bit of an outlier
# because mean and stddev are used instead of median/mad
-# - read assumed ADCP soundspeed from data file, instead of assuming 1500m/s
-# ? use instrument tilt in sidelobe editing?
$antsSummary = 'calculate vertical velocities from LADCP & CTD time series';
@@ -203,6 +203,16 @@
# - implemented outGrid_firstBin eq '*' (also lastBin)
# Jan 27, 2016: - BUG: outGrid_lastBin eq '*' did not work
# - removed large ref-lr w warning
+# Feb 14, 2016: - fiddled with ping-coherent residuals code, fixing
+# 0-2 potential (minor?) bugs related to outGrid_{first,last}Bin;
+# output of new code checked against old code => identical
+# Feb 16, 2016: - BUG: per-bin-residual QC could cause division by zero
+# Feb 19, 2016: - added -l (disable time-lag filtering)
+# Mar 7, 2016: - added error message when -h is neither number nor file
+# - BUG: -ve depth error message referred to obsolete -d
+# - BUG: dn field name did not use zero filling for year number
+# Mar 8, 2016: - removed L0 water-depth-difference warning
+# - added test for 1500m/s sound speed
# HISTORY END
# CTD REQUIREMENTS
@@ -291,7 +301,7 @@
#------
$antsParseHeader = 0;
-&antsUsage('3:4a:b:c:e:g:h:i:k:m:n:o:p:qs:t:uv:Vw:x:',0,
+&antsUsage('3:4a:b:c:e:g:h:i:k:lm:n:o:p:qs:t:uv:Vw:x:',0,
'[print software -V)ersion] [-v)erbosity <level[1]>]',
'[-q)uick (no single-ping denoising)]',
'[require -4)-beam solutions] [apply beamvel-m)ask <file> if it exists]',
@@ -302,6 +312,7 @@
'[-i)nitial CTD time offset <guestimate> [-u)se as final]]',
'[calculate -n) <lags,lags[10,100]>] [lag -w)indow <sz,sz[240s,20s]>] [lag-p)iece <CTD_elapsed_min|+[,...]>]',
'[require top-3) lags to account for <frac[0.6]> of all]',
+ '[disable time-l)ag filtering]',
'[pressure-sensor -a)cceleration correction <residual/CTD_w_tt>]',
'[-o)utput bin <resolution[20m]>] [-k) require <min[20]> samples]',
'[e-x)ecute <perl-expr>]',
@@ -655,8 +666,14 @@
error("$0: implausibly short cast ($cast_duration seconds)\n")
unless ($cast_duration > 600);
+for (my($ens)=$firstGoodEns; $ens<=$lastGoodEns; $ens++) {
+ croak("$LADCP_file: unexpected sound speed (%d m/s != 1500 m/s) at ensemble \#%d\n",
+ $LADCP{ENSEMBLE}[$ens]->{SPEED_OF_SOUND},$LADCP{ENSEMBLE}[$ens]->{NUMBER})
+ unless ($LADCP{ENSEMBLE}[$ens]->{SPEED_OF_SOUND} == 1500);
+}
+
my($year) = $LADCP{ENSEMBLE}[$firstGoodEns]->{YEAR} % 100;
-&antsAddParams("dn$year",$LADCP{ENSEMBLE}[$firstGoodEns]->{DAYNO});
+&antsAddParams(sprintf('dn%02d',$year),$LADCP{ENSEMBLE}[$firstGoodEns]->{DAYNO});
$LADCP{MEAN_DT} = $cast_duration / ($lastGoodEns-$firstGoodEns-1);
@@ -953,7 +970,7 @@
&& numberp($CTD{DEPTH}[$scan])) {
$LADCP{ENSEMBLE}[$ens]->{REFLR_W_NOSSCORR} = $LADCP{ENSEMBLE}[$ens]->{REFLR_W};
$LADCP{ENSEMBLE}[$ens]->{REFLR_W} *= $CTD{SVEL}[$scan]/1500; # correct for sound-speed variations at source
- error(sprintf("\n$0: negative depth (%.1fm) in CTD file at elapsed(CTD) = %.1fs (use -d?)\n",
+ error(sprintf("\n$0: negative depth (%.1fm) in CTD file at elapsed(CTD) = %.1fs\n",
$CTD{DEPTH}[$scan],$CTD{ELAPSED}[$scan]))
unless ($CTD{DEPTH}[$scan] >= 0);
$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH} = $CTD{DEPTH}[$scan];
@@ -1022,7 +1039,9 @@
$water_depth = &antsFileScanParam(WDF,'water_depth');
close(WDF);
undef($water_depth) unless numberp($water_depth);
- }
+ } else {
+ croak("$0: -h $opt_h defines neither number nor existing file\n");
+ }
}
die("assertion failed (water_depth = $water_depth)")
@@ -1036,8 +1055,8 @@
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(0,sprintf("Large instrument vs. backscatter-derived water-depth difference (%.1fm)\n",$dd))
- if ($dd > 10);
+# warning(0,sprintf("Large instrument vs. backscatter-derived water-depth difference (%.1fm)\n",$dd))
+# if ($dd > 10);
}
if (!$SS_use_BT && !defined($water_depth) && defined($water_depth_BT)) {
warning(1,"using water_depth from ADCP BT data\n");
@@ -1320,8 +1339,8 @@
'max_elapsed',$LADCP{ENSEMBLE}[$realLastGoodEns]->{CTD_ELAPSED});
#------------------------------------------------------------------------------------------------------
-# remove ping-coherent errors
-# - error sources:
+# remove ping-coherent residuals
+# - potential error sources:
# 1) acceleration-dependence of Paroscientific pressure measurements; O(10cm/s) [IWISE 28]
# 2) residual CTD/LADCP mismatch errors; O(1cm/s) [Thurnherr, CWTMC 2011]
# 3) ADCP short-term variability; O(1cm/s) for vertical?
@@ -1334,7 +1353,7 @@
#------------------------------------------------------------------------------------------------------
unless ($opt_q) {
- progress("Removing ping-coherent errors...\n");
+ progress("Removing ping-coherent residuals...\n");
for ($ens=$firstGoodEns; $ens<=$lastGoodEns; $ens++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
@@ -1351,9 +1370,10 @@
: $UPCAST{MEDIAN_W}[$bi]));
}
$LADCP{ENSEMBLE}[$ens]->{MEDIAN_RESIDUAL_W} = median(@residuals); # NB: can be nan!
+ next unless numberp($LADCP{ENSEMBLE}[$ens]->{MEDIAN_RESIDUAL_W});
for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) { # remove from ocean velocities
- next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]); # also ensures that MEDIAN_RESIDUAL_W is not nan
+ next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin] -=
$LADCP{ENSEMBLE}[$ens]->{MEDIAN_RESIDUAL_W};
$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W12}[$bin] -=
@@ -1364,9 +1384,9 @@
if numberp($LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W34}[$bin]);
}
- $LADCP{ENSEMBLE}[$ens]->{REFLR_W} -= # reflr_ocean_w is calculated below
- $LADCP{ENSEMBLE}[$ens]->{MEDIAN_RESIDUAL_W} # NB: this can be set to nan here
- if numberp($LADCP{ENSEMBLE}[$ens]->{REFLR_W})
+ $LADCP{ENSEMBLE}[$ens]->{REFLR_W} -= # NB: this can be nan here
+ $LADCP{ENSEMBLE}[$ens]->{MEDIAN_RESIDUAL_W}
+ if numberp($LADCP{ENSEMBLE}[$ens]->{REFLR_W});
}
#------------------------------------------------------------
@@ -1374,7 +1394,7 @@
progress("\tre-binning profile data...\n");
for (my($bi)=0; $bi<=$#{$DNCAST{ENSEMBLE}}; $bi++) { # bin data
- for (my($i)=0; $i<@{$DNCAST{W}[$bi]}; $i++) {
+ 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};
@@ -1466,8 +1486,8 @@
$uc_n++;
}
}
- local($dc_bres_rms) = sqrt($dc_sumsq/$dc_n);
- local($uc_bres_rms) = sqrt($uc_sumsq/$uc_n);
+ local($dc_bres_rms) = ($dc_n > 0) ? sqrt($dc_sumsq/$dc_n) : nan;
+ local($uc_bres_rms) = ($uc_n > 0) ? sqrt($uc_sumsq/$uc_n) : nan;
progress("\nWriting binned residuals to ");
--- a/LADCP_w_postproc Tue Feb 02 14:55:24 2016 +0000
+++ b/LADCP_w_postproc Tue Mar 08 14:09:57 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: Sat Jan 30 14:50:53 2016
+# dlm: Mon Mar 7 21:03:31 2016
# (c) 2015 A.M. Thurnherr
-# uE-Info: 640 0 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 199 71 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
$antsSummary = 'edit and re-grid LADCP vertical-velocity samples';
@@ -55,6 +55,11 @@
# Jan 25, 2016: - added software version %PARAM
# Jan 26, 2016: - added -v
# Jan 30, 2016: - added -w
+# Feb 14, 2016: - BUG: all bins were off by 1! (-v, inherited limits)
+# Mar 1, 2016: - added required ADCP sampling %PARAMs to dual-head
+# output
+# Mar 7, 2016: - BUG: correlation stats were defined/used for single-head data
+# - removed good_bins() from library as -v allows more control
($ANTS) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
($WCALC) = ($0 =~ m{^(.*)/[^/]*$});
@@ -151,7 +156,6 @@
#
# output_resolution(dz) output_resolution(40)
# bad_range[_dc|_uc](field,min_val,max_val) bad_range_uc('depth',3500,3600)
-# good_bins(first_good,last_good) good_bins(2,'*')
#
#----------------------------------------------------------------------
@@ -192,14 +196,7 @@
push(@brDUc,0);
}
-sub good_bins($$)
-{
- push(@gbFirst,shift); push(@gbLast,shift);
- $gbFirst[$#gbFirst] = 1 if ($gbFirst[$#gbFirst] eq '*');
- $gbLast[$#gbLast] = 9e99 if ($gbLast[$#gbLast] eq '*');
-}
-
-
+#----------------------------------------------------------------------
sub isBad()
{
@@ -376,19 +373,42 @@
'UL_last_valid_bin',&antsRequireParam('outgrid_lastbin'));
}
- for (my($bi)=0; $bi<=$#dcw1; $bi++) { # calc DL median profile
+ for (my($bi)=0; $bi<=$#dcw1; $bi++) { # calc DL median profile (before reading UL data)
$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}))
- unless (round($pulse_length) == round($P{ADCP_pulse_length}));
- $blanking_dist = sprintf('%d & %d',round($blanking_dist),($P{ADCP_blanking_distance}))
- unless (round($blanking_dist) == round($P{ADCP_blanking_distance}));
+ #
+ # ADCP bin length, pulse length, and blanking distance for dual head casts
+ # with inconsistent values:
+ # bin length: use smaller value, which will lead to smaller spectral correction
+ # pulse length: same
+ # blanking distance: use smaller value, which is conservative e.g. for filters for ringing
+ #
+ my($warned);
+ unless (round($bin_length) == round($P{ADCP_bin_length})) {
+ unless ($warned) {
+ &antsInfo("WARNING: inconsistent ADCP sampling parameters --- using conservative values");
+ $warned = 1;
+ }
+ $bin_length = min($bin_length,$P{ADCP_bin_length});
+ }
+ unless (round($pulse_length) == round($P{ADCP_pulse_length})) {
+ unless ($warned) {
+ &antsInfo("WARNING: inconsistent ADCP sampling parameters --- using conservative values");
+ $warned = 1;
+ }
+ $pulse_length = min($pulse_length,$P{ADCP_pulse_length});
+ }
+ unless (round($blanking_dist) == round($P{ADCP_blanking_distance})) {
+ unless ($warned) {
+ &antsInfo("WARNING: inconsistent ADCP sampling parameters --- using conservative values");
+ $warned = 1;
+ }
+ $blanking_dist = min($blanking_dist,$P{ADCP_blanking_distance});
+ }
- $PROF = $STN = $ID = $id; $RUN = antsRequireParam('run_label');
+ $PROF = $STN = $ID = $id; $RUN = antsRequireParam('run_label'); # set variables for editing
undef(@rngMin); undef(@rngMax); undef(@bins);
unless ($return = do "./EditParams") { # man perlfunc
croak("./EditParams: $@\n") if ($@);
@@ -408,11 +428,11 @@
} # 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);
+ next if ($ants_[0][$bF]<$fvBin || # => apply them
+ $ants_[0][$bF]>$lvBin);
} 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} || # => use values from [LADCP_w_ocean]
+ $ants_[0][$bF]>$P{outgrid_lastbin});
}
$filt++,next if &isBad(); # additional editing
@@ -477,6 +497,8 @@
&antsAddParams('outgrid_dz',$opt_o,'outgrid_minsamp',$opt_k);
&antsAddParams('min_depth',round($min_depth),'max_depth',round($max_depth));
&antsAddParams('water_depth',$P{water_depth},'water_depth.sig',$P{water_depth.sig});
+ &antsAddParams('ADCP_bin_length',$bin_length,'ADCP_pulse_length',$pulse_length);
+ &antsAddParams('ADCP_blanking_distance',$blanking_dist);
undef($antsOldHeaders);
} else {
@antsNewLayout = ('depth','hab',
@@ -510,28 +532,32 @@
}
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);
+ if ($dual_head) {
+ GMT_psxy('-W2/100/100/255'); # surface layer limit
+ print(GMT "-0.1 $opt_s\n0.07 $opt_s\n");
+ if ($dc_R < 0.3 || !numberp($dc_R)) { # correlation statistics
+ &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);
+ 0.02,0.9*$ymax,100*$dc_esig,100*$dc_rsig) if ($dual_head);
+ if ($dual_head) {
+ 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);
+ 0.02,0.95*$ymax,100*$uc_esig,100*$uc_rsig) if ($dual_head);
GMT_setR($R);
--- a/LADCP_wspec Tue Feb 02 14:55:24 2016 +0000
+++ b/LADCP_wspec Tue Mar 08 14:09:57 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: Mon Jan 25 09:33:54 2016
+# dlm: Tue Mar 1 21:27:41 2016
# (c) 2012 A.M. Thurnherr
-# uE-Info: 39 0 NIL 0 0 72 10 2 4 NIL ofnI
+# uE-Info: 26 29 NIL 0 0 72 10 2 4 NIL ofnI
#======================================================================
$antsSummary = 'calculate VKE window spectra from LADCP profiles';
@@ -20,6 +20,10 @@
# Oct 12, 2015: - require ANTSlibs V6.2 for release
# Oct 13, 2015: - adapted to [version.pl]
# Jan 25, 2016: - added software version %PARAM
+# Mar 1, 2016: - made trailing message much less frequent
+# - BUG: program croaked gracelessly or entered infinite loop
+# when presented with insufficient (no valid) input
+# data
($ANTSBIN) = (`which list` =~ m{^(.*)/[^/]*$});
($ANTSLIB) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
@@ -104,15 +108,19 @@
&antsInstallBufFull(0); # read entire file
&antsIn();
-while ($ants_[0][$zfnr] < $opt_s) { # remove surface layer
+while (@ants_ && $ants_[0][$zfnr] < $opt_s) { # remove surface layer
shift(@ants_);
}
+croak("$0: insufficient data (no valid records found)\n")
+ unless (@ants_);
-for ($trimmed=0; !numberp($ants_[$#ants_][$dcwfnr]) && !numberp($ants_[$#ants_][$ucwfnr]); $trimmed++) {
+for ($trimmed=0; @ants_ && !numberp($ants_[$#ants_][$dcwfnr]) && !numberp($ants_[$#ants_][$ucwfnr]); $trimmed++) {
pop(@ants_);
}
&antsInfo("$trimmed trailing non-numeric records trimmed")
- if ($trimmed);
+ if ($trimmed > 1); # 1 is very common
+croak("$0: insufficient data (no valid records found)\n")
+ unless (@ants_);
$dz = &antsXCheck($zfnr,0,$#ants_,1.01); # calc dT; 1% jitter
&antsAddParams("wspec::input_depth_resolution",$dz);
--- a/time_lag.pl Tue Feb 02 14:55:24 2016 +0000
+++ b/time_lag.pl Tue Mar 08 14:09:57 2016 +0000
@@ -1,9 +1,9 @@
#======================================================================
# T I M E _ L A G . P L
# doc: Fri Dec 17 21:59:07 2010
-# dlm: Tue Jan 26 15:04:54 2016
+# dlm: Mon Mar 7 18:31:34 2016
# (c) 2010 A.M. Thurnherr
-# uE-Info: 68 33 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 71 68 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -66,6 +66,9 @@
# - BUG: time-lag plot was not produced when final lag piece had problems
# Jan 25, 2016: - search-radius-doubling heuristic had typo
# - added %PARAMs
+# Feb 19, 2016: - added support for -l
+# - added warning
+# Mar 7, 2016: - BUG: editing did not work correctly in all cases
# DIFFICULT STATIONS:
# NBP0901#131 this requires the search-radius doubling heuristic
@@ -331,20 +334,23 @@
# of uncertain time lags
#--------------------------------------------------
- if ($scan_increment == 1) {
+ if ($scan_increment == 1 && !$opt_l) {
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]]);
+ push(@elim,($fg[$i] == 0) ? $LADCP{ENSEMBLE}[$firstGoodEns]->{ELAPSED} : $elapsed[$fg[$i]],
+ ($lg[$i] == $n_windows-1) ? $LADCP{ENSEMBLE}[$lastGoodEns]->{ELAPSED} : $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);
+ my($pct) = round(100*$nerm/($last_ens-$first_ens+1));
+ progress("\t\t$nerm ensembles removed ($pct%% of total), leaving %d run(s)\n",scalar(@elim)/2);
+ warning(1,"time-lag editing removed large fraction of samples (%d%% of total)\n",$pct)
+ if ($pct > 30);
}
#------------------------------------------------------