#!/usr/bin/perl
#======================================================================
# L A D C P _ W _ C T D
# doc: Mon Nov 3 17:34:19 2014
# dlm: Mon Apr 29 18:00:44 2019
# (c) 2014 A.M. Thurnherr
# uE-Info: 111 15 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
$antsSummary = 'pre-process SBE 9plus CTD data for LADCP_w';
# HISTORY:
# Nov 3, 2014: - created
# Nov 4, 2014: - improved
# Nov 6, 2014: - BUG: sound speed was not calculated correctly
# - added -a
# - added conductivity & temperature editing
# Nov 7, 2014: - loosened outlier editing
# - added no-valid-data error message
# - modified binning criterion to allow any sampling
# frequency (not just divisors of 24)
# Apr 17, 2015: - added in-situ temperature to output
# Apr 23, 2015: - had to losen unrealistic-soundspeed criterion
# because of DUK2#45
# May 13, 2015: - BUG: did not work with casts beginning at depth
# - added errors when editing edits too much
# May 14, 2015: - BUG: depth was wrong for partial-depth casts
# - BUG: $badval had not been handled correctly
# - BUG: editing errors were too tight
# Jun 16, 2015: - added STDOUT redirection (default output file) on tty stdout
# - added -p)lot option
# Jun 17, 2015: - added cond to std output
# - added conductivity spike filter (different from pressure)
# Jun 18, 2015: - renamed from SBE_w
# Jun 19, 2015: - added press & salin to make output suitable for u/v processing
# - removed a couple of data assertions
# Jul 20, 2015: - allow for non-numeric -i
# Aug 17, 2015: - modified plots
# Sep 21, 2015: - added lat/lon to output if available in input
# Sep 22, 2015: - added error on non-numeric station information
# - changed plot output name
# Sep 23, 2015: - added -o
# - improved pressure editing
# - modified editing algorithm, primarily reducing temperature resolution
# from 0.1deg to 1deg and finding contiguous range up to single gap in histogram
# (requring halving of resolution)
# Oct 12, 2015: - require ANTSlibs V6.2 for release
# 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
# Mar 13, 2016: - added simple ASCII file format
# - BUG: Layout error when input had no lat/lon
# - added lon to simple ASCII format
# Mar 16, 2016: - adapted to gmt5
# Mar 19, 2016: - added support for $VERB environment var
# - added $libSBE_quiet flag
# Mar 30, 2016: - added station number of ASCII format
# Mar 31, 2016: - changed version %PARAM
# May 26, 2016: - renamed w[.raw] to CTD_w[.raw]
# - added winch_w
# Jul 29, 2016: - change w_CTD plot label to reflect sign convention
# Oct 3, 2017: - cosmetics
# - added -q to suppress plots when generating 1Hz files
# Dec 9, 2017: - added $antsSuppressCommonOptions = 1;
# Dec 17, 2017: - added dependencies
# Mar 8, 2018: - made error messages work with -v
# - renamed -l)owpass option to -c)utoff
# - added new -l)atitude
# Mar 9, 2018: - changed -l) to location (include lon as well)
# May 22, 2018: - added code to remove initial at surface interval
# Jul 24, 2018: - improved CTD load error message
# Oct 4, 2018: - minor improvement in log
# - BUG: profiles collected with LOTS of wave heave (AAIW 2005, 004)
# failed the contiguous pressure test at 2dbar resolution
# => increased to 3dbar
# Oct 5, 2018: - BUG: previous bug "fix" did not solve a real problem, because
# problem with AAIW data is lack of pressure resolution
# => returned to 2dbar
# - added plotting errors
# - improved log message
# - BUG: initial in-air scans were not handled correctly (nscans not updated)
# Mar 25, 2019: - changed error message to allow creating 1Hz from 4Hz file (SBE19)
# - BUG: ITS was not set. How is this possible?????
# Apr 21, 2019: - modified code to allow production of 24Hz files (previous code required
# min 2 samples per bin, allowing for max 12Hz sampling rate)
# NOTES:
# w_CTD is positive during the downcast to make the sign of the apparent
# water velocity consistent with w_ocean
($ANTS) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
($WCALC) = ($0 =~ m{^(.*)/[^/]*$});
$WCALC = '.' if ($WCALC eq '');
require "$WCALC/version.pl";
require "$ANTS/ants.pl";
require "$ANTS/fft.pl";
require "$ANTS/libstats.pl";
require "$ANTS/libconv.pl";
require "$ANTS/libGMT.pl";
require "$ANTS/libSBE.pl";
require "$ANTS/libEOS83.pl";
&antsAddParams('LADCP_w_CTD::version',$VERSION);
$antsParseHeader = 0; # usage
$antsSuppressCommonOptions = 1;
&antsUsage('ac:i:l:orp:qs:v:w:',1,
'[-v)erbosity <level[0]>]',
'[use -a)lternate sensor pair]',
'[-r)etain all data (no editing)] [allow infinite -o)utliers]',
'[-s)ampling <rate[6Hz]>]',
'[lowpass w_CTD -c)utoff <limit[2s]>] [-w)inch-speed <granularity[10s]>]',
'[profile -i)d <id>] [station -l)ocation <lat/lon>]',
'[-p)lot_basenames <[%03d_w_CTD.ps],[%03d_sspd.ps]>]',
'[-q)uiet (no plots)]',
'<SBE CNV file>');
&antsFloatOpt(\$opt_c,2); # default low-pass cutoff for w_CTD
&antsCardOpt(\$opt_s,6); # default output sampling rate (Hz)
&antsFloatOpt(\$opt_w,10); # winch velocity granularity
&antsCardOpt(\$opt_v,$ENV{VERB}); # support VERB env variable
$CNVfile = $ARGV[0]; # open CNV file
open(F,&antsFileArg());
&antsAddDeps($CNVfile);
&antsActivateOut(); # activate ANTS file
#----------------------------------------------------------------------
# Read Data
#----------------------------------------------------------------------
sub _croak(@)
{
print(STDERR "\n") if ($opt_v);
croak(@_);
}
print(STDERR "Reading $CNVfile...") if ($opt_v);
chomp($rec = <F>);
_croak("$CNVfile: no data\n")
unless defined($rec);
if ($rec =~ /^\*/) { # SBE CNV file
$libSBE_quiet = 1; # suppress diagnostic messages
($nfields,$nscans,$sampint,$badval,$ftype,$lat,$lon) = # decode SBE header
SBE_parseHeader(F,0,0); # SBE field names, no time check
# _croak("$CNVfile: unexpected sampling interval $sampint s\n")
# unless (abs($sampint-1/24) < 1e-5);
_croak("$CNVfile: insufficient time resolution ($sampint s) for ${opt_s}Hz time series\n")
if (round(1/$sampint/$opt_s) < 1);
if (defined($opt_l)) { # set/override station location with -l
my($slat,$slon) = split('[,/]',$opt_l);
$lat = GMT2deg($slat);
$lon = GMT2deg($slon);
_croak("$0: cannot decode -l $opt_l\n")
unless numberp($lat) && numberp($lon);
}
_croak("$CNVfile: no latitude in header => -l required\n")
unless numberp($lat);
&antsAddParams('lat',$lat);
&antsAddParams('lon',$lon);
$pressF = fnrNoErr('prdM'); # pressure; SBE19plus
$pressF = fnr('prDM') unless defined($pressF); # SBE 9plus
if ($opt_a) { # temp/cond alternate sensor pair
$tempF = fnr('t190C');
&antsAddParams('ITS',90);
$condF = fnrNoErr('c1S/m');
if (defined($condF)) {
$condHistRes = 20; # 0.2 S/m bins
} else {
$condF = fnr('c1mS/cm');
$condHistRes = 2; # 2.0 mS/cm bins
}
} else { # primary sensor pair
$tempF = fnrNoErr('tv290C'); # SBE 19plus
$tempF = fnr('t090C') unless defined($tempF); # SBE 9plus
&antsAddParams('ITS',90);
$condF = fnrNoErr('c0S/m');
if (defined($condF)) {
$condHistRes = 20;
} else {
$condF = fnr('c0mS/cm');
$condHistRes = 2;
}
}
$latF = &fnrNoErr('lat'); # GPS data if available (to make files useful for u/v processing)
$lonF = &fnrNoErr('lon');
&antsInstallBufFull(0); # read entire CNV file
&SBEin(F,$ftype,$nfields,$nscans,$badval);
if (@ants_ != $nscans) {
if ($opt_v > 1) {
printf(STDERR "\n\nWARNING: $CNVfile has wrong number of scans in header (%d instead of %d)\n",$nscans,scalar(@ants_));
} else {
printf(STDERR "WARNING: $CNVfile has wrong number of scans in header (%d instead of %d)\n",$nscans,scalar(@ants_));
}
$nscans = @ants_;
}
} else { # simple CSV ASCII file format:
($lat,$lon,$station) = split(',',$rec);
_croak("$CNVfile: ASCII file format error (1st rec must be lat,lon[,id])\n") # header: lat,lon[,station]
unless numberp($lat) && numberp($lon) &&
$lat>=-90 && $lat<=90 &&
$lon>=-360 && $lon<=360;
&antsAddParams('station',$station) if defined($station);
&antsAddParams('lat',$lat,'lon',$lon);
$sampint = 1/24; $condHistRes = 2; # assumptions: 24Hz
&antsAddParams('ITS',90,'cond.unit','mS/cm'); # ITS-90, mS/cm
$pressF = 1; $tempF = 2; $condF = 3; # Layout: scan press temp cond
for ($nscans=0; <F>; $nscans++) {
chomp;
@{$ants_[$nscans]} = split(','); # CSV format
_croak(sprintf("$CNVfile: unexpected scan #$ants_[$nscans][0] (%d expected)\n",$nscans+1))
unless ($ants_[$nscans][0] == $nscans+1);
$ants_[$nscans][$pressF] = nan unless defined($ants_[$nscans][$pressF]); # missing values
$ants_[$nscans][$tempF] = nan unless defined($ants_[$nscans][$tempF]);
$ants_[$nscans][$condF] = nan unless defined($ants_[$nscans][$condF]);
}
}
printf(STDERR "\n\t%d scans (%d seconds)",$nscans,int($nscans*$sampint))
if ($opt_v > 1);
printf(STDERR "\n") if ($opt_v);
#----------------------------------------------------------------------
# Redirect STDOUT to %.6Hz & create %_w_CTD.ps,%_sspd.ps if STDOUT is a tty
#----------------------------------------------------------------------
$id = defined($opt_i) ? $opt_i : &antsParam('station');
_croak("$CNVfile: no station information in header => -i required\n")
unless defined($id);
_croak("$CNVfile: non-numeric station information <$id> in header => -i required\n")
unless numberp($id);
if (-t STDOUT) {
if (numberp($id)) {
$opt_p = '%03d_w_CTD.ps,%03d_sspd.ps'
unless defined($opt_p);
$outfile = sprintf('%03d.%dHz',$id,$opt_s);
} else {
$opt_p = '%s_w_CTD.ps,%s_sspd.ps'
unless defined($opt_p);
$outfile = sprintf('%s.%dHz',$id,$opt_s);
}
open(STDOUT,">$outfile") || die("$outfile: $!\n");
}
undef($opt_p) if $opt_q; # suppress all plots on -q
#----------------------------------------------------------------------
# Edit Data
# - pressure outliers & spikes
# - conductivity outliers & spikes
# - max_bin+1 is returned to ensure full range
#----------------------------------------------------------------------
sub validRange($)
{
my($guess_bin) = @_;
my($min_bin,$max_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);
}
unless ($opt_r) {
print(STDERR "Editing Data...") if ($opt_v);
#----------------------------------------
# trim initial scans with
# - nan pressure
# - nan conductivity
# - conductivity <= 10 mS/cm
#----------------------------------------
my($trimmed) = 0; # trim leading nan pressures
shift(@ants_),$trimmed++
# ,printf(STDERR "-> p=$ants_[0][$pressF] c=$ants_[0][$condF]\n")
until !@ants_ ||
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));
_croak("\n$CNVfile: no valid scans (wrong conductivity units?)\n")
unless (@ants_);
$nscans -= $trimmed;
printf(STDERR "\n\t%d initial in-air scans trimmed",$trimmed) if ($opt_v > 1);
my($lvp) = $ants_[0][$pressF];
my($lvc) = $ants_[0][$condF];
#-------------------------------------------------------------------------
# edit pressure outliers outside contiguous range
# - 2dbar resolution increased to 3dbar because of 2005 AAIW profile 004
# - histogram shifted by 100dbar to allow for negative values
#-------------------------------------------------------------------------
my($press_rez) = 2;
my($outliers) = my($modeSamp) = 0; my($modeBin,$min,$max); local(@hist);
for (my($s)=0; $s<$nscans; $s++) {
next unless ($ants_[$s][$pressF]>=-100 && $ants_[$s][$pressF]<=6500);
my($b) = ($ants_[$s][$pressF]+100) / $press_rez;
$hist[$b]++;
next unless ($hist[$b] > $modeSamp);
$modeSamp = $hist[$b]; $modeBin = $b;
}
printf(STDERR "\n\tvalid pressure guess: %d dbar (%d samples)",$press_rez*$modeBin-100,$modeSamp)
if ($opt_v > 1);
($min,$max) = validRange($modeBin);
$min = $press_rez*$min-100; $max = $press_rez*$max-100;
for (my($s)=0; $s<$nscans; $s++) {
if ($ants_[$s][$pressF] > $max) { $outliers++; $ants_[$s][$pressF] = nan; }
if ($ants_[$s][$pressF] < $min) { $outliers++; $ants_[$s][$pressF] = nan; }
}
&antsAddParams("pressure_outliers",sprintf("%d",$outliers));
printf(STDERR "\n\tcontinuous pressure range: %d..%d dbar (%d outliers removed)",$min,$max,$outliers) if ($opt_v > 1);
_croak("$CNVfile: pressure editing removed too many 'outliers'\n")
unless ($opt_o || $outliers < 100);
#----------------------------------------------------
# edit conductivity outliers outside contiguous range
#----------------------------------------------------
$outliers = $modeSamp = 0;
undef(@hist);
for (my($s)=0; $s<$nscans; $s++) {
next unless ($ants_[$s][$condF] > 0);
my($b) = $ants_[$s][$condF]*$condHistRes; # 1/10 S/m histogram resolution (1 mS/cm)
$hist[$b]++;
next unless ($hist[$b] > $modeSamp);
$modeSamp = $hist[$b]; $modeBin = $b;
}
($min,$max) = validRange($modeBin);
$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; }
}
&antsAddParams("conductivity_outliers",sprintf("%d",$outliers));
printf(STDERR "\n\tcontinuous conductivity range: %.1f..%.1f S/m (%d outliers removed)",$min,$max,$outliers) if ($opt_v > 1);
_croak("$CNVfile: conductivity editing removed too many 'outliers'\n")
unless ($opt_o || $outliers/$nscans < 0.4);
#----------------------------------------------------
# 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] >= -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/0.5-10,$modeSamp)
# if ($opt_v > 1);
($min,$max) = validRange($modeBin);
$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; }
}
&antsAddParams("temperature_outliers",sprintf("%d",$outliers));
printf(STDERR "\n\tcontinuous temperature range: %.1f..%.1f degC (%d outliers removed)",$min,$max,$outliers)
if ($opt_v > 1);
_croak("$CNVfile: temperature editing removed too many 'outliers'\n")
unless ($opt_o || $outliers/$nscans < 0.4);
#----------------------------------------
# edit pressure spikes based on gradients
#----------------------------------------
for (my($s)=1; $s<$nscans; $s++) { # calculate pressure gradients (across gaps)
if (numberp($ants_[$s][$pressF])) {
$dp[$s-1] = $ants_[$s][$pressF] - $lvp;
$lvp = $ants_[$s][$pressF];
} else {
$dp[$s-1] = nan;
}
}
my($ns1,$ns2) = (0,0);
for (my($s)=0; $s<$nscans-2; $s++) { # consecutive large pressure gradients of opposite sign
if (($dp[$s]*$dp[$s+1] < 0) && # tests return false if either of the dps is not defined
(abs($dp[$s]) > 10) &&
(abs($dp[$s+1]) > 10)) {
$ants_[$s+1][$pressF] = nan;
$dp[$s] = $dp[$s+1] = undef;
$ns1++;
}
}
for (my($s)=0; $s<$nscans-3; $s++) { # 3 consecutive large pressure gradients of opposite sign
if (($dp[$s]>2 && $dp[$s+1]<-4 && $dp[$s+2]>2) ||
($dp[$s]<-2 && $dp[$s+1]>4 && $dp[$s+2]<-2)) {
$ants_[$s+1][$pressF] = $ants_[$s+2][$pressF] = nan;
$dp[$s] = $dp[$s+1] = $dp[$s+2] = undef;
$ns2+=2;
}
}
&antsAddParams("pressure_spikes_removed",sprintf("%d+%d",$ns1,$ns2));
printf(STDERR "\n\t%d+%d pressure spikes removed",$ns1,$ns2) if ($opt_v>1 && $ns1+$ns2>0);
#--------------------------------------------------
# edit conductivity spikes based on large gradients
# - $lvc = $ants_[0] is guaranteed numeric here
#--------------------------------------------------
my($nv) = my($ns) = my($last_dcond) = my($lvcs) = 0;
my($run_start) = my($run_dcond) = 0;
for (my($s)=1; $s<$nscans; $s++) { # calculate conductivity gradients (across gaps)
next unless numberp($ants_[$s][$condF]);
my($dcond) = $ants_[$s][$condF] - $lvc; # integrate gradient across runs
# print(STDERR "ldc/dc/rdc: $last_dcond/$dcond/$run_dcond\n");
if ($last_dcond*$dcond >= 0) { # run is continuing (gradient does not change sign)
$run_dcond += $dcond;
} else { # run has ended
# print(STDERR "new run at $lvcs (run_dcond = $last_dcond)\n");
$run_start = $lvcs;
$run_dcond = $dcond;
}
if (abs($run_dcond) <= 0.02) { # small integrated gradient => okay
$lvc = $ants_[$s][$condF]; # update stored previous values
$lvcs = $s;
$last_dcond = $dcond;
next; # process next scan
}
# print(STDERR "large gradient ($ants_[$run_start][$condF]-$ants_[$s][$condF], $run_dcond) run $run_start-$s\n");
my($i);
my($max_spike_length) = 24;
for ($i=$s; $i<=$run_start+$max_spike_length && $i<$nscans-1; $i++) { # large gradient => check whether values return within 10s
next unless (numberp($ants_[$i][$condF]) &&
numberp($ants_[$i+1][$condF]));
last if ((abs($ants_[$i][$condF]-$ants_[$run_start][$condF]) <= 0.005) && # 2 vals to avoid large gradients
(abs($ants_[$i+1][$condF]-$ants_[$run_start][$condF]) <= 0.005));
}
if ($i>$run_start+$max_spike_length || $i==$nscans-1) { # values don't return => leave data alone
# print(STDERR "values don't return => new run at $s\n");
$run_start = $s; $run_dcond = 0;
$lvc = $ants_[$run_start][$condF]; # start new run
$lvcs = $run_start;
$last_dcond = 0;
next; # process next scan
}
$ns++;
for (my($j)=$run_start+1; $j<$i; $j++) { # values return => remove bad spike
$ants_[$j][$condF] = nan;
$nv++;
}
# print(STDERR "values return at $i ($ants_[$i][$condF]) => deleted $run_start+1-$i; new run\n");
$run_start = $i; $run_dcond = 0;
$lvc = $ants_[$run_start][$condF];
$lvcs = $run_start;
$last_dcond = 0;
}
&antsAddParams("conductivity_spikes_removed",$ns);
printf(STDERR "\n\t%d conductivity values removed from %d spikes",$nv,$ns) if ($opt_v > 1);
printf(STDERR "\n") if ($opt_v);
} # if $opt_r
#----------------------------------------------------------------------
# Correcting for pressure bias
#----------------------------------------------------------------------
print(STDERR "Correcting for pressure bias...") if ($opt_v);
my($minP) = 9e99;
for (my($s)=0; $s<$nscans; $s++) {
$minP = $ants_[$s][$pressF]
if numberp($ants_[$s][$pressF]) && ($ants_[$s][$pressF] < $minP);
}
_croak("$CNVfile: no valid CTD pressure data below 25dbar\n")
unless ($minP < 9e99);
if ($minP < 25) {
&antsAddParams('pressure_bias',$minP);
printf(STDERR "\n\tsubtracting %.1f dbar",$minP) if ($opt_v > 1);
for (my($s)=0; $s<$nscans; $s++) {
$ants_[$s][$pressF] -= $minP
if numberp($ants_[$s][$pressF]);
}
} else {
printf(STDERR "\n\tpartial-depth cast below %.1f dbar (no correction applied)",$minP) if ($opt_v > 1);
}
printf(STDERR "\n") if ($opt_v);
#----------------------------------------------------------------------
# Removing Initial At-Surface Data
#----------------------------------------------------------------------
print(STDERR "Removing intial at-surface data...") if ($opt_v);
my($trimmed) = 0;
while ($ants_[$trimmed][$pressF] < 0.5) { $trimmed++; }
for (my($r)=$trimmed; $r<$nscans; $r++) {
$ants_[$r][$elapsedF] -= $ants_[$trimmed][$elapsedF];
}
splice(@ants_,0,$trimmed);
$nscans -= $trimmed;
&antsAddParams('surface_data_trimmed',int($trimmed*$sampint));
&antsAddParams('cast_duration',int($nscans*$sampint));
printf(STDERR "\n\t%d seconds of data trimmed",int($trimmed*$sampint)) if ($opt_v > 1);
printf(STDERR "\n") if ($opt_v);
#----------------------------------------------------------------------
# Binning data
#----------------------------------------------------------------------
my($sps) = round(1 / $sampint / $opt_s);
print(STDERR "Creating ${opt_s}Hz time series ($sps samples per bin)...") if ($opt_v);
&antsAddParams('sampling_interval',1/$opt_s);
&antsAddParams('sampling_frequency',$opt_s);
my(@press,@temp,@cond);
my($sp,$np,$st,$nt,$sc,$nc);
$sp = $st = $sc = $np = $nt = $nc = 0;
$slat = $slon = $nGPS = 0;
for (my($rec)=1,my($s)=0; $s<$nscans; $s++) {
if ($s*$sampint > $rec/$opt_s) {
$rec++;
push(@press,$np>0?$sp/$np:nan);
push(@temp, $nt>0?$st/$nt:nan);
push(@cond, $nc>0?$sc/$nc:nan);
$sp = $st = $sc = $np = $nt = $nc = 0;
if (defined($latF)) {
push(@lat,$nGPS>0?$slat/$nGPS:nan);
push(@lon,$nGPS>0?$slon/$nGPS:nan);
$slat = $slon = $nGPS = 0;
}
}
$sp+=$ants_[$s][$pressF],$np++ if numberp($ants_[$s][$pressF]);
$st+=$ants_[$s][$tempF],$nt++ if numberp($ants_[$s][$tempF]);
$sc+=$ants_[$s][$condF],$nc++ if numberp($ants_[$s][$condF]);
if (defined($latF) && numberp($ants_[$s][$latF])) {
$slat += $ants_[$s][$latF];
$slon += $ants_[$s][$lonF];
$nGPS++;
}
}
printf(STDERR "\n") if ($opt_v);
#----------------------------------------------------------------------
# Calculating derived quantities
#----------------------------------------------------------------------
print(STDERR "Calculating vertical package velocity & sound speed...") if ($opt_v);
my($maxP) = -9e99; my($atBtm);
my($min_sspd) = 9e99; my($max_sspd) = -9e99;
for (my($r)=0; $r<@press; $r++) {
$maxP=$press[$r],$atBtm=$r if ($press[$r] > $maxP);
$elapsed[$r] = $r/$opt_s;
$depth[$r] = &depth($press[$r],$lat);
# print(STDERR "depth[$r] = depth($press[$r],$lat) = $depth[$r]\n");
$salin[$r] = &salin($cond[$r],$temp[$r],$press[$r]);
$sspd[$r] = &sVel($salin[$r],$temp[$r],$press[$r]);
$min_sspd = $sspd[$r] if ($sspd[$r] < $min_sspd);
$max_sspd = $sspd[$r] if ($sspd[$r] > $max_sspd);
}
$w[0] = nan;
for (my($r)=1; $r<@depth-1; $r++) {
$w[$r] = numbersp($depth[$r-1],$depth[$r+1])
? ($depth[$r+1] - $depth[$r-1]) * $opt_s
: nan;
}
push(@w,nan);
printf(STDERR "\n") if ($opt_v);
#----------------------------------------------------------------------
# Low-pass filter velocity data
# - interpolate missing vertical velocities first
#----------------------------------------------------------------------
if ($opt_c > 0) {
print(STDERR "Low-pass filtering vertical package velocity...") if ($opt_v);
&antsAddParams('w_lowpass_cutoff',$opt_c);
my($trimmed) = 0;
shift(@w),shift(@depth),shift(@elapsed),shift(@sspd),$trimmed++
until numberp($w[0]);
my($interpolated) = 0;
for ($r=1; $r<@w; $r++) {
next if numberp($w[$r]);
my($lv) = $r-1;
for ($nv=$r+1; $nv<@depth && !numberp($w[$nv]); $nv++) {}
if ($nv < @depth) {
while ($r < $nv) {
$w[$r] = $w[$lv] + ($r-$lv)/($nv-$lv) * ($w[$nv]-$w[$lv]);
$interpolated++;
$r++;
}
} else {
$trimmed += @w-$r;
splice(@w,$r); splice(@depth,$r);
splice(@elapsed,$r); splice(@sspd,$r);
}
}
&antsAddParams('w_interpolated',$interpolated);
printf(STDERR "\n\t%d/%d vertical velocities trimmed/interpolated",$trimmed,$interpolated) if ($opt_v > 1);
#--------------------
# Zero Pad Data
#--------------------
for ($pot=1; $pot<@w; $pot<<=1) {} # determine power of two
for ($r=0; $r<@w; $r++) { # copy data
$fftbuf[2*$r] = $w[$r];
$fftbuf[2*$r+1] = 0;
}
printf(STDERR "\n\tzero-padded %d scans",$pot-$r) if ($opt_v > 1);
while ($r < $pot) { # pad with zeroes
$fftbuf[2*$r] = $fftbuf[2*$r+1] = 0;
$r++;
}
#--------------------
# Low-Pass Filter
#--------------------
@fco = &FOUR1(-1,@fftbuf); # forward FFT
$n = @fco/2;
for (my($ip)=2; $ip<=$n; $ip+=2) { # +ve freq fco
my($in) = 2*$n-$ip; # -ve freq fco
my($f) = $ip/2/$n*$opt_s; # frequency
$fco[$ip] = $fco[$ip+1] = $fco[$in] = $fco[$in+1] = 0
if ($f > 1/$opt_c); # low-pass filter
}
@w_lp = &FOUR1(1,@fco); # inverse FFT
printf(STDERR "\n") if ($opt_v);
} else {
@w_lp = @w;
}
#----------------------------------------
# Estimate winch speed
#----------------------------------------
print(STDERR "Estimating winch velocity...") if ($opt_v);
&antsAddParams('winch_velocity_granularity',$opt_w);
my($from_r) = 0; my($to_r); # step 1: bin average in time
for (my($from_r)=my($to_r)=0; $from_r<@elapsed; $from_r=$to_r+1) {
my($sumw) = $w_lp[2*$from_r]/@w_lp; my($n) = 1;
for ($to_r=$from_r+1; $to_r<@elapsed && $elapsed[$to_r]-$elapsed[$from_r]<$opt_w; $to_r++) {
$sumw += $w_lp[2*$to_r]/@w_lp; $n++;
}
$winch[$from_r] = $sumw/$n;
}
my($pwinch) = $winch[0];
for (my($to_r),my($from_r)=0; $from_r<@elapsed; ) { # step 2: fill after median filtering
for ($to_r=$from_r+1; $to_r<@elapsed && !defined($winch[$to_r]); $to_r++) {}
my($nwinch) = $to_r<@elapsed ? $winch[$to_r] : $winch[$from_r];
my($winch) = median($pwinch,$winch[$from_r],$nwinch);
$pwinch = $winch[$from_r]; $winch[$from_r] = $winch;
while (++$from_r < $to_r) { $winch[$from_r] = $winch[$from_r-1]; }
}
printf(STDERR "\n") if ($opt_v);
#----------------------------------------
# Plot Data
#----------------------------------------
if (defined($opt_p)) {
print(STDERR "Plotting data...\n") if ($opt_v);
my(@pfmt) = split(',',$opt_p);
_croak("$0: cannot decode -p $opt_p\n")
unless (@pfmt == 2);
my($xmin) = $elapsed[0]/60;
my($xmax) = $elapsed[$#elapsed]/60;
my($ymin) = -3; my($ymax) = 3;
my($plotsize) = '13c';
_croak(sprintf("%s: invalid region of interest (-R$xmin/$xmax/$ymin/$ymax)\n",sprintf($pfmt[0],$id)))
unless ($xmax > $xmin && $ymax > $ymin);
GMT_begin(sprintf($pfmt[0],$id),"-JX${plotsize}","-R$xmin/$xmax/$ymin/$ymax",'-X6 -Y4 -P');
GMT_psxy('-W1,coral');
for ($r=0; $r<@w; $r++) {
printf(GMT "%f %f\n",$elapsed[$r]/60,$w_lp[2*$r]/@w_lp);
GMT_psxy('-W1,SeaGreen') if ($r == $atBtm);
}
GMT_psxy('-W1,magenta');
for ($r=0; $r<@w; $r++) {
printf(GMT "%f %f\n",$elapsed[$r]/60,$winch[$r]);
}
GMT_psbasemap('-Bg60a30f5:"Elapsed Time [min]":/g1a1f0.1:"Downward Package Velocity [ms@+-1@+]":WeSn');
GMT_unitcoords();
GMT_pstext('-F+f14,Helvetica,coral+jBR'); print(GMT "0.98 0.96 dc\n");
GMT_pstext('-F+f14,Helvetica,SeaGreen+jBR'); print(GMT "0.98 0.93 uc\n");
GMT_pstext('-F+f14,Helvetica,magenta+jBR'); print(GMT "0.98 0.89 winch\n");
GMT_pstext('-F+f14,Helvetica,blue+jTL -N');
if (defined($outfile)) { printf(GMT "0.01 1.06 $outfile\n",$id); }
else { printf(GMT "0.01 1.06 %03d\n",$id); }
GMT_end();
my($xmin) = round($min_sspd-3,5);
my($xmax) = round($max_sspd+3,5);
my($ymin) = 0; my($ymax) = round($depth[$atBtm]+70,100);
my($plotsize) = '13c';
_croak(sprintf("%s: invalid region of interest (-R$xmin/$xmax/$ymin/$ymax)\n",sprintf($pfmt[1],$id)))
unless ($xmax > $xmin && $ymax > $ymin);
GMT_begin(sprintf($pfmt[1],$id),"-JX${plotsize}/-${plotsize}","-R$xmin/$xmax/$ymin/$ymax",'-X6 -Y4 -P');
GMT_psbasemap('-Bg10a10f1:"Speed of Sound [m/s]":/g1000a500f100:"Depth [m]":WeSn');
GMT_psxy('-W2,coral');
for ($r=0; $r<@w; $r++) {
printf(GMT "%f %f\n",$sspd[$r],$depth[$r]);
GMT_psxy('-W1.5,SeaGreen') if ($r == $atBtm);
}
GMT_unitcoords();
GMT_pstext('-F+f14,Helvetica,coral+jTR'); print(GMT "0.98 0.02 dc\n");
GMT_pstext('-F+f14,Helvetica,SeaGreen+jTR'); print(GMT "0.98 0.06 uc\n");
GMT_pstext('-F+f14,Helvetica,blue+jTL -N');
if (defined($outfile)) { printf(GMT "0.01 -0.06 $outfile\n",$id); }
else { printf(GMT "0.01 -0.06 %03d\n",$id); }
GMT_end();
}
#----------------------------------------------------------------------
print(STDERR "Writing output...\n") if ($opt_v);
@antsNewLayout = ('elapsed','press','temp','cond','depth','salin','sspd','w_CTD.raw','w_CTD','w_winch','lat','lon');
for ($r=0; $r<@w; $r++) {
&antsOut($elapsed[$r],$press[$r],$temp[$r],$cond[$r],$depth[$r],$salin[$r],
$sspd[$r],$w[$r],$w_lp[2*$r]/@w_lp,$winch[$r],
$lat[$r],$lon[$r]);
}
exit(0); # don't flush @ants_