--- a/LADCPproc
+++ b/LADCPproc
@@ -2,9 +2,9 @@
#======================================================================
# L A D C P P R O C
# doc: Thu Sep 16 20:36:10 2010
-# dlm: Sun Jul 27 16:23:50 2014
+# dlm: Tue Aug 5 14:39:26 2014
# (c) 2010 A.M. Thurnherr
-# uE-Info: 103 69 NIL 0 0 72 11 2 4 NIL ofnI
+# uE-Info: 410 82 NIL 0 0 72 10 2 4 NIL ofnI
#======================================================================
# NOTES:
@@ -101,6 +101,9 @@
# May 20, 2014: - merged laptop with whoosher versions (folded in Feb 22 change)
# Jul 27, 2014: - renamed LADCPproc.UHcode to LADCPproc.shearmethod, because the code has
# diverged more and more from the UH implementation
+# - added -v to allow calculation of package velocity
+# Aug 5, 2014: - BUG: (bad one, too): ref_lr_w called from mk_prof had edited some of the
+# horizontal velocity data, which were nevertheless used later on!!!
($ANTS) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
($PERL_TOOLS) = (`which mkProfile` =~ m{^(.*)/[^/]*$});
@@ -125,7 +128,7 @@
$antsSummary = "$version -- process LADCP data to get shear, time series";
$antsParseHeader = 0;
-&antsUsage('24a:b:c:df:g:i:kl:n:o:p:rs:t:u:w:z',2,
+&antsUsage('24a:b:c:df:g:i:kl:n:o:p:rs:t:u:v:w:z',2,
'[use -2)dary CTD sensor pair]',
'[require -4)-beam LADCP solutions]',
'[use -r)DI bottom-track data]',
@@ -135,6 +138,7 @@
'[-i)nitial LADCP time lag <guestimate>]',
'[-l)ag LADCP <by>] [auto-lag -w)indow <size[120s]>] [-n) <auto-lag windows[20]]',
'[correct echo amplitude -u)sing D[eines99]|V[isbeck04]|T[hurnherr11]|n[ocorr]',
+ '[ocean -v)elocity <file> for calculating package velocity]',
'[-d)iagnostic screen output] [-z)oom through problems]',
'output: [shear-p)rofile <file>] [-t)ime series <file>] [-f)lag <file>] [-b)ottom-track <file>]',
' [-a)coustic backscatter <dts-file] [bottom-trac-k) profs]',
@@ -171,7 +175,7 @@
croak("$0: cannot decode -c $opt_c\n")
unless numberp($CC_offset) && numberp($CC_cos_fac) && numberp($CC_sin_fac);
}
-
+
$LADCP_file = &antsFileArg();
$CTD_file = &antsFileArg();
@@ -221,6 +225,34 @@
&antsAddParams('LADCPproc::vertical_resolution',$GRID_DZ);
#----------------------------------------------------------------------
+# Step 2a: Load Ocean Velocity Profile
+#----------------------------------------------------------------------
+
+if (defined($opt_v)) { # load velocity profile
+ print(STDERR "Reading ocean-velocity profile ($opt_v)...");
+ if (open(OVF,$opt_v)) {
+ @ovl = &antsFileLayout(OVF);
+ $ovdF = &localFnr('depth',@ovl);
+ $ovuF = &localFnr('u',@ovl);
+ $ovvF = &localFnr('v',@ovl);
+ my(@ov) = &antsFileIn(OVF);
+ my($first_depth) = $ov[$ovdF]; my($last_depth);
+ do {
+ push(@ovu,$ov[$ovuF]);
+ push(@ovv,$ov[$ovvF]);
+ $last_depth = $ov[$ovdF];
+ } while (@ov = &antsFileIn(OVF));
+ close(OVF);
+ croak("$opt_v: incompatible depth grid\n")
+ unless (($last_depth-$first_depth) == $#ovu*$GRID_DZ);
+ printf(STDERR "\n\t%d velocities",scalar(@ovv)) if ($opt_d);
+ } else {
+ printf(STDERR "\n\n\t\tWARNING: $opt_v: $!\n");
+ }
+ print(STDERR "\n");
+}
+
+#----------------------------------------------------------------------
# Step 3: Read CTD data
#----------------------------------------------------------------------
@@ -328,6 +360,7 @@
@{$LADCP{ENSEMBLE}[$ens]->{BEAM_VEL}[$bin]} = @{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]};
@{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]} =
velInstrumentToEarth(\%LADCP,$ens,velBeamToInstrument(\%LADCP,@{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]}));
+# print(STDERR "<$ens,$bin>: @{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]}\n");
@{$LADCP{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$bin]} = # fake it to fool ref_lr_w
(0,0,0,defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W]) ? 100 : 0);
}
@@ -373,8 +406,8 @@
print(STDERR "\t\tconstructing profile time series...")
if ($opt_d);
-($LADCP_start,$LADCP_end,$LADCP_bottom,$w_gap_time,$zErr,$maxz) =
- mk_prof(\%LADCP,0,undef,1,6,70,0.1,$LADCP_max_gap);
+($LADCP_start,$LADCP_end,$LADCP_bottom,$w_gap_time,$zErr,$maxz) = # NB: chose parameters to avoid editing by
+ mk_prof(\%LADCP,0,undef,1,6,0,0.1,$LADCP_max_gap,0); # ref_lr_w
croak("\n$LADCP_file: no good ensembles found\n")
unless defined($LADCP_start);
@@ -559,6 +592,8 @@
print(STDERR "Calculating shear profiles & producing time-depth-series (.tds) output...");
@antsNewLayout = ('ens','elapsed','CTD_depth','downcast','CTD_lat','CTD_lon','depth','u_z','v_z','w_z');
+push(@antsNewLayout,'pkg_u','pkg_v')
+ if defined($opt_v);
#--------------------------------------------------------------------------------
# callback routine to output .tds data, called once each for down-/upcasts after
@@ -580,23 +615,31 @@
'LADCPproc::min_depth',depthOfGI($mingi),'LADCPproc::max_depth',depthOfGI($#ens_vals));
for (my($gi)=0; $gi<@ush_vals; $gi++) {
for (my($i)=0; $i<@{$ush_vals[$gi]}; $i++){
- &antsOut($ens_vals[$gi][$i],
- $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{ELAPSED_TIME}+$CTD{first_elapsed}-$opt_l,
- $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{DEPTH},$downcast,
- $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{CTD_LAT},
- $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{CTD_LON},
- depthOfGI($gi),$ush_vals[$gi][$i],$vsh_vals[$gi][$i],$wsh_vals[$gi][$i]);
+ my(@out) = ($ens_vals[$gi][$i],
+ $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{ELAPSED_TIME}+$CTD{first_elapsed}-$opt_l,
+ $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{DEPTH},$downcast,
+ $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{CTD_LAT},
+ $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{CTD_LON},
+ depthOfGI($gi),$ush_vals[$gi][$i],$vsh_vals[$gi][$i],$wsh_vals[$gi][$i]);
+ push(@out,$LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{PACKAGE_VELOCITY}[$U],
+ $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{PACKAGE_VELOCITY}[$V])
+ if defined($opt_v);
+ &antsOut(@out);
}
}
} else {
for (my($gi)=$#ush_vals; $gi>=0; $gi--) {
for (my($i)=0; $i<@{$ush_vals[$gi]}; $i++) {
- &antsOut($ens_vals[$gi][$i],
- $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{ELAPSED_TIME}+$CTD{first_elapsed}-$opt_l,
- $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{DEPTH},$downcast,
- $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{CTD_LAT},
- $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{CTD_LON},
- depthOfGI($gi),$ush_vals[$gi][$i],$vsh_vals[$gi][$i],$wsh_vals[$gi][$i]);
+ my(@out) = ($ens_vals[$gi][$i],
+ $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{ELAPSED_TIME}+$CTD{first_elapsed}-$opt_l,
+ $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{DEPTH},$downcast,
+ $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{CTD_LAT},
+ $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{CTD_LON},
+ depthOfGI($gi),$ush_vals[$gi][$i],$vsh_vals[$gi][$i],$wsh_vals[$gi][$i]);
+ push(@out,$LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{PACKAGE_VELOCITY}[$U],
+ $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{PACKAGE_VELOCITY}[$V])
+ if defined($opt_v);
+ &antsOut(@out);
}
}
}
--- a/LADCPproc.BT
+++ b/LADCPproc.BT
@@ -1,9 +1,9 @@
#======================================================================
# L A D C P P R O C . B T
# doc: Wed Oct 20 21:05:37 2010
-# dlm: Thu Sep 19 13:17:41 2013
+# dlm: Tue Aug 5 14:38:08 2014
# (c) 2010 A.M. Thurnherr
-# uE-Info: 119 96 NIL 0 0 72 10 2 4 NIL ofnI
+# uE-Info: 205 0 NIL 0 0 72 10 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -15,6 +15,10 @@
# - changed from echo amplitude to Sv
# May 16, 2012: - added support for -r)DI BT data
# Sep 19, 2013: - added support for $BT_range_method
+# Aug 5, 2014: - BUG: invalid velocities were used; see July 15 bug fix
+# in [LADCPproc.shearmethod]
+# - artifically removed BT profile data apparently below
+# seabed but that had passed previous tests
my($BEAM1) = 0;
my($BEAM2) = 1;
@@ -174,7 +178,7 @@
if ($opt_k) {
for (my($bin)=$BT_bin_start-1; $bin<$LADCP{N_BINS}; $bin++) {
- next if ($edit_flags[$ens][$bin]);
+ next if ($edit_flags[$ens][$bin] || $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W]==0);
printf(BTF "%d %d %d %f %f %f %f %f %f %f %f %f %f %f\n",
$LADCP{ENSEMBLE}[$ens]->{NUMBER},
depthOfBin($ens,$bin),$LADCP{ENSEMBLE}[$ens]->{DEPTH},
@@ -191,13 +195,14 @@
}
for (my($bin)=$BT_bin_start-1; $bin<$LADCP{N_BINS}; $bin++) {
- next if ($edit_flags[$ens][$bin]);
- my($gi) = depthOfBin($ens,$bin) / $GRID_DZ;
+ next if ($edit_flags[$ens][$bin] || $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W]==0);
+ my($dob) = depthOfBin($ens,$bin);
+ next if ($dob > $water_depth);
+ my($gi) = int($dob / $GRID_DZ);
push(@{$BTu_vals[$gi]},$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$U]-$CTD_u);
push(@{$BTv_vals[$gi]},$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$V]-$CTD_v);
push(@{$BTw_vals[$gi]},$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W]-$CTD_w);
}
-
}
sub getBTprof($$)
--- a/LADCPproc.backscatter
+++ b/LADCPproc.backscatter
@@ -1,9 +1,9 @@
#======================================================================
# L A D C P P R O C . B A C K S C A T T E R
# doc: Wed Oct 20 13:02:27 2010
-# dlm: Thu Mar 27 18:54:02 2014
+# dlm: Tue Aug 5 13:32:13 2014
# (c) 2010 A.M. Thurnherr
-# uE-Info: 25 52 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 30 71 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -23,6 +23,11 @@
# for shallow casts)
# Mar 21, 2014: - adapted to new [LADCPproc.utils]
# Mar 27, 2014: - adapted to depthOfBinAlongBeam()
+# Jul 27, 2014: - moved depthOfGI() to [LADCPproc.utils]
+# Aug 5, 2014: - BUG: find_backscatter_seabed() discarded everything if
+# LADCP bin 1 had the backscatter max at the edge of
+# the search domain, which could easily happen when
+# the bottom stop was a long way from the seabed
my($BEAM1) = 0;
my($BEAM2) = 1;
@@ -243,8 +248,6 @@
} # for $end
} # sub
-sub depthOfGI($) { return $_[0]*$GRID_DZ + $GRID_DZ/2; } # depth corresponding to particular grid index
-
sub find_backscatter_seabed($)
{
my($water_depth) = @_;
@@ -252,7 +255,8 @@
my($search_below) = max(0,$water_depth-$BT_begin_search_above);
my($mdgi) = int($search_below/$GRID_DZ); # grid index to begin search
- print(STDERR "\n\t\tlooking for seabed below $search_below m (gi >= $mdgi)") if ($opt_d);
+ printf(STDERR "\n\t\tlooking for seabed below %d m (gi = [%d..%d])",$search_below,$mdgi,scalar(@nSv))
+ if ($opt_d);
print(STDERR "\n\t\tseabed-max grid indices:") if ($opt_d);
@@ -266,19 +270,13 @@
$min = $avg if ($avg < $min);
$max = $avg, $gimax = $gi if ($avg > $max);
}
- if ($gimax == $firstvalid) { # should be robust except maybe for huge bins
- printf(STDERR "\n\t\tdata from below-seabed bins (%d-%d) discarded",
- $bin+1,$LADCP{N_BINS}) if ($opt_d);
- last;
- }
- if ($max-$min>10 && $gimax!=$lastvalid) { # ignore boundary maxima & scatter
- printf(STDERR " %d",$gimax-$mdgi) if ($opt_d);
+ if ($max-$min>10 && $gimax!=$firstvalid && $gimax!=$lastvalid) { # ignore boundary maxima & scatter
+ printf(STDERR " %d",$gimax) if ($opt_d);
push(@wdepth_gi,$gimax);
}
}
return (depthOfGI(avg(@wdepth_gi)),stddev(@wdepth_gi)*$GRID_DZ);
-
}
1;
--- a/LADCPproc.shearmethod
+++ b/LADCPproc.shearmethod
@@ -1,9 +1,9 @@
#======================================================================
# L A D C P P R O C . S H E A R M E T H O D
# doc: Fri Sep 17 20:27:53 2010
-# dlm: Sun Jul 27 16:25:56 2014
+# dlm: Sun Jul 27 19:44:41 2014
# (c) 2010 A.M. Thurnherr & E. Firing
-# uE-Info: 314 0 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 338 65 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# PERLified functions from E. Firing's [merge.c], modified by A.M. Thurnherr
@@ -54,6 +54,7 @@
# Jul 27, 2014: - renamed from [LADCPproc.UHmethod] because the code has
# diverged quite a bit from the original mplementation. However,
# the shear editing core remains as in the UH code
+# - added calculation of PACKAGE_VELOCITY (-v)
#======================================================================
# VELOCITY EDITING
@@ -314,25 +315,37 @@
#==============================================================================
#----------------------------------------------------------------------
-# uv_to_shear(ens)
+# uv_to_shear(ens,calc_pkg_vel)
# - sets @sh_i0, @sh_i1, @dsh, @ush, @vsh, @wsh for a given ensemble
+# - also sets $LADCP{ENSEMBLE}[$ens]->{PACKAGE_VELOCITY} on -v
#----------------------------------------------------------------------
-sub uv_to_shear($)
+sub uv_to_shear($$)
{
- my($ens) = @_;
+ my($ens,$calc_pkg_vel) = @_;
my($nvel) = 0;
+ my(@pu,@pv);
@sh_i0 = @sh_i1 = ();
- for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
+ for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) { # select valid velocities
next if ($edit_flags[$ens][$bin] || $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W]==0);
$nvel++;
push(@sh_i1,$bin) if (@sh_i0);
push(@sh_i0,$bin) if ($bin < $LADCP{N_BINS}-1);
+ if (defined($calc_pkg_vel)) {
+ my($ovgi) = int(&depthOfBin($ens,$bin) / $GRID_DZ);
+ push(@pu,$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$U]-$ovu[$ovgi]);
+ push(@pv,$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$V]-$ovv[$ovgi]);
+ }
+ }
+
+ if (defined($calc_pkg_vel)) { # calculate package velocity
+ $LADCP{ENSEMBLE}[$ens]->{PACKAGE_VELOCITY}[$U] = avg(@pu);
+ $LADCP{ENSEMBLE}[$ens]->{PACKAGE_VELOCITY}[$V] = avg(@pv);
}
@dsh = ();
- for (my($i)=0; $i<@sh_i1; $i++) { # calc and bin shears
+ for (my($i)=0; $i<@sh_i1; $i++) { # calc and bin shears
my($d0) = &depthOfBin($ens,$sh_i0[$i]);
my($d1) = &depthOfBin($ens,$sh_i1[$i]);
next unless ($d0>=0 && $d1>=0);
@@ -415,10 +428,10 @@
local(@sh_i0,@sh_i1);
local(@dsh,@ush,@vsh,@wsh,@ens);
- uv_to_shear($ens);
+ uv_to_shear($ens,0);
if ($edit_shear) {
set_shear_flag($ens);
- $nvel += uv_to_shear($ens);
+ $nvel += uv_to_shear($ens,defined($opt_v));
$nsh += @dsh;
}
--- a/LADCPproc.utils
+++ b/LADCPproc.utils
@@ -1,27 +1,23 @@
#======================================================================
# L A D C P P R O C . U T I L S
# doc: Fri Mar 21 15:16:59 2014
-# dlm: Thu Mar 27 19:36:29 2014
+# dlm: Sun Jul 27 17:00:57 2014
# (c) 2014 A.M. Thurnherr
-# uE-Info: 57 38 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 14 49 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
# Mar 21, 2014: - created
# - added rangeToBin()
# Mar 27, 2014: - added rangeToBinAlongBeam()
-
-sub rangeToBin($$)
-{
- my($ens,$bin) = @_;
- my($sscorr) = $LADCP{ENSEMBLE}[$ens]->{CTD_SVEL} / $LADCP{ENSEMBLE}[$ens]->{SPEED_OF_SOUND};
- return $sscorr * ($LADCP{DISTANCE_TO_BIN1_CENTER} + $bin*$LADCP{BIN_LENGTH}) / cos(rad($LADCP{BEAM_ANGLE}));
-}
+# Jul 27, 2014: - improved comments
+# - moved depthOfGI() here from [LADCPproc.backscatter]
#----------------------------------------------------------------------
-# - in the UH code
-# - the distance to the first bin is not soundspeed-corrected
-# - instrument tilt is not considered
+# calculate depth of particular bin in particular ensemble
+# in contrast to the original UH code:
+# - the distance to the first bin is soundspeed-corrected
+# - instrument tilt is considered
#----------------------------------------------------------------------
sub depthOfBin($$)
@@ -40,6 +36,24 @@
$LADCP{ENSEMBLE}[$ens]->{DEPTH} + &dzToBin($ens,$bin);
}
+#----------------------------------------------------------------------
+# calculate along-beam distance between transducer and center of
+# particular bin in particular ensemble
+# - used for acoustic backscatter correction
+#----------------------------------------------------------------------
+
+sub rangeToBin($$)
+{
+ my($ens,$bin) = @_;
+ my($sscorr) = $LADCP{ENSEMBLE}[$ens]->{CTD_SVEL} / $LADCP{ENSEMBLE}[$ens]->{SPEED_OF_SOUND};
+ return $sscorr * ($LADCP{DISTANCE_TO_BIN1_CENTER} + $bin*$LADCP{BIN_LENGTH}) / cos(rad($LADCP{BEAM_ANGLE}));
+}
+
+#----------------------------------------------------------------------------
+# calculate depth of particular bin of particular beam in particular ensemble
+# - used to map acoustic backscatter of different beams correctly when
+# instrument tilt is large
+#----------------------------------------------------------------------------
sub depthOfBinAlongBeam($$$)
{
@@ -67,6 +81,21 @@
$LADCP{ENSEMBLE}[$ens]->{DEPTH} + &dzToBinAlongBeam($ens,$bin,$beam);
}
+#---------------------------------------------------------------------------
+# return center depth corresponding to particular grid index (in shear grid)
+#---------------------------------------------------------------------------
+
+sub depthOfGI($) { return $_[0]*$GRID_DZ + $GRID_DZ/2; } # depth corresponding to particular grid index
+
+#----------------------------------------------------------------------
+# return ocean velocity (u,v,w) at a given depth
+#----------------------------------------------------------------------
+
+sub oceanVel($)
+{
+
+}
+
1;