--- a/LADCPintsh
+++ b/LADCPintsh
@@ -2,9 +2,9 @@
#======================================================================
# L A D C P I N T S H
# doc: Thu Oct 14 21:22:50 2010
-# dlm: Wed Feb 16 15:13:46 2011
+# dlm: Thu Jul 7 06:30:50 2011
# (c) 2010 A.M. Thurnherr & E. Firing
-# uE-Info: 401 0 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 32 28 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
$antsSummary = 'integrate LADCP shear';
@@ -29,20 +29,25 @@
# samples (i.e. primarily when chosing a small -r)
# Dec 9, 2010: - allowed for empty BT file
# Feb 16, 2011: - BUG: gaps in shear data were not handled correctly in baroclinic solution
+# Jul 7, 2011: - added -m
($ANTS) = (`which list` =~ m{^(.*)/[^/]*$});
require "$ANTS/ants.pl";
require "$ANTS/libstats.pl";
-&antsUsage('b:dr:s:w:u:',0,
+&antsUsage('b:dm:r:s:w:u:',0,
'[-d)ebug]',
- '[reference with -b)ottom-track <file>]',
+ '[reference with -b)ottom-track <file> [-m)in BT samp[10]]]',
'[-u)plooker <file>]',
'[min -r)aw <shear samp[10]>] [min -s)moothed <shear samp[20]>]',
'[-w)eighted avg <below[1200]>]',
'[LADCP shear file]');
+croak("$0: -m requires -b\n")
+ if defined($opt_m) && !defined($opt_b);
+
&antsCardOpt(\$opt_r,10); ## minimum number of samples for shear
+&antsCardOpt(\$opt_m,10); # minimum number of samples for BT data
&antsCardOpt(\$opt_s,20); ## minimum number of samples for smoothed shear
&antsCardOpt(\$opt_w,1200); ## probably 1200 for midlatitude-equatorial, 2400 for high lat
@@ -311,7 +316,7 @@
while (my(@BTr) = &antsFileIn(BTF)) {
my($gi) = int($BTr[$BTdF] / $DZ);
- next unless ($BTr[$BTndF] >= $opt_r);
+ next unless ($BTr[$BTndF] >= $opt_m);
$BT_nsamp[$gi] = $BTr[$BTndF];
$BT_u[$gi] = $BTr[$BTuF];
$BT_v[$gi] = $BTr[$BTvF];
--- 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: Wed Jun 15 15:47:25 2011
+# dlm: Fri Jul 8 03:33:21 2011
# (c) 2010 A.M. Thurnherr & E. Firing
-# uE-Info: 479 55 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 491 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
$antsSummary = 'process LADCP data to get shear, time series';
@@ -37,6 +37,11 @@
# Jan 22, 2011: - added -g) lat,lon
# - added -c)ompass corr
# Jun 15, 2011: - added mean backscatter profile to default output
+# Jul 7, 2011: - added support for $BT_* processing parameters
+# - replaced old per-bin acoustic backscatter profile output by
+# acoustic backscatter depth-time-series
+# - disabled seabed and acoustic backscatter code when not required (e.g. UL)
+# - made non-diagnostic output terser
($ANTS) = (`which list` =~ m{^(.*)/[^/]*$});
($PERL_TOOLS) = (`which mkProfile` =~ m{^(.*)/[^/]*$});
@@ -55,7 +60,7 @@
require "$PERL_TOOLS/RDI_Utils.pl";
$antsParseHeader = 0;
-&antsUsage('24ab:c:df:g:i:kl:n:o:ps:t:w:',2,
+&antsUsage('24a:b:c:df:g:i:kl:n:o:ps:t:w:',2,
'[use -2)dary CTD sensor pair]',
'[require -4)-beam LADCP solutions]',
'[-s)etup <file>] [-g)ps <lat,lon>]',
@@ -66,7 +71,7 @@
'[-l)ag LADCP <by>] [auto-lag -w)indow <size[120s]>] [-n) <auto-lag windows[20]]',
'[-d)iagnostic output]',
'output: [-t)ime series <file>] [-f)lag <file>] [-b)ottom-track <file>]',
- ' [per-bin -a)coustic backscatter profiles] [bottom-trac-k) profs]',
+ ' [-a)coustic backscatter <depth-time-series file] [bottom-trac-k) profs]',
'<RDI file> <SeaBird file>');
$RDI_Coords::minValidVels = 4 if ($opt_4);
@@ -107,7 +112,8 @@
print(STDERR "Reading LADCP data ($LADCP_file)...");
readData($LADCP_file,\%LADCP);
-printf(STDERR "\n\t%d ensembles\n",scalar(@{$LADCP{ENSEMBLE}}));
+printf(STDERR "\n\t%d ensembles",scalar(@{$LADCP{ENSEMBLE}})) if ($opt_d);
+print(STDERR "\n");
#----------------------------------------------------------------------
# Step 2: Set Processing Parameters
@@ -129,10 +135,11 @@
$wbin_start = 2 unless ($wbin_start > 2);
$ubin_start = 2 unless ($ubin_start > 2);
$shbin_start = 2 unless ($shbin_start > 2);
+ $BT_bin_start = 2 unless ($BT_bin_start > 2);
}
&antsAddParams('ADCP_orientation',
- $dta->{ENSEMBLE}[0]->{XDUCER_FACING_UP} ? 'uplooker' : 'downlooker');
+ $LADCP{ENSEMBLE}[0]->{XDUCER_FACING_UP} ? 'uplooker' : 'downlooker');
$SHEAR_PREGRID_DZ = 20;
$GRID_DZ = defined($opt_o) ? $opt_o : 5;
@@ -153,14 +160,15 @@
print(STDERR "Reading CTD data ($CTD_file)...");
readCTD($CTD_file,\%CTD);
-printf(STDERR "\n\t%d scans\n",scalar(@{$CTD{press}}));
+printf(STDERR "\n\t%d scans",scalar(@{$CTD{press}})) if ($opt_d);
+print(STDERR "\n");
#----------------------------------------------------------------------
# Step 4: Pre-Process CTD & LADCP Data
#----------------------------------------------------------------------
printf(STDERR "Pre-processing data...");
-printf(STDERR "\n\tCTD...");
+printf(STDERR "\n\tCTD...") if ($opt_d);
#------------------------
# clean CTD pressure data
@@ -187,13 +195,11 @@
printf(STDERR "\n\t\tmax pressure [%ddbar] at scan#%d",$CTD{maxpress},$CTD{atbottom})
if $opt_d;
-print(STDERR "\n");
-
#----------------------------------------------------------------------
# Step 4b: Pre-Process LADCP Data
#----------------------------------------------------------------------
-print(STDERR "\tLADCP...");
+print(STDERR "\n\tLADCP...") if ($opt_d);
#-------------------------------------------
# transform to earth coordinates if required
@@ -287,8 +293,8 @@
print(STDERR "Associating CTD data with LADCP ensembles...");
-for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
- my($lastSvel);
+for (my($min_depth)=9e99,my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
+ my($lastSvel);
my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} - $opt_l) / $CTD{sampint});
if ($r < 0 && $ens == $LADCP_start) {
$r = int(($LADCP{ENSEMBLE}[++$ens]->{ELAPSED_TIME} - $opt_l) / $CTD{sampint})
@@ -299,13 +305,18 @@
}
if ($r > $#{$CTD{press}}) {
printf(STDERR "\n\tCTD data end while instrument is still in water => truncating %ds of LADCP data",
- $LADCP{ENSEMBLE}[$LADCP_end]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME});
+ $LADCP{ENSEMBLE}[$LADCP_end]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME})
+ if ($opt_d);
$LADCP_end = $ens - 1;
last;
}
my($dr);
for ($dr=0; !numberp($CTD{press}[$r+$dr]); $dr--) {}
$LADCP{ENSEMBLE}[$ens]->{DEPTH} = depth($CTD{press}[$r+$dr],$CTD{lat});
+ if ($LADCP{ENSEMBLE}[$ens]->{DEPTH} < $min_depth) {
+ $min_depth = $LADCP{ENSEMBLE}[$ens]->{DEPTH};
+ $LADCP_top = $ens;
+ }
$LADCP{ENSEMBLE}[$ens]->{CTD_W} = $CTD{w}[$r];
$LADCP{ENSEMBLE}[$ens]->{CTD_TEMP} = $CTD{temp}[$r];
$LADCP{ENSEMBLE}[$ens]->{CTD_SVEL} = sVel($CTD{salin}[$r],$CTD{temp}[$r],$CTD{press}[$r+$dr]);
@@ -323,7 +334,8 @@
}
}
-&antsAddParams('bottom_depth',round($LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH}),
+&antsAddParams('top_depth',,round($LADCP{ENSEMBLE}[$LADCP_top]->{DEPTH}),
+ 'bottom_depth',round($LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH}),
'start_date',$LADCP{ENSEMBLE}[$LADCP_start]->{DATE},
'start_time',$LADCP{ENSEMBLE}[$LADCP_start]->{TIME},
'bottom_date',$LADCP{ENSEMBLE}[$LADCP_bottom]->{DATE},
@@ -337,53 +349,62 @@
# Step 6: Calculate Acoustic Backscatter Profile
#----------------------------------------------------------------------
-print(STDERR "Finding seabed in acoustic backscatter profiles...");
-
+print(STDERR "Calculating acoustic backscatter profiles...");
mk_backscatter_profs($LADCP_start,$LADCP_end);
-($water_depth,$sig_water_depth) =
- find_backscatter_seabed($LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH});
-
-$min_hab = $water_depth - $LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH};
-printf(STDERR "\n\twater depth = %d(+-%.1f)m",$water_depth,$sig_water_depth);
-printf(STDERR "\n\tclosest approach = %dmab",$min_hab);
-
print(STDERR "\n");
#----------------------------------------------------------------------
# Step 7: Find Seabed
#----------------------------------------------------------------------
-print(STDERR "Finding seabed in BT data...");
+if ($LADCP{ENSEMBLE}[$LADCP_start]->{XDUCER_FACING_DOWN}) {
+
+ print(STDERR "Finding seabed...");
-($BT_water_depth,$sig_BT_water_depth) =
- find_seabed(\%LADCP,$LADCP_bottom,$LADCP{BEAM_COORDINATES});
-
-if (defined($BT_water_depth)) {
- $min_hab = $BT_water_depth - $LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH};
- printf(STDERR "\n\twater depth = %d(+-%.1f)m",$BT_water_depth,$sig_BT_water_depth);
+ print(STDERR "\n\tin acoustic backscatter profiles...") if ($opt_d);
+ ($water_depth,$sig_water_depth) =
+ find_backscatter_seabed($LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH});
+
+ $min_hab = $water_depth - $LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH};
+ printf(STDERR "\n\twater depth = %d(+-%.1f)m",$water_depth,$sig_water_depth);
printf(STDERR "\n\tclosest approach = %dmab",$min_hab);
-# $water_depth = $BT_water_depth; # assume BT data are better
-# $sig_water_depth = $sig_BT_water_depth; # (at least they are higher vertical resolution)
+
+ print(STDERR "\n\tin RDI BT data...") if ($opt_d);
+ ($BT_water_depth,$sig_BT_water_depth) =
+ find_seabed(\%LADCP,$LADCP_bottom,$LADCP{BEAM_COORDINATES});
+
+ if (defined($BT_water_depth)) {
+ $min_hab = $BT_water_depth - $LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH};
+ printf(STDERR "\n\twater depth = %d(+-%.1f)m",$BT_water_depth,$sig_BT_water_depth);
+ printf(STDERR "\n\tclosest approach = %dmab",$min_hab);
+# $water_depth = $BT_water_depth; # assume BT data are better
+# $sig_water_depth = $sig_BT_water_depth; # (at least they are higher vertical resolution)
+ }
+
+ unless (defined($water_depth)) {
+ print(STDERR "\n\tno seabed found\n");
+ print(STDERR "\n\tunknown water depth => PPI editing disabled\n")
+ if ($opt_d);
+ $clip_margin = 0;
+ }
+
+ print(STDERR "\n");
+
+} else { # UPLOOKER
+ $water_depth = $sig_water_depth = $min_hab = nan;
+ print(STDERR "Uplooker data => PPI editing disabled\n")
+ if ($opt_d);
}
-unless (defined($water_depth)) {
- print(STDERR "\n\tno seabed found\n");
- print(STDERR "\n\tunknown water depth => PPI editing disabled\n")
- if ($opt_d);
- $clip_margin = 0;
-}
-
-print(STDERR "\n");
-
#----------------------------------------------------------------------
# Step 8: Edit Data
#----------------------------------------------------------------------
-print(STDERR "Calculating shear profiles...\n");
+print(STDERR "Calculating shear profiles...");
$LADCP_start = 1 if ($LADCP_start == 0); # ensure that there is previous ensemble
-print(STDERR "\tdowncast...");
+print(STDERR "\n\tdowncast...") if ($opt_d);
edit_velocity($LADCP_start,$LADCP_bottom); # downcast
calc_shear($LADCP_start,$LADCP_bottom,$SHEAR_PREGRID_DZ,0); # pre-grid shear @SHEAR_PREGRID_DZm resolution
calc_shear($LADCP_start,$LADCP_bottom,$GRID_DZ,1); # calculate final gridded shear profile
@@ -393,7 +414,7 @@
@dc_vsh_mu = @vsh_mu; @dc_vsh_sig = @vsh_sig;
@dc_wsh_mu = @wsh_mu; @dc_wsh_sig = @wsh_sig;
-print(STDERR "\n\tupcast...");
+print(STDERR "\n\tupcast...") if ($opt_d);
edit_velocity($LADCP_end,$LADCP_bottom); # upcast
calc_shear($LADCP_end,$LADCP_bottom,$SHEAR_PREGRID_DZ,0);
calc_shear($LADCP_end,$LADCP_bottom,$GRID_DZ,1);
@@ -403,7 +424,7 @@
@uc_vsh_mu = @vsh_mu; @uc_vsh_sig = @vsh_sig;
@uc_wsh_mu = @wsh_mu; @uc_wsh_sig = @wsh_sig;
-print(STDERR "\n\tcombined...");
+print(STDERR "\n\tcombined...") if ($opt_d);
for (my($gi)=0; $gi<@dc_ush_mu; $gi++) {
if ($dc_sh_n[$gi]>0 && $uc_sh_n[$gi]>0) {
$sh_n[$gi] = $dc_sh_n[$gi] + $uc_sh_n[$gi];
@@ -434,12 +455,14 @@
# Step 9: Get bottom track profile
#----------------------------------------------------------------------
-print(STDERR "Getting BT profile...");
-getBTprof($LADCP_start,$LADCP_end);
-print(STDERR "\n");
+if ($LADCP{ENSEMBLE}[$LADCP_start]->{XDUCER_FACING_DOWN}) {
+ print(STDERR "Getting BT profile...");
+ getBTprof($LADCP_start,$LADCP_end);
+ print(STDERR "\n");
+}
#----------------------------------------------------------------------
-# Step 10: Write Output
+# Step 10: Write Output Files
#----------------------------------------------------------------------
print(STDERR "Writing shear profiles...");
@@ -460,7 +483,15 @@
'water_depth',round($water_depth),'water_depth.sig',round($sig_water_depth),
'min_hab',round($min_hab),
'clip_margin',$clip_margin,'first_clip_bin',$first_clip_bin,
- 'Svbin_start',$Svbin_start,'Svbin_end',$Svbin_end);
+ 'Svbin_start',$Svbin_start,'Svbin_end',$Svbin_end,
+ 'BT_bin_start',$BT_bin_start,'BT_bin_search_above',$BT_bin_search_above,
+ 'BT_max_bin_spread',$BT_max_bin_spread,'BT_max_w_difference',$BT_max_w_difference,
+);
+if (defined($BT_min_depth)) {
+ &antsAddParams('BT_min_depth',$BT_min_depth,'BT_max_depth',$BT_max_depth);
+} else {
+ &antsAddParams('BT_max_depth_error',$BT_max_depth_error);
+}
for (my($gi)=0; $gi<@ush_mu; $gi++) {
&antsOut(depthOfGI($gi), # depth in center of bin
@@ -483,33 +514,49 @@
print(STDERR "\n");
-#----------------------------------------------------------------------
+#---------------------------------------
+# Acoustic backscatter depth-time-series
+#---------------------------------------
if (defined($opt_a)) {
- print(STDERR "Writing per-bin acoustic backscatter profiles...");
+ print(STDERR "Writing acoustic backscatter depth-time-series to <$opt_a>...");
+
+
+ @antsNewLayout = ('ens','elapsed','CTD_depth','depth','bin','downcast',
+ 'Sv_beam1','Sv_beam2','Sv_beam3','Sv_beam4','Sv.median');
+ &antsOut('EOF');
+ close(STDOUT);
+ open(STDOUT,">$opt_a") || croak("$opt_a: $!\n");
- for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
- my($fn) = sprintf("bin%02d.Sv",$bin);
- print(STDERR " $fn");
-
- @antsNewLayout = ('depth','Sv');
- &antsOut('EOF');
- $antsCurParams = $commonParams;
- close(STDOUT);
- open(STDOUT,">$fn") || croak("$fn: $!\n");
-
- for (my($gi)=0; $gi<@sSv; $gi++) {
- &antsOut(depthOfGI($gi),
- $nSv[$gi][$bin] ? $sSv[$gi][$bin]/ $nSv[$gi][$bin] : nan);
- }
+ $antsCurParams = $commonParams;
+ &antsAddParams('min_elapsed',$LADCP{ENSEMBLE}[$LADCP_start]->{ELAPSED_TIME},
+ 'max_elapsed',$LADCP{ENSEMBLE}[$LADCP_end]->{ELAPSED_TIME},
+ 'min_depth',$LADCP{ENSEMBLE}[$LADCP_top]->{XDUCER_FACING_UP} ?
+ &depthOfBin($LADCP_top,$LADCP{N_BINS}-1) : $LADCP{ENSEMBLE}[$LADCP_top]->{DEPTH},
+ 'max_depth',$LADCP{ENSEMBLE}[$LADCP_bottom]->{XDUCER_FACING_UP} ?
+ $LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH} : &depthOfBin($LADCP_bottom,$LADCP{N_BINS}-1)
+ );
+
+ for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
+ for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
+ &antsOut($LADCP{ENSEMBLE}[$ens]->{NUMBER},
+ $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME},
+ $LADCP{ENSEMBLE}[$ens]->{DEPTH},
+ &depthOfBin($ens,$bin),$bin,
+ ($ens <= $LADCP_bottom) ? 1 : 0,
+ @{$LADCP{ENSEMBLE}[$ens]->{SV}[$bin]},
+ median(@{$LADCP{ENSEMBLE}[$ens]->{SV}[$bin]})
+ );
+ }
}
+
print(STDERR "\n");
}
#----------------------------------------------------------------------
if (defined($opt_t)) {
- print(STDERR "Writing time series to $opt_t...");
+ print(STDERR "Writing time series to <$opt_t>...");
@antsNewLayout = ('ens','elapsed','depth','CTD_w','LADCP_w');
&antsOut('EOF');
@@ -530,7 +577,7 @@
#----------------------------------------------------------------------
if (defined($opt_b)) {
- print(STDERR "Writing bottom-track data to $opt_b...");
+ print(STDERR "Writing bottom-track data to <$opt_b>...");
@antsNewLayout = ('depth','u','v','w','u.sig','v.sig','w.sig','ndata');
&antsOut('EOF');
@@ -550,7 +597,7 @@
#----------------------------------------------------------------------
if (defined($opt_f)) {
- print(STDERR "Writing data flags to $opt_f...");
+ print(STDERR "Writing data flags to <$opt_f>...");
@antsNewLayout = ('ens');
for (my($i)=1; $i<=$LADCP{N_BINS}; $i++) {
--- a/LADCPproc.BT
+++ b/LADCPproc.BT
@@ -1,52 +1,65 @@
#======================================================================
# L A D C P P R O C . B T
# doc: Wed Oct 20 21:05:37 2010
-# dlm: Mon Jan 10 12:07:17 2011
+# dlm: Fri Jul 8 03:24:56 2011
# (c) 2010 A.M. Thurnherr
-# uE-Info: 11 28 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 148 69 NIL 0 0 72 10 2 4 NIL ofnI
#======================================================================
# HISTORY:
# Oct 20, 2010: - created
# Jan 10, 2010: - -o => -k
+# Jul 7, 2010: - added $DEBUG
+# - added BTrangeFlag
+# - added $BT processing parameters
+# - changed from echo amplitude to Sv
my($BEAM1) = 0;
my($BEAM2) = 1;
my($BEAM3) = 2;
my($BEAM4) = 3;
-my($nBTfound,$nBTdepthFlag,$nBTvalidVelFlag,$nBTwFlag) = (0,0,0,0);
+my($nBTfound,$nBTrangeFlag,$nBTdepthFlag,$nBTvalidVelFlag,$nBTwFlag) = (0,0,0,0,0);
+
+my($DEBUG) = 0;
sub binBTprof($)
{
my($ens) = @_;
- my(@ea_max) = (0,0,0,0); my(@ea_max_bin) = (nan,nan,nan,nan);
- for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
- $ea_max[$BEAM1] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][$BEAM1],
- $ea_max_bin[$BEAM1] = $bin
- if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][$BEAM1] > $ea_max[$BEAM1]);
- $ea_max[$BEAM2] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][$BEAM2],
- $ea_max_bin[$BEAM2] = $bin
- if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][$BEAM2] > $ea_max[$BEAM2]);
- $ea_max[$BEAM3] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][$BEAM3],
- $ea_max_bin[$BEAM3] = $bin
- if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][$BEAM3] > $ea_max[$BEAM3]);
- $ea_max[$BEAM4] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][$BEAM4],
- $ea_max_bin[$BEAM4] = $bin
- if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][$BEAM4] > $ea_max[$BEAM4]);
+ my(@Sv_max) = (-9e99,-9e99,-9e99,-9e99); my(@Sv_max_bin) = (nan,nan,nan,nan);
+ for (my($bin)=$BT_bin_start-1; $bin<$LADCP{N_BINS}; $bin++) {
+ if (defined($BT_min_depth)) { # manually supplied bottom depth range
+ my($dob) = &depthOfBin($ens,$bin);
+ next unless ($dob >= $BT_min_depth && $dob <= $BT_max_depth);
+ }
+ $Sv_max[$BEAM1] = $LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM1],
+ $Sv_max_bin[$BEAM1] = $bin
+ if ($LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM1] > $Sv_max[$BEAM1]);
+ $Sv_max[$BEAM2] = $LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM2],
+ $Sv_max_bin[$BEAM2] = $bin
+ if ($LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM2] > $Sv_max[$BEAM2]);
+ $Sv_max[$BEAM3] = $LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM3],
+ $Sv_max_bin[$BEAM3] = $bin
+ if ($LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM3] > $Sv_max[$BEAM3]);
+ $Sv_max[$BEAM4] = $LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM4],
+ $Sv_max_bin[$BEAM4] = $bin
+ if ($LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM4] > $Sv_max[$BEAM4]);
}
-# print(STDERR "@ea_max | @ea_max_bin\n");
-
- return unless (max(@ea_max_bin)-min(@ea_max_bin) <= 3); # inconsistent range (&, impliclity, large tilt)
- # SHOULD CONSIDER RELAXING LATER
+ print(STDERR "@Sv_max | @Sv_max_bin\n") if ($DEBUG);
$nBTfound++;
- my($range_bin) = round(avg(@ea_max_bin));
-# printf(STDERR "water_depth = $water_depth; BT peak depth = %d\n",depthOfBin($ens,$range_bin));
+ $nBTrangeFlag++,return # inconsistent range (&, impliclity, large tilt)
+ unless (max(@Sv_max_bin)-min(@Sv_max_bin) <= $BT_max_bin_spread);
+
+ my($range_bin) = round(avg(@Sv_max_bin));
+ printf(STDERR "water_depth = $water_depth; BT peak depth = %d in bin $range_bin\n",depthOfBin($ens,$range_bin))
+ if ($DEBUG);
+
$nBTdepthFlag++,return # BT range inconsistent with water depth
- unless (abs($water_depth-depthOfBin($ens,$range_bin)) < 20+$sig_water_depth);
+ unless defined($BT_min_depth) ||
+ (abs($water_depth-depthOfBin($ens,$range_bin)) < $sig_water_depth + $BT_max_depth_error);
# try bin of max plus one above and below
# this does not really work because, often, only one of the bins has valid velocities
@@ -54,7 +67,8 @@
my($w2) = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin][$W];
my($w3) = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+1][$W];
-# printf(STDERR "w123 = %.1f,%.1f,%.1f\n",$w1,$w2,$w3);
+ printf(STDERR "w123 = %.1f,%.1f,%.1f\n",$w1,$w2,$w3)
+ if ($DEBUG);
$w1 = 9e99 unless numberp($w1); # no valid velocities
$w2 = 9e99 unless numberp($w1);
@@ -78,12 +92,13 @@
$CTD_w = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin][$W];
}
- $nBTwFlag++,return if (abs($CTD_w-$LADCP{ENSEMBLE}[$ens]->{W}) > 0.03); # velocity error is too great
+ $nBTwFlag++,return if (abs($CTD_w-$LADCP{ENSEMBLE}[$ens]->{W}) > $BT_max_w_difference);
-# printf(STDERR "good BT [%5.2f %5.2f %5.2f] found at ens $ens\n",$CTD_u,$CTD_v,$CTD_w);
+ printf(STDERR "good BT [%5.2f %5.2f %5.2f] found at ens $ens\n",$CTD_u,$CTD_v,$CTD_w)
+ if ($DEBUG);
if ($opt_k) {
- for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
+ for (my($bin)=$BT_bin_start-1; $bin<$LADCP{N_BINS}; $bin++) {
next if ($edit_flags[$ens][$bin]);
printf(BTF "%d %d %d %f %f %f %f %f %f %f %f %f %f %f\n",
$LADCP{ENSEMBLE}[$ens]->{NUMBER},
@@ -100,7 +115,7 @@
print(BTF "nan nan nan nan\n");
}
- for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
+ 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;
push(@{$BTu_vals[$gi]},$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$U]-$CTD_u);
@@ -120,15 +135,19 @@
}
for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
- next unless ($water_depth-$LADCP{ENSEMBLE}[$ens]->{DEPTH} < 300);
+ next unless ($water_depth-$LADCP{ENSEMBLE}[$ens]->{DEPTH} < $BT_begin_search_above);
binBTprof($ens);
}
if ($opt_d) {
print(STDERR "\n\t$nBTfound BT ensembles found\n");
+ print(STDERR "\t\t$nBTrangeFlag flagged bad because of inconsistent range to seabed\n");
print(STDERR "\t\t$nBTdepthFlag flagged bad because of wrong bottom depth\n");
- print(STDERR "\t\t$nBTvalidVelFlag flagged bad because of no valid velocities\n");
- print(STDERR "\t\t$nBTwFlag flagged bad because of incorrect vertical velocity");
+ print(STDERR "\t\t$nBTvalidVelFlag flagged bad because of lack of valid velocities\n");
+ print(STDERR "\t\t$nBTwFlag flagged bad because of incorrect vertical velocities");
+ printf(STDERR "\n\t=> %d velocities from %d BT ensembles used",
+ scalar(@BTu_vals),
+ $nBTfound-$nBTrangeFlag-$nBTdepthFlag-$nBTvalidVelFlag-$nBTwFlag);
}
@BTu = @BTv = @BTw = ();
--- 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: Wed Jun 15 15:31:43 2011
+# dlm: Thu Jul 7 08:55:45 2011
# (c) 2010 A.M. Thurnherr
-# uE-Info: 67 0 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 96 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -12,6 +12,13 @@
# when run with uplooker data
# Jun 15, 2011: - added calculation of backscatter profiles from
# subset of bins
+# Jul 7, 2011: - use $BEAM? vars to clarify code
+# - save SV values to use in BT code
+
+my($BEAM1) = 0;
+my($BEAM2) = 1;
+my($BEAM3) = 2;
+my($BEAM4) = 3;
#----------------------------------------------------------------------
# Volume Scattering Coefficient, following Deines (IEEE 1999)
@@ -53,14 +60,14 @@
my(@Er) = (1e99,1e99,1e99,1e99); # echo intensity reference level
for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
- $Er[0] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][0]
- if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][0] < $Er[0]);
- $Er[1] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][1]
- if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][1] < $Er[1]);
- $Er[2] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][2]
- if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][2] < $Er[2]);
- $Er[3] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][3]
- if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][3] < $Er[3]);
+ $Er[$BEAM1] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][$BEAM1]
+ if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][$BEAM1] < $Er[$BEAM1]);
+ $Er[$BEAM2] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][$BEAM2]
+ if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][$BEAM2] < $Er[$BEAM2]);
+ $Er[$BEAM3] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][$BEAM3]
+ if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][$BEAM3] < $Er[$BEAM3]);
+ $Er[$BEAM4] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][$BEAM4]
+ if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][$BEAM4] < $Er[$BEAM4]);
}
print(STDERR "\n\t\@per-beam noise levels = @Er") if ($opt_d);
@@ -69,22 +76,23 @@
my($gi) = int(&depthOfBin($ens,$bin) / $GRID_DZ);
next if ($gi < 0);
my($range_to_bin) = &dzToBin($ens,$bin) / cos(rad($LADCP{BEAM_ANGLE}));
- my($Sv) = Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},
- $LADCP{TRANSMITTED_PULSE_LENGTH},
- $Er[0],$range_to_bin,
- $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][0])/4 +
- Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},
- $LADCP{TRANSMITTED_PULSE_LENGTH},
- $Er[1],$range_to_bin,
- $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][1])/4 +
- Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},
- $LADCP{TRANSMITTED_PULSE_LENGTH},
- $Er[2],$range_to_bin,
- $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][2])/4 +
- Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},
- $LADCP{TRANSMITTED_PULSE_LENGTH},
- $Er[3],$range_to_bin,
- $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][3])/4;
+ $LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM1] = Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},
+ $LADCP{TRANSMITTED_PULSE_LENGTH},
+ $Er[$BEAM1],$range_to_bin,
+ $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][$BEAM1]);
+ $LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM2] = Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},
+ $LADCP{TRANSMITTED_PULSE_LENGTH},
+ $Er[$BEAM2],$range_to_bin,
+ $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][$BEAM2]);
+ $LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM3] = Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},
+ $LADCP{TRANSMITTED_PULSE_LENGTH},
+ $Er[$BEAM3],$range_to_bin,
+ $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][$BEAM3]);
+ $LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM4] = Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},
+ $LADCP{TRANSMITTED_PULSE_LENGTH},
+ $Er[$BEAM4],$range_to_bin,
+ $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][$BEAM4]);
+ my($Sv) = avg(@{$LADCP{ENSEMBLE}[$ens]->{SV}[$bin]});
$sSv[$gi][$bin] += $Sv;
$nSv[$gi][$bin]++;
--- a/LADCPproc.bestLag
+++ b/LADCPproc.bestLag
@@ -1,9 +1,9 @@
#======================================================================
# L A D C P P R O C . B E S T L A G
# doc: Tue Sep 28 21:58:48 2010
-# dlm: Wed Jan 5 19:44:09 2011
+# dlm: Fri Jul 8 02:54:29 2011
# (c) 2010 A.M. Thurnherr
-# uE-Info: 67 28 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 55 101 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# TODO:
@@ -16,6 +16,8 @@
# Jan 5, 2011: - changed first guess from 80% down to 10% down
# - added LADCP time lag to %PARAMs
# - added support of -i
+# Jul 7, 2011: - added code to remove window-mean of w before lagging to
+# make it work in regions of crazy ocean w (IWISE 16007)
sub interp_LADCP_w($$)
{
@@ -38,15 +40,26 @@
my($best);
my($bestmad) = 9e99; # mean absolute deviation
for (my($Llag)=-int($opt_w/2); $Llag<int($opt_w/2); $Llag++) {
- my($sad) = my($nad) = 0;
+ my($mCw,$mLw,$nw) = (0,0,0); # first calc means
for (my($Ci)=0; $Ci<$opt_w; $Ci++) {
my($Li) = $Ci + $Llag;
next if ($Li<0 || $Li>=$opt_w);
next unless numberp($CTD{w}[$ws+$Ci]) && numberp($LADCP_w[$ws+$Li]);
- $sad += abs($CTD{w}[$ws+$Ci] - $LADCP_w[$ws+$Li]);
+ $mCw += $CTD{w}[$ws+$Ci];
+ $mLw += $LADCP_w[$ws+$Li];
+ $nw++;
+ }
+ next unless ($nw > 0);
+ $mCw /= $nw; $mLw /= $nw;
+
+ my($sad) = my($nad) = 0; # calc mad with means removed
+ for (my($Ci)=0; $Ci<$opt_w; $Ci++) {
+ my($Li) = $Ci + $Llag;
+ next if ($Li<0 || $Li>=$opt_w);
+ next unless numberp($CTD{w}[$ws+$Ci]) && numberp($LADCP_w[$ws+$Li]);
+ $sad += abs($CTD{w}[$ws+$Ci]-$mCw - ($LADCP_w[$ws+$Li]-$mLw));
$nad++;
}
- next unless ($nad > 0);
if ($sad/$nad < $bestmad) {
$best = $Llag;
$bestmad = $sad/$nad;
--- a/LADCPproc.defaults
+++ b/LADCPproc.defaults
@@ -1,9 +1,9 @@
#======================================================================
# L A D C P P R O C . D E F A U L T S
# doc: Fri Sep 17 09:44:21 2010
-# dlm: Wed Jun 15 15:15:41 2011
+# dlm: Fri Jul 8 00:06:07 2011
# (c) 2010 A.M. Thurnherr
-# uE-Info: 23 48 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 91 16 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# default parameters for [mkShearProf]
@@ -15,12 +15,14 @@
# - the default version in the ANTS bin dir is always loaded
# - if there is a version in the current processing directory it is loaded
# afterwards
-# - for additional notes, see [mkShearProf]
+# - per-station parameters can be chosed based on $LADCP_file and $CTD_file
+# - for additional notes, see [LADCPproc]
# HISTORY:
# Sep 17, 2010: - created
# Dec 9, 2010: - added doc for ASCII CTD file support
# Jun 15, 2011: - added Svbin_start, Svbin_end
+# Jul 7, 2011: - added BT processing parameters
#----------------------------------------------------------------------
# ASCII CTD file support
@@ -62,6 +64,48 @@
$Svbin_end = 9;
#----------------------------------------------------------------------
+# Bottom Track Parameters
+#----------------------------------------------------------------------
+
+# First bin to consider when looking for seabed to calculated post-
+# processed BT. NB: For consistency with EF's code, bin 1 is the
+# first bin.
+
+$BT_bin_start = 1;
+
+# CTD depth above seabed where BT search begins
+
+$BT_begin_search_above = 300;
+
+# Maximum allowed spread of bin number where max echo is found. Large values
+# imply sloping seabed and/or large instrument tilt. The default value of
+# 3 is probably too tight. The range spread calculation should be modified
+# to take instrument tilt into consideration.
+
+$BT_max_bin_spread = 3;
+
+# For tricky BT cases, the acoustic backscatter from the seabed can be plotted
+# and the bottom range set manually. This has been found to be very powerful
+# in case of 2011_IWISE yoyo profile 160, where the bottom was around 150m
+# from the lower turning points and where the method of Deines [1999] clearly
+# does not work.
+
+#$BT_min_depth = xxx;
+#$BT_max_depht = xxx;
+
+# Maximum difference between water depth and average distance of echo max.,
+# if $BT_min_depth and $BT_max_depth are not set. The stddev of the detected
+# water depth is added To this number.
+
+$BT_max_depth_error = 20;
+
+# Maximum allowed difference between reference-layer w and BT w. Note
+# that this must be relaxed in regions of significant near-bottom
+# vertical velocity, e.g. due to cross-slope currents.
+
+$BT_max_w_difference = 0.03;
+
+#----------------------------------------------------------------------
# Shear Processing Parameters
#----------------------------------------------------------------------
deleted file mode 100755
--- a/splitCNV
+++ /dev/null
@@ -1,48 +0,0 @@
-#!/usr/bin/perl
-#======================================================================
-# S P L I T C N V
-# doc: Mon Oct 25 22:24:21 2010
-# dlm: Fri Oct 29 20:04:42 2010
-# (c) 2010 A.M. Thurnherr
-# uE-Info: 16 29 NIL 0 0 72 2 2 4 NIL ofnI
-#======================================================================
-
-# split SeaBird CNV file using scan values found with yoyo
-# e.g. splitCNV P403_015_1sec.cnv `importCNV -s P403_015_1sec.cnv | yoyo -QFscan -ut`
-
-# HISTORY:
-# Oct 25, 2010: - created
-# Oct 26, 2010: - BUG: *.00 file always only contained the 1st record
-# Oct 29, 2010: - cosmetics
-
-# scan number must be the first field of each record
-
-die("Usage: $0 <CNV-file> <val> [...]\n")
- unless (@ARGV >= 2 && -f $ARGV[0]);
-
-chomp($basename = `basename $ARGV[0] .cnv`);
-
-$state = 0; # init finite state machine
-$next = 0; # next extension number
-$trg = 2; # NB: yoyo -ut includes first & last scans
-
-open(IN,$ARGV[0]);
-while (<IN>) {
- if ($state == 0) { # state 0: reading header
- push(@HDR,$_);
- $state = 1 if /^\*END\*/;
- next;
- }
- if ($state == 1) { # state 1: begin new file
- close(OUT);
- printf(STDERR "writing $basename.%02d\n",$next);
- open(OUT,sprintf(">$basename.%02d",$next++));
- foreach my $h (@HDR) { print(OUT $h); }
- $state = 2;
- }
- if ($state == 2) { # state 2: copy data until target scan
- print(OUT);
- my(@f) = split;
- $trg++,$state=1 if ($f[0] == $ARGV[$trg]);
- }
-}
new file mode 100755
--- /dev/null
+++ b/splitCNV.olderVersion
@@ -0,0 +1,48 @@
+#!/usr/bin/perl
+#======================================================================
+# S P L I T C N V
+# doc: Mon Oct 25 22:24:21 2010
+# dlm: Fri Oct 29 20:04:42 2010
+# (c) 2010 A.M. Thurnherr
+# uE-Info: 16 29 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# split SeaBird CNV file using scan values found with yoyo
+# e.g. splitCNV P403_015_1sec.cnv `importCNV -s P403_015_1sec.cnv | yoyo -QFscan -ut`
+
+# HISTORY:
+# Oct 25, 2010: - created
+# Oct 26, 2010: - BUG: *.00 file always only contained the 1st record
+# Oct 29, 2010: - cosmetics
+
+# scan number must be the first field of each record
+
+die("Usage: $0 <CNV-file> <val> [...]\n")
+ unless (@ARGV >= 2 && -f $ARGV[0]);
+
+chomp($basename = `basename $ARGV[0] .cnv`);
+
+$state = 0; # init finite state machine
+$next = 0; # next extension number
+$trg = 2; # NB: yoyo -ut includes first & last scans
+
+open(IN,$ARGV[0]);
+while (<IN>) {
+ if ($state == 0) { # state 0: reading header
+ push(@HDR,$_);
+ $state = 1 if /^\*END\*/;
+ next;
+ }
+ if ($state == 1) { # state 1: begin new file
+ close(OUT);
+ printf(STDERR "writing $basename.%02d\n",$next);
+ open(OUT,sprintf(">$basename.%02d",$next++));
+ foreach my $h (@HDR) { print(OUT $h); }
+ $state = 2;
+ }
+ if ($state == 2) { # state 2: copy data until target scan
+ print(OUT);
+ my(@f) = split;
+ $trg++,$state=1 if ($f[0] == $ARGV[$trg]);
+ }
+}