#!/usr/bin/perl
#======================================================================
# 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
# (c) 2010 A.M. Thurnherr
# uE-Info: 103 69 NIL 0 0 72 11 2 4 NIL ofnI
#======================================================================
# NOTES:
# - the shear-method editing in this code is based on Eric Firing's merge.c
# - as described in [LADCPproc.backscatter], there are three different codes
# for correcting echo amplitudes for attenuation loss & beam spreading
# - comments starting with ## are taken verbatim from the original
# implementations
# - the basic idea of the time lagging implemented in this code is similar
# to the one implemented in Martin Visbeck's bestlag.m
# - for SeaBird files:
# - CTD elapsed time is estimated from recno * CTD{sampint}
# - first elapsed time is added on output
# - output time:
# - is from CTD to allow yoyo processing without loss of time information
# - to obtain LADCP time, add %LADCP_time_lag to CTD time
# - %LADCP_time_lag can be used for -l
# - -i is different: it should be set to the number that is added to
# LADCP_elapsed to get the CTD time, i.e.
# -%LADCP_time_lag
# - CTD{elapsed} is undefined for records before instrument is in the water
# - ITS-90 temp field in degC required
# - salin field prequired
# - pressure field in dbar required
# HISTORY:
# Sep 16, 2010: - incepted
# Oct 13, 2010: - first working version
# Oct 14, 2010: - renamed from LADCPshear
# Oct 19, 2010: - added -a)coustic backscatter profiles
# Oct 20, 2010: - added -2)dary CTD sensors
# Oct 23, 2010: - added magnetic-declination correction
# Oct 26, 2010: - added tilt calculation
# Dec 9, 2010: - added support for ASCII CTD files
# Dec 10, 2010: - change -w) default to 120s
# - changed nshear output to 0 from nan when there are no samples
# Dec 27, 2010: - changed sign of -l to accept lag output from [LADCP_w]
# Jan 10, 2011: - -o => -k added new -o
# - added code to fill CTD sound vel gaps
# 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
# Jul 11, 2011: - changed default output to .tds and added -p)rofile option
# Jul 14, 2011: - added -u)se D[eines1999][V[isbeck2004]]
# Jul 15, 2011: - changed output bin# to 1-based in -a output
# - added T[hurnherr11] -u)se option
# - added $CTD{first_elapsed}
# Jul 27, 2011: - replaced ndata by nsamp
# Feb 5, 2012: - added profile max depth consistency check
# Feb 19, 2912: - added elapsed time to shear profile output
# - replaced "nshear" output field by "nsamp"
# - BUG: bottom of profiles was incorrect when dc max depth > uc max depth
# - BUG: profile depth consistency check did not work for partial-depth yoyo casts
# Apr 10, 2012: - changed default backscatter correction to Deines (1999)
# - improved and relaxed depth consistency check
# Apr 11, 2012: - BUG: double comma that did not affect the output
# - BUG: code had assumed, but not ensured, that first CTD scan is valid
# Apr 13, 2012: - removed -s argument from dependencies
# Apr 17, 2012: - BAD BUG: magdec code call was bad and did not return correct value. ever.
# Apr 18, 2012: - replaced Sv.n by Sv.nsamp
# May 16, 2012: - adapted to ANTSlib V5.0
# - added support for -r)DI BT data
# - added $LADCP{MEAN_DT} from vertical-velocity code
# May 17, 2012: - significantly relaxed profile-range consistency check (from 20% to 50%)
# which is required for Dan Torres' shallow casts
# - added version info
# May 18, 2012: - reduced "implausibly short cast" threshold to 5 min because of Dan's shallow casts
# Jun 13, 2012: - added CTD_depth to tds output
# Oct 28, 2012: - added -z
# Mar 16, 2013: - BUG: usage message had a -p)PI flag
# Jun 5, 2013: - added $bad_beam support
# Jun 25, 2013: - added %PARAMS used for spectral correction
# - adapted to new ::-PARAM convention
# Sep 25, 2013: - BUG: %PARAM magnetic_declination did not have LADCPproc:: prefix
# - added CTD lat/lon info to most output files (but not BT)
# - BUG: moved %water_depth to common %PARAMs
# Feb 22, 2014: - modified diagnostics output
# Mar 19, 2014: - moved code to set LADCP_time_lag %PARAM into main prog so it is
# set, even when -l is used
# - added pitch, roll, hdg to -t output
# Mar 20, 2014: - BUG: wrong number of samples were recorded when upcast had no
# valid data whatsoever
# - added support for $LADCP_max_gap
# Mar 21, 2014: - added [LADCPproc.utils]
# - added $ignore_tilt
# Mar 24, 2014: - added $pitch_offset, $roll_offset
# Mar 27, 2014: - BUG: Sv output did not reflect true bin depths
# Apr 25, 2014: - added LADCP_errvel to -t output
# 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
($ANTS) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
($PERL_TOOLS) = (`which mkProfile` =~ m{^(.*)/[^/]*$});
($LADCPPROC) = ($0 =~ m{^(.*)/[^/]*$});
require "$ANTS/ants.pl";
require "$ANTS/libEOS83.pl";
require "$ANTS/libstats.pl";
require "$LADCPPROC/LADCPproc.version";
require "$LADCPPROC/LADCPproc.utils";
require "$LADCPPROC/LADCPproc.loadCTD";
require "$LADCPPROC/LADCPproc.bestLag";
require "$LADCPPROC/LADCPproc.BT";
require "$LADCPPROC/LADCPproc.backscatter";
require "$LADCPPROC/LADCPproc.shearmethod";
require "$PERL_TOOLS/RDI_BB_Read.pl";
require "$PERL_TOOLS/RDI_Coords.pl";
require "$PERL_TOOLS/RDI_Utils.pl";
$ANTSLIBS = $LADCPPROC; # for -L libraries
$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,
'[use -2)dary CTD sensor pair]',
'[require -4)-beam LADCP solutions]',
'[use -r)DI bottom-track data]',
'[-s)etup <file>] [-g)ps <lat,lon>]',
'[-c)ompass-corr <offset,cos-fac,sin-fac>]',
'[-o)utput grid <resolution[5m]>]',
'[-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]',
'[-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]',
'<RDI file> <SeaBird file>');
$RDI_Coords::minValidVels = 4 if ($opt_4);
&antsFloatOpt($opt_l);
&antsCardOpt(\$opt_w,120);
# old default of -w 30 does not work if there are significant ambiguity-velocity
# problems, as is the case, e.g., with 2010_DIMES_US2 station 142
# old default of -w 60 did not work for DIMES_UK2 station 4 (DL), possibly again
# related to ambiguity velocity
&antsCardOpt(\$opt_n,20);
#&antsFileOpt($opt_s); # DON'T, AS THIS WILL ADD AN UNWANTED DEPENDCY
&antsFloatOpt($opt_i);
&antsCardOpt(\$opt_o,5);
if (defined($opt_u)) {
croak("$0: cannot decode -u $opt_u\n")
unless ($opt_u =~ /^[dDvVtTn]/);
} else {
$opt_u = 'Deines99';
}
if (defined($opt_g)) {
($CTD{stn_lat},$CTD{stn_lon}) = split(',',$opt_g);
croak("$0: cannot decode -g $opt_g\n")
unless numberp($CTD{stn_lat}) && numberp($CTD{stn_lon});
}
if (defined($opt_c)) {
($CC_offset,$CC_cos_fac,$CC_sin_fac) = split(',',$opt_c);
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();
&antsAddParams('LADCPproc::LADCP_file',$LADCP_file,'LADCPproc::CTD_file',$CTD_file);
&antsActivateOut();
#----------------------------------------------------------------------
# Step 1: Read LADCP Data
#----------------------------------------------------------------------
print(STDERR "Reading LADCP data ($LADCP_file)...");
readData($LADCP_file,\%LADCP);
printf(STDERR "\n\t%d ensembles",scalar(@{$LADCP{ENSEMBLE}})) if ($opt_d);
print(STDERR "\n");
&antsAddParams('LADCPproc::bin_length',$LADCP{BIN_LENGTH});
#----------------------------------------------------------------------
# Step 2: Set Processing Parameters
#----------------------------------------------------------------------
print(STDERR "Setting processing parameters...\n");
printf(STDERR "\tloading $LADCPPROC/LADCPproc.defaults...\n");
require "$LADCPPROC/LADCPproc.defaults";
if (defined($opt_s)) {
print(STDERR "\tloading $opt_s...\n");
require $opt_s;
}
if ($LADCP{BLANKING_DISTANCE} == 0) {
print(STDERR "\t\tBLANKING_DISTANCE == 0 => excluding all data from bin 1\n")
if ($opt_d);
$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('LADCPproc::instrument_orientation',
$LADCP{ENSEMBLE}[0]->{XDUCER_FACING_UP} ? 'UL' : 'DL');
$SHEAR_PREGRID_DZ = 20;
$GRID_DZ = $opt_o;
&antsAddParams('LADCPproc::vertical_resolution',$GRID_DZ);
#----------------------------------------------------------------------
# Step 3: Read CTD data
#----------------------------------------------------------------------
print(STDERR "Reading CTD data ($CTD_file)...");
readCTD($CTD_file,\%CTD);
printf(STDERR "\n\t%d scans (%.1fs sample interval)",scalar(@{$CTD{press}}),$CTD{sampint}) if ($opt_d);
print(STDERR "\n");
my($year) = substr($LADCP{ENSEMBLE}[0]->{DATE},6,4);
my($month) = substr($LADCP{ENSEMBLE}[0]->{DATE},0,2);
my($day ) = substr($LADCP{ENSEMBLE}[0]->{DATE},3,2);
my($magdec,$maginc,$h_strength,$v_strength) = split('\s+',`magdec $CTD{stn_lon} $CTD{stn_lat} $year $month $day`);
croak("$0: unknown magnetic declination\n")
unless defined($magdec);
&antsAddParams('LADCPproc::magnetic_declination',$magdec);
#----------------------------------------------------------------------
# Step 4: Pre-Process CTD & LADCP Data
#----------------------------------------------------------------------
printf(STDERR "Pre-processing data...");
printf(STDERR "\n\tCTD...") if ($opt_d);
#------------------------
# clean CTD pressure data
#------------------------
my($pSpikes) = 0;
for (my($r)=1; $r<@{$CTD{press}}; $r++) {
$pSpikes++,$CTD{press}[$r]=nan
if (abs($CTD{press}[$r]-$CTD{press}[$r-1])/$CTD{sampint} > 2);
}
print(STDERR "\n\t\t$pSpikes CTD pressure spikes removed")
if ($pSpikes>0 && $opt_d);
#-----------------------------------
# trim "preamble" without valid data
#-----------------------------------
my($r) = 0;
until (numberp($CTD{press}[$r]) && numberp($CTD{temp}[$r]) && numberp($CTD{salin}[$r])) { $r++ }
splice(@{$CTD{press}},0,$r);
splice(@{$CTD{temp}},0,$r);
splice(@{$CTD{salin}},0,$r);
print(STDERR "\n\t\t$r leading CTD scans without valid data removed")
if ($r>0 && $opt_d);
#--------------------------------------------------
# calculate w and find deepest & shallowest records
#--------------------------------------------------
$CTD{maxpress} = -9e99;
$CTD{minpress} = 9e99;
for (my($r)=1; $r<@{$CTD{press}}-1; $r++) {
$CTD{w}[$r] = 0.99*($CTD{press}[$r+1] - $CTD{press}[$r-1]) / (2*$CTD{sampint});
if ($CTD{press}[$r] > $CTD{maxpress}) {
$CTD{maxpress} = $CTD{press}[$r];
$CTD{atbottom} = $r;
}
if ($CTD{press}[$r] < $CTD{minpress}) {
$CTD{minpress} = $CTD{press}[$r];
$CTD{attop} = $r;
}
}
printf(STDERR "\n\t\tmin/max pressure [%d/%ddbar] at scans#%d/%d",
$CTD{minpress},$CTD{maxpress},$CTD{attop},$CTD{atbottom})
if $opt_d;
#----------------------------------------------------------------------
# Step 4b: Pre-Process LADCP Data
#----------------------------------------------------------------------
print(STDERR "\n\tLADCP...") if ($opt_d);
#-------------------------------------------------
# transform to earth coordinates if required
# - discard data from a particular beam if requested
# - save beam_vels for later (e.g. wake) editing
#-------------------------------------------------
$U = 0; # velocity indices
$V = 1;
$W = 2;
$E = 3;
$LADCP{HEADING_BIAS} = -$magdec;
if ($LADCP{BEAM_COORDINATES}) {
if ($opt_d) {
print(STDERR "\n\t\ttransforming beam to Earth coordinates...");
print(STDERR "\n\t\t\tdiscarding data from beam $bad_beam...")
if ($bad_beam);
}
for (my($ens)=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
$LADCP{ENSEMBLE}[$ens]->{TILT} = $ignore_tilt ? 0
: &angle_from_vertical($LADCP{ENSEMBLE}[$ens]->{PITCH}+$pitch_offset,
$LADCP{ENSEMBLE}[$ens]->{ROLL}+$roll_offset);
for (my($bin)=0; $bin<$LADCP{N_BINS}; $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]->{BEAM_VEL}[$bin]} = @{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]};
@{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]} =
velInstrumentToEarth(\%LADCP,$ens,velBeamToInstrument(\%LADCP,@{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]}));
@{$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);
}
@{$LADCP{ENSEMBLE}[$ens]->{BT_VELOCITY}} =
velInstrumentToEarth(\%LADCP,$ens,velBeamToInstrument(\%LADCP,@{$LADCP{ENSEMBLE}[$ens]->{BT_VELOCITY}}));
}
$LADCP{BEAM_COORDINATES} = 0;
$LADCP{EARTH_COORDINATES} = 1;
unless ($opt_4) {
print(STDERR "\n\t\t\t3-beam solutions: $RDI_Coords::threeBeam_1 $RDI_Coords::threeBeam_2 $RDI_Coords::threeBeam_3 $RDI_Coords::threeBeam_4\n")
if ($opt_d);
&antsAddParams('LADCPproc::3_beam_solutions',"$RDI_Coords::threeBeam_1 $RDI_Coords::threeBeam_2 $RDI_Coords::threeBeam_3 $RDI_Coords::threeBeam_4");
}
} elsif ($LADCP{EARTH_COORDINATES}) {
if ($opt_d) {
if ($opt_c) {
printf(STDERR "\n\t\tcalculating tilt and correcting for compass error and magnetic declination of %.1f deg...\n",$magdec);
} else {
printf(STDERR "\n\t\tcalculating tilt and correcting for magnetic declination of %.1f deg...\n",$magdec);
}
}
for (my($ens)=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
$LADCP{ENSEMBLE}[$ens]->{TILT} = $ignore_tilt ? 0
: &angle_from_vertical($LADCP{ENSEMBLE}[$ens]->{PITCH}+$pitch_offset,
$LADCP{ENSEMBLE}[$ens]->{ROLL}+$roll_offset);
my($hdg) = rad($LADCP{ENSEMBLE}[$ens]->{HEADING});
$LADCP{HEADING_BIAS} = ($CC_offset + $CC_cos_fac*cos($hdg) + $CC_sin_fac*sin($hdg)) - $magdec
if ($opt_c);
for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
@{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]} =
velApplyHdgBias(\%LADCP,$ens,@{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]});
}
}
} else {
croak("$0: can only handle beam or earth coordinates\n")
}
#------------------------------------------------------
# construct a depth-vs-time "profile" from the raw data
#------------------------------------------------------
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);
croak("\n$LADCP_file: no good ensembles found\n")
unless defined($LADCP_start);
my($cast_duration) = $LADCP{ENSEMBLE}[$LADCP_end]->{UNIX_TIME} -
$LADCP{ENSEMBLE}[$LADCP_start]->{UNIX_TIME};
croak("$0: implausibly short cast ($cast_duration seconds)\n")
unless ($cast_duration > 300);
$LADCP{MEAN_DT} = $cast_duration / ($LADCP_end-$LADCP_start-1);
if ($opt_d) {
printf(STDERR "\n\t\t\tStart of cast : %s (#%5d) at %6.1fm\n",
$LADCP{ENSEMBLE}[$LADCP_start]->{TIME},
$LADCP{ENSEMBLE}[$LADCP_start]->{NUMBER},
$LADCP{ENSEMBLE}[$LADCP_start]->{DEPTH});
printf(STDERR "\t\t\tBottom of cast : %s (#%5d) at %6.1fm\n",
$LADCP{ENSEMBLE}[$LADCP_bottom]->{TIME},
$LADCP{ENSEMBLE}[$LADCP_bottom]->{NUMBER},
$LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH});
printf(STDERR "\t\t\tEnd of cast : %s (#%5d) at %6.1fm\n",
$LADCP{ENSEMBLE}[$LADCP_end]->{TIME},
$LADCP{ENSEMBLE}[$LADCP_end]->{NUMBER},
$LADCP{ENSEMBLE}[$LADCP_end]->{DEPTH});
printf(STDERR "\t\t\tAvg ping interval: %.1fs",$LADCP{MEAN_DT});
}
print(STDERR "\n");
my($LADCP_prof_range) = $LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH} -
min($LADCP{ENSEMBLE}[$LADCP_start]->{DEPTH},$LADCP{ENSEMBLE}[$LADCP_end]->{DEPTH});
my($CTD_prof_range) = $CTD{maxpress} - $CTD{minpress};
croak(sprintf("$0: LADCP profile depth range (%dm) inconsistent with max CTD pressure range (%ddbar) [-z to disable]\n",
$LADCP_prof_range,$CTD_prof_range))
if (!$opt_z && abs($LADCP_prof_range-0.99*$CTD_prof_range) >
max(0.2*0.99*$CTD_prof_range,abs($LADCP{ENSEMBLE}[$LADCP_end]->{DEPTH}),50));
#----------------------------------------------------------------
# call edit function
# - e.g. for wake editing
# - this is done here, because during step 5 the velocities are
# sound-speed corrected, i.e.
#----------------------------------------------------------------
&edit_LADCP_vels()
if (exists(&edit_LADCP_vels));
#----------------------------------------------------------------------
# Step 5: Add CTD to LADCP Data & correct velocities for sound speed
# - {DEPTH} field is overwritten!
#----------------------------------------------------------------------
print(STDERR "Matching CTD to LADCP time series...");
print(STDERR "\n") if defined($opt_l);
$opt_l = defined($opt_l) ? -$opt_l : &lagLADCP2CTD();
&antsAddParams('LADCPproc::LADCP_time_lag',-$opt_l);
print(STDERR "Associating CTD data with LADCP ensembles...");
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})
while ($r < 0);
printf(STDERR "\n\tCTD data begin with instrument already in water => skipping %ds of LADCP data",
$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$LADCP_start]->{ELAPSED_TIME});
$LADCP_start = $ens;
}
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})
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{stn_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]);
if (numberp($LADCP{ENSEMBLE}[$ens]->{CTD_SVEL})) {
$lastSvel = $LADCP{ENSEMBLE}[$ens]->{CTD_SVEL};
} else {
$LADCP{ENSEMBLE}[$ens]->{CTD_SVEL} = $lastSvel;
}
my($sscorr) = $LADCP{ENSEMBLE}[$ens]->{CTD_SVEL} / $LADCP{ENSEMBLE}[$ens]->{SPEED_OF_SOUND};
for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
next unless defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W]);
$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$U] *= $sscorr;
$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$V] *= $sscorr;
$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W] *= $sscorr;
}
$LADCP{ENSEMBLE}[$ens]->{CTD_LAT} = $CTD{lat}[$r];
$LADCP{ENSEMBLE}[$ens]->{CTD_LON} = $CTD{lon}[$r];
}
&antsAddParams('LADCPproc::min_depth',round($LADCP{ENSEMBLE}[$LADCP_top]->{DEPTH}),
'LADCPproc::max_depth',round($LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH}),
'LADCPproc::start_date',$LADCP{ENSEMBLE}[$LADCP_start]->{DATE},
'LADCPproc::start_time',$LADCP{ENSEMBLE}[$LADCP_start]->{TIME},
'LADCPproc::start_lat',$LADCP{ENSEMBLE}[$LADCP_start]->{CTD_LAT},
'LADCPproc::start_lon',$LADCP{ENSEMBLE}[$LADCP_start]->{CTD_LON},
'LADCPproc::bottom_date',$LADCP{ENSEMBLE}[$LADCP_bottom]->{DATE},
'LADCPproc::bottom_time',$LADCP{ENSEMBLE}[$LADCP_bottom]->{TIME},
'LADCPproc::bottom_lat',$LADCP{ENSEMBLE}[$LADCP_bottom]->{CTD_LAT},
'LADCPproc::bottom_lon',$LADCP{ENSEMBLE}[$LADCP_bottom]->{CTD_LON},
'LADCPproc::end_date',$LADCP{ENSEMBLE}[$LADCP_end]->{DATE},
'LADCPproc::end_time',$LADCP{ENSEMBLE}[$LADCP_end]->{TIME},
'LADCPproc::end_lat',$LADCP{ENSEMBLE}[$LADCP_end]->{CTD_LAT},
'LADCPproc::end_lon',$LADCP{ENSEMBLE}[$LADCP_end]->{CTD_LON});
print(STDERR "\n");
#----------------------------------------------------------------------
# Step 6: Calculate Acoustic Backscatter Profile
#----------------------------------------------------------------------
print(STDERR "Calculating acoustic backscatter profiles...");
mk_backscatter_profs($LADCP_start,$LADCP_end);
print(STDERR "\n");
#----------------------------------------------------------------------
# Step 7: Find Seabed
#----------------------------------------------------------------------
if ($LADCP{ENSEMBLE}[$LADCP_start]->{XDUCER_FACING_DOWN}) {
print(STDERR "Finding seabed...");
print(STDERR "\n\tin acoustic backscatter profiles...") if ($opt_d);
($water_depth,$sig_water_depth) =
find_backscatter_seabed($LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH});
printf(STDERR "\n\t\twater depth = %d(+-%.1f)m",$water_depth,$sig_water_depth)
if ($opt_d);
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)) {
printf(STDERR "\n\t\twater depth = %d(+-%.1f)m",$BT_water_depth,$sig_BT_water_depth)
if ($opt_d);
} else {
printf(STDERR "\n\t\tno seabed found") if ($opt_d);
}
if ($opt_r) { # use RDI BT data
$water_depth = $BT_water_depth;
$sig_water_depth = $sig_BT_water_depth;
}
$min_hab = $water_depth - $LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH};
if (defined($water_depth)) {
printf(STDERR "\n\twater depth = %d(+-%.1f)m",$water_depth,$sig_water_depth);
printf(STDERR "\n\tclosest approach = %dmab",$min_hab);
} else {
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);
}
&antsAddParams('LADCPproc::water_depth',round($water_depth),
'LADCPproc::water_depth.sig',round($sig_water_depth));
#-----------------------------------------------------------------------
# Step 8: Edit Data & also produce depth-time-series output via callback
#-----------------------------------------------------------------------
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');
#--------------------------------------------------------------------------------
# callback routine to output .tds data, called once each for down-/upcasts after
# shear has been gridded. Note that only the nominal (bin center) depth is known
# for each sample at this stage of processing.
#-------------------------------------------------------------------------------
sub outTDseries($)
{
my($downcast) = @_; # also use local [uvw]sh_vals[][]
if ($downcast) {
my($mingi);
for ($mingi=0; $mingi<@ush_vals; $mingi++) {
last if @{$ush_vals[$mingi]};
}
&antsAddParams('LADCPproc::min_ens',$LADCP_start,'LADCPproc::min_elapsed',$LADCP{ENSEMBLE}[$LADCP_start]->{ELAPSED_TIME}+$CTD{first_elapsed}-$opt_l,
'LADCPproc::max_ens',$LADCP_end,'LADCPproc::max_elapsed',$LADCP{ENSEMBLE}[$LADCP_end]->{ELAPSED_TIME}+$CTD{first_elapsed}-$opt_l,
'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]);
}
}
} 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]);
}
}
}
} # sub outDSseries
#--------------------------------------------------------------------------------
# end of callback routine
#--------------------------------------------------------------------------------
$commonParams = $antsCurParams; # PARAMs added up to now are common to all output files
&antsAddParams('LADCPproc::ubin_start',$ubin_start,'LADCPproc::ubin_end',$ubin_end, # record processing params
'LADCPproc::wbin_start',$wbin_start,'LADCPproc::wbin_end',$wbin_end,
'LADCPproc::shbin_start',$shbin_start,'LADCPproc::shbin_end',$shbin_end,
'LADCPproc::w_ref_bin',$w_ref_bin,'LADCPproc::w_dif',$w_dif,
'LADCPproc::wake_hd_dif',$wake_hd_dif,'LADCPproc::wake_ang_min',$wake_ang_min,
'LADCPproc::min_wake_w',$min_wake_w,'LADCPproc::n_wake_bins',$n_wake_bins,
'LADCPproc::e_max',$e_max,'LADCPproc::min_cor',$min_cor,
'LADCPproc::max_shdev',$max_shdev,'LADCPproc::max_shdev_sum',$max_shdev_sum,
'LADCPproc::min_hab',round($min_hab),'LADCPproc::PPI_editing_enabled',$PPI_editing_enabled,
'LADCPproc::clip_margin',$clip_margin,'LADCPproc::first_clip_bin',$first_clip_bin,
'LADCPproc::Svbin_start',$Svbin_start,'LADCPproc::Svbin_end',$Svbin_end,
'LADCPproc::BT_bin_start',$BT_bin_start,'LADCPproc::BT_bin_search_above',$BT_bin_search_above,
'LADCPproc::BT_max_bin_spread',$BT_max_bin_spread,'LADCPproc::BT_max_w_difference',$BT_max_w_difference,
);
if (defined($BT_min_depth)) { # BT-related params
&antsAddParams('LADCPproc::BT_min_depth',$BT_min_depth,'LADCPproc::BT_max_depth',$BT_max_depth);
} else {
&antsAddParams('LADCPproc::BT_max_depth_error',$BT_max_depth_error);
}
$fullParams = $antsCurParams;
$LADCP_start = 1 if ($LADCP_start == 0); # ensure that there is previous ensemble
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
@dc_sh_n = @sh_n; # save downcast results
@dc_ush_mu = @ush_mu; @dc_ush_sig = @ush_sig;
@dc_vsh_mu = @vsh_mu; @dc_vsh_sig = @vsh_sig;
@dc_wsh_mu = @wsh_mu; @dc_wsh_sig = @wsh_sig;
@dc_esh_mu = @esh_mu;
@dc_lash_mu = @lash_mu; @dc_losh_mu = @losh_mu;
print(STDERR "\n\tupcast...") if ($opt_d);
@sh_n=@ush_mu=@ush_sig=@vsh_mu=@vsh_sig=@wsh_mu=@wsh_sig=@esh_mu=@lash_mu=@losh_mu=undef;
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);
@uc_sh_n = @sh_n; # save upcast results
@uc_ush_mu = @ush_mu; @uc_ush_sig = @ush_sig;
@uc_vsh_mu = @vsh_mu; @uc_vsh_sig = @vsh_sig;
@uc_wsh_mu = @wsh_mu; @uc_wsh_sig = @wsh_sig;
@uc_esh_mu = @esh_mu;
@uc_lash_mu = @lash_mu; @uc_losh_mu = @losh_mu;
print(STDERR "\n\tcombined...") if ($opt_d);
my($nsh) = (@dc_ush_mu > @uc_ush_mu) ? scalar(@dc_ush_mu) : scalar(@uc_ush_mu);
for (my($gi)=0; $gi<$nsh; $gi++) {
if (($gi<@dc_ush_mu&&$dc_sh_n[$gi]>0) && ($gi<@uc_ush_mu&&$uc_sh_n[$gi]>0)) {
$sh_n[$gi] = $dc_sh_n[$gi] + $uc_sh_n[$gi];
$ush_mu[$gi] = ($dc_sh_n[$gi]*$dc_ush_mu[$gi] + $uc_sh_n[$gi]*$uc_ush_mu[$gi]) / $sh_n[$gi];
$vsh_mu[$gi] = ($dc_sh_n[$gi]*$dc_vsh_mu[$gi] + $uc_sh_n[$gi]*$uc_vsh_mu[$gi]) / $sh_n[$gi];
$wsh_mu[$gi] = ($dc_sh_n[$gi]*$dc_wsh_mu[$gi] + $uc_sh_n[$gi]*$uc_wsh_mu[$gi]) / $sh_n[$gi];
$esh_mu[$gi] = ($dc_sh_n[$gi]*$dc_esh_mu[$gi] + $uc_sh_n[$gi]*$uc_esh_mu[$gi]) / $sh_n[$gi];
$ush_sig[$gi] = sqrt(($dc_sh_n[$gi]*$dc_ush_sig[$gi]**2 + $uc_sh_n[$gi]*$uc_ush_sig[$gi]**2) / $sh_n[$gi]);
$vsh_sig[$gi] = sqrt(($dc_sh_n[$gi]*$dc_vsh_sig[$gi]**2 + $uc_sh_n[$gi]*$uc_vsh_sig[$gi]**2) / $sh_n[$gi]);
$wsh_sig[$gi] = sqrt(($dc_sh_n[$gi]*$dc_wsh_sig[$gi]**2 + $uc_sh_n[$gi]*$uc_wsh_sig[$gi]**2) / $sh_n[$gi]);
} elsif ($gi < @dc_ush_mu && $dc_sh_n[$gi] > 0) {
$sh_n[$gi] = $dc_sh_n[$gi];
$ush_mu[$gi] = $dc_ush_mu[$gi]; $vsh_mu[$gi] = $dc_vsh_mu[$gi]; $wsh_mu[$gi] = $dc_wsh_mu[$gi]; $esh_mu[$gi] = $dc_esh_mu[$gi];
$ush_sig[$gi] = $dc_ush_sig[$gi]; $vsh_sig[$gi] = $dc_vsh_sig[$gi]; $wsh_sig[$gi] = $dc_wsh_sig[$gi];
} elsif ($gi < @uc_ush_mu && $uc_sh_n[$gi] > 0) {
$sh_n[$gi] = $uc_sh_n[$gi];
$ush_mu[$gi] = $uc_ush_mu[$gi]; $vsh_mu[$gi] = $uc_vsh_mu[$gi]; $wsh_mu[$gi] = $uc_wsh_mu[$gi]; $esh_mu[$gi] = $uc_esh_mu[$gi];
$ush_sig[$gi] = $uc_ush_sig[$gi]; $vsh_sig[$gi] = $uc_vsh_sig[$gi]; $wsh_sig[$gi] = $uc_wsh_sig[$gi];
} else {
$sh_n[$gi] = 0;
$ush_mu[$gi] = $vsh_mu[$gi] = $wsh_mu[$gi] = $esh_mu[$gi] = nan;
$ush_sig[$gi] = $vsh_sig[$gi] = $wsh_sig[$gi] = nan;
}
}
print(STDERR "\n");
#----------------------------------------------------------------------
# Step 9: Get bottom track profile
#----------------------------------------------------------------------
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 Profiles if requested
#----------------------------------------------------------------------
if (defined($opt_p)) {
print(STDERR "Writing shear profiles...");
@antsNewLayout = ('depth','dc_elapsed','dc_lat','dc_lon','dc_nsamp','dc_u_z','dc_u_z.sig','dc_v_z','dc_v_z.sig','dc_w_z','dc_w_z.sig',
'uc_elapsed','uc_lat','uc_lon','uc_nsamp','uc_u_z','uc_u_z.sig','uc_v_z','uc_v_z.sig','uc_w_z','uc_w_z.sig',
'elapsed','nsamp','u_z','u_z.sig','v_z','v_z.sig','w_z','w_z.sig','Sv','Sv.nsamp');
&antsOut('EOF');
close(STDOUT);
open(STDOUT,">$opt_p") || croak("$opt_p: $!\n");
$antsCurParams = $fullParams;
for (my($gi)=0; $gi<$nsh; $gi++) {
&antsOut(depthOfGI($gi), # depth in center of bin
$dc_esh_mu[$gi],$dc_lash_mu[$gi],$dc_losh_mu[$gi], # downcast
numberp($dc_sh_n[$gi])?$dc_sh_n[$gi]:0,
$dc_ush_mu[$gi],$dc_ush_sig[$gi],
$dc_vsh_mu[$gi],$dc_vsh_sig[$gi],
$dc_wsh_mu[$gi],$dc_wsh_sig[$gi],
$uc_esh_mu[$gi],$uc_lash_mu[$gi],$uc_losh_mu[$gi], # upcast
numberp($uc_sh_n[$gi])?$uc_sh_n[$gi]:0,
$uc_ush_mu[$gi],$uc_ush_sig[$gi],
$uc_vsh_mu[$gi],$uc_vsh_sig[$gi],
$uc_wsh_mu[$gi],$uc_wsh_sig[$gi],
$esh_mu[$gi],$sh_n[$gi], # combined
$ush_mu[$gi],$ush_sig[$gi],
$vsh_mu[$gi],$vsh_sig[$gi],
$wsh_mu[$gi],$wsh_sig[$gi],
$nSv_prof[$gi]?$sSv_prof[$gi]/$nSv_prof[$gi]:nan,
$nSv_prof[$gi],
);
}
print(STDERR "\n");
} # if -p
#---------------------------------------
# Acoustic backscatter depth-time-series
#---------------------------------------
if (defined($opt_a)) {
print(STDERR "Writing acoustic backscatter depth-time-series to <$opt_a>...");
@antsNewLayout = ('ens','elapsed','CTD_depth','CTD_lat','CTD_lon','depth','bin','beam','downcast','Sv');
&antsOut('EOF');
close(STDOUT);
open(STDOUT,">$opt_a") || croak("$opt_a: $!\n");
$antsCurParams = $commonParams;
&antsAddParams('LADCPproc::min_elapsed',$LADCP{ENSEMBLE}[$LADCP_start]->{ELAPSED_TIME}+$CTD{first_elapsed}-$opt_l,
'LADCPproc::max_elapsed',$LADCP{ENSEMBLE}[$LADCP_end]->{ELAPSED_TIME}+$CTD{first_elapsed}-$opt_l,
'LADCPproc::min_depth',$LADCP{ENSEMBLE}[$LADCP_top]->{XDUCER_FACING_UP} ?
&depthOfBin($LADCP_top,$LADCP{N_BINS}-1) : $LADCP{ENSEMBLE}[$LADCP_top]->{DEPTH},
'LADCPproc::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++) {
for (my($beam)=0; $beam<4; $beam++) {
&antsOut($LADCP{ENSEMBLE}[$ens]->{NUMBER},
$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}+$CTD{first_elapsed}-$opt_l,
$LADCP{ENSEMBLE}[$ens]->{DEPTH},
$LADCP{ENSEMBLE}[$ens]->{CTD_LAT},$LADCP{ENSEMBLE}[$ens]->{CTD_LON},
&depthOfBinAlongBeam($ens,$bin,$beam),$bin+1,$beam+1,
($ens <= $LADCP_bottom) ? 1 : 0,
$LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$beam],
);
}
}
}
print(STDERR "\n");
}
#----------------------------------------------------------------------
if (defined($opt_t)) {
print(STDERR "Writing time series to <$opt_t>...");
@antsNewLayout = ('ens','elapsed','depth','CTD_lat','CTD_lon','CTD_w','LADCP_w','LADCP_errvel','pitch','roll','hdg');
&antsOut('EOF');
$antsCurParams = $commonParams;
close(STDOUT);
open(STDOUT,">$opt_t") || croak("$opt_t: $!\n");
for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
&antsOut($LADCP{ENSEMBLE}[$ens]->{NUMBER},
$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}+$CTD{first_elapsed}-$opt_l,
$LADCP{ENSEMBLE}[$ens]->{DEPTH},
$LADCP{ENSEMBLE}[$ens]->{CTD_LAT},
$LADCP{ENSEMBLE}[$ens]->{CTD_LON},
$LADCP{ENSEMBLE}[$ens]->{CTD_W},
$LADCP{ENSEMBLE}[$ens]->{W},
$LADCP{ENSEMBLE}[$ens]->{ERR_VEL},
$LADCP{ENSEMBLE}[$ens]->{PITCH},
$LADCP{ENSEMBLE}[$ens]->{ROLL},
$LADCP{ENSEMBLE}[$ens]->{HEADING},
);
}
print(STDERR "\n");
}
#----------------------------------------------------------------------
if (defined($opt_b)) {
print(STDERR "Writing bottom-track data to <$opt_b>...");
@antsNewLayout = ('depth','u','v','w','u.sig','v.sig','w.sig','nsamp');
&antsOut('EOF');
$antsCurParams = $commonParams;
close(STDOUT);
open(STDOUT,">$opt_b") || croak("$opt_b: $!\n");
my($skipped);
for (my($gi)=0; $gi<@BT_nsamp; $gi++) {
$skipped = 1 if ($BT_nsamp[$gi] > 0);
next unless ($skipped);
&antsOut(depthOfGI($gi),$BTu[$gi],$BTv[$gi],$BTw[$gi],$BTu_sig[$gi],$BTv_sig[$gi],$BTw_sig[$gi],$BT_nsamp[$gi]);
}
print(STDERR "\n");
}
#----------------------------------------------------------------------
if (defined($opt_f)) {
print(STDERR "Writing data flags to <$opt_f>...");
@antsNewLayout = ('ens');
for (my($i)=1; $i<=$LADCP{N_BINS}; $i++) {
$antsNewLayout[$i] = "bin$i";
}
&antsOut('EOF');
$antsCurParams = $commonParams;
close(STDOUT);
open(STDOUT,">$opt_f") || croak("$opt_f: $!\n");
&antsPrintHeaders(STDOUT,@antsNewLayout);
for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
printf('%4d ',$LADCP{ENSEMBLE}[$ens]->{NUMBER});
for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
printf("%02x ",$edit_flags[$ens][$bin]);
}
print($opt_R);
}
print(STDERR "\n");
}
&antsExit();