author | A.M. Thurnherr <athurnherr@yahoo.com> |
Tue, 02 Feb 2016 14:55:24 +0000 | |
changeset 34 | e550db661c17 |
parent 33 | 866c881b3a4a |
child 35 | 54b8bb450e5f |
permissions | -rwxr-xr-x |
18 | 1 |
#!/usr/bin/perl |
2 |
#====================================================================== |
|
29 | 3 |
# L A D C P _ W _ C T D |
18 | 4 |
# doc: Mon Nov 3 17:34:19 2014 |
34 | 5 |
# dlm: Mon Jan 25 09:32:32 2016 |
18 | 6 |
# (c) 2014 A.M. Thurnherr |
34 | 7 |
# uE-Info: 51 0 NIL 0 0 72 2 2 4 NIL ofnI |
18 | 8 |
#====================================================================== |
9 |
||
10 |
$antsSummary = 'pre-process SBE 9plus CTD data for LADCP_w'; |
|
11 |
||
12 |
# HISTORY: |
|
13 |
# Nov 3, 2014: - created |
|
14 |
# Nov 4, 2014: - improved |
|
15 |
# Nov 6, 2014: - BUG: sound speed was not calculated correctly |
|
16 |
# - added -a |
|
17 |
# - added conductivity & temperature editing |
|
18 |
# Nov 7, 2014: - loosened outlier editing |
|
19 |
# - added no-valid-data error message |
|
20 |
# - modified binning criterion to allow any sampling |
|
21 |
# frequency (not just divisors of 24) |
|
28
b07b23485336
after folding in EGU flight mods (empirical Sv correction)
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
26
diff
changeset
|
22 |
# Apr 17, 2015: - added in-situ temperature to output |
29 | 23 |
# Apr 23, 2015: - had to losen unrealistic-soundspeed criterion |
24 |
# because of DUK2#45 |
|
25 |
# May 13, 2015: - BUG: did not work with casts beginning at depth |
|
26 |
# - added errors when editing edits too much |
|
27 |
# May 14, 2015: - BUG: depth was wrong for partial-depth casts |
|
28 |
# - BUG: $badval had not been handled correctly |
|
29 |
# - BUG: editing errors were too tight |
|
30 |
# Jun 16, 2015: - added STDOUT redirection (default output file) on tty stdout |
|
31 |
# - added -p)lot option |
|
32 |
# Jun 17, 2015: - added cond to std output |
|
33 |
# - added conductivity spike filter (different from pressure) |
|
34 |
# Jun 18, 2015: - renamed from SBE_w |
|
35 |
# Jun 19, 2015: - added press & salin to make output suitable for u/v processing |
|
36 |
# - removed a couple of data assertions |
|
37 |
# Jul 20, 2015: - allow for non-numeric -i |
|
31 | 38 |
# Aug 17, 2015: - modified plots |
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
39 |
# Sep 21, 2015: - added lat/lon to output if available in input |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
40 |
# Sep 22, 2015: - added error on non-numeric station information |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
41 |
# - changed plot output name |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
42 |
# Sep 23, 2015: - added -o |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
43 |
# - improved pressure editing |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
44 |
# - modified editing algorithm, primarily reducing temperature resolution |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
45 |
# from 0.1deg to 1deg and finding contiguous range up to single gap in histogram |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
46 |
# (requring halving of resolution) |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
47 |
# Oct 12, 2015: - require ANTSlibs V6.2 for release |
33 | 48 |
# Oct 13, 2015: - added code to deal with CNV files with wrong number of scans in header |
49 |
# Oct 13, 2015: - adapted to [version.pl] |
|
34 | 50 |
# Jan 25, 2016: - added software version %PARAM |
18 | 51 |
|
33 | 52 |
($ANTS) = (`which ANTSlib` =~ m{^(.*)/[^/]*$}); |
53 |
($WCALC) = ($0 =~ m{^(.*)/[^/]*$}); |
|
54 |
$WCALC = '.' if ($WCALC eq ''); |
|
18 | 55 |
|
33 | 56 |
require "$WCALC/version.pl"; |
18 | 57 |
require "$ANTS/ants.pl"; |
58 |
require "$ANTS/fft.pl"; |
|
29 | 59 |
require "$ANTS/libGMT.pl"; |
18 | 60 |
require "$ANTS/libSBE.pl"; |
61 |
require "$ANTS/libEOS83.pl"; |
|
34 | 62 |
&antsAddParams('LADCP_w_CTD',"Version $VERSION"); |
18 | 63 |
|
64 |
$antsParseHeader = 0; # usage |
|
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
65 |
&antsUsage('ai:l:orp:s:v:',1, |
18 | 66 |
'[-v)erbosity <level[0]>]', |
67 |
'[use -a)lternate sensor pair]', |
|
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
68 |
'[-r)etain all data (no editing)] [allow infinite -o)utliers]', |
18 | 69 |
'[-s)ampling <rate[6Hz]>]', |
70 |
'[-l)owpass w <cutoff[2s]>]', |
|
29 | 71 |
'[profile -i)d <id>]', |
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
72 |
'[-p)lot_basenames <[%03d_w_CTD.ps],[%03d_sspd.ps]>]', |
18 | 73 |
'<SBE CNV file>'); |
74 |
||
75 |
&antsFloatOpt(\$opt_l,2); # defaults |
|
76 |
&antsCardOpt(\$opt_s,6); |
|
77 |
&antsCardOpt(\$opt_v,0); |
|
78 |
||
79 |
$CNVfile = $ARGV[0]; # open CNV file |
|
80 |
open(F,&antsFileArg()); |
|
28
b07b23485336
after folding in EGU flight mods (empirical Sv correction)
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
26
diff
changeset
|
81 |
&antsActivateOut(); # activate ANTS file |
18 | 82 |
|
83 |
#---------------------------------------------------------------------- |
|
84 |
# Read Data |
|
85 |
#---------------------------------------------------------------------- |
|
86 |
||
87 |
print(STDERR "Reading $CNVfile...") if ($opt_v); |
|
88 |
||
89 |
($nfields,$nscans,$sampint,$badval,$ftype,$lat,$lon) = # decode SBE header |
|
90 |
SBE_parseHeader(F,0,0); # SBE field names, no time check |
|
91 |
||
92 |
croak("$CNVfile: unexpected sampling interval $sampint\n") |
|
93 |
unless (abs($sampint-1/24) < 1e-5); |
|
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
94 |
croak("$CNVfile: no latitude in header\n") |
18 | 95 |
unless numberp($lat); |
96 |
||
97 |
$pressF = fnr('prDM'); |
|
98 |
||
99 |
if ($opt_a) { # temp/cond alternate sensor pair |
|
100 |
$tempF = fnr('t190C'); |
|
101 |
$condF = fnrNoErr('c1S/m'); |
|
102 |
if (defined($condF)) { |
|
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
103 |
$condHistRes = 20; # 0.2 S/m bins |
18 | 104 |
} else { |
105 |
$condF = fnr('c1mS/cm'); |
|
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
106 |
$condHistRes = 2; # 2.0 mS/cm bins |
18 | 107 |
} |
108 |
} else { # primary sensor pair |
|
109 |
$tempF = fnr('t090C'); |
|
110 |
$condF = fnrNoErr('c0S/m'); |
|
111 |
if (defined($condF)) { |
|
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
112 |
$condHistRes = 20; |
18 | 113 |
} else { |
114 |
$condF = fnr('c0mS/cm'); |
|
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
115 |
$condHistRes = 2; |
18 | 116 |
} |
117 |
} |
|
118 |
||
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
119 |
$latF = &fnrNoErr('lat'); # GPS data if available (to make files useful for u/v processing) |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
120 |
$lonF = &fnrNoErr('lon'); |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
121 |
|
18 | 122 |
&antsInstallBufFull(0); # read entire CNV file |
29 | 123 |
&SBEin(F,$ftype,$nfields,$nscans,$badval); |
18 | 124 |
|
33 | 125 |
if (@ants_ != $nscans) { |
126 |
if ($opt_v > 1) { |
|
127 |
printf(STDERR "\n\nWARNING: $CNVfile has wrong number of scans in header\n"); |
|
128 |
} else { |
|
129 |
printf(STDERR "WARNING: $CNVfile has wrong number of scans in header\n"); |
|
130 |
} |
|
131 |
$nscans = @ants_; |
|
132 |
} |
|
133 |
||
18 | 134 |
printf(STDERR "\n\t%d scans",$nscans) if ($opt_v > 1); |
135 |
printf(STDERR "\n") if ($opt_v); |
|
136 |
||
137 |
#---------------------------------------------------------------------- |
|
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
138 |
# Redirect STDOUT to %.6Hz & create %_w_CTD.ps,%_sspd.ps if STDOUT is a tty |
29 | 139 |
#---------------------------------------------------------------------- |
140 |
||
141 |
$id = defined($opt_i) ? $opt_i : &antsParam('station'); |
|
142 |
croak("$CNVfile: no station information in header => -i required\n") |
|
143 |
unless defined($id); |
|
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
144 |
croak("$CNVfile: non-numeric station information <$id> in header => -i required\n") |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
145 |
unless numberp($id); |
29 | 146 |
|
147 |
if (-t STDOUT) { |
|
148 |
if (numberp($id)) { |
|
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
149 |
$opt_p = '%03d_w_CTD.ps,%03d_sspd.ps' |
29 | 150 |
unless defined($opt_p); |
151 |
$outfile = sprintf('%03d.%dHz',$id,$opt_s); |
|
152 |
} else { |
|
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
153 |
$opt_p = '%s_w_CTD.ps,%s_sspd.ps' |
29 | 154 |
unless defined($opt_p); |
155 |
$outfile = sprintf('%s.%dHz',$id,$opt_s); |
|
156 |
} |
|
157 |
open(STDOUT,">$outfile") || die("$outfile: $!\n"); |
|
158 |
} |
|
159 |
||
160 |
#---------------------------------------------------------------------- |
|
18 | 161 |
# Edit Data |
162 |
# - pressure outliers & spikes |
|
163 |
# - conductivity outliers & spikes |
|
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
164 |
# - max_bin+1 is returned to ensure full range |
18 | 165 |
#---------------------------------------------------------------------- |
166 |
||
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
167 |
sub validRange($) |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
168 |
{ |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
169 |
my($guess_bin) = @_; |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
170 |
my($min_bin,$max_bin); |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
171 |
|
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
172 |
croak("assertion failed") unless ($hist[$guess_bin]); |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
173 |
for ($max_bin=$guess_bin; $hist[$max_bin]; $max_bin++) { } |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
174 |
for ($min_bin=$guess_bin; $hist[$min_bin]; $min_bin--) { } |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
175 |
return ($min_bin,$max_bin+1); |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
176 |
} |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
177 |
|
18 | 178 |
unless ($opt_r) { |
179 |
print(STDERR "Editing Data...") if ($opt_v); |
|
180 |
||
181 |
#---------------------------------------- |
|
182 |
# trim initial records with nan pressures |
|
183 |
#---------------------------------------- |
|
184 |
my($trimmed) = 0; # trim leading nan pressures |
|
185 |
shift(@ants_),$trimmed++ |
|
29 | 186 |
until numberp($ants_[0][$pressF]) && numberp($ants_[0][$condF]); |
18 | 187 |
printf(STDERR "\n\t%d initial records with nan pressure and/or conductivity trimmed",$trimmed) if ($opt_v > 1); |
188 |
my($lvp) = $ants_[0][$pressF]; |
|
29 | 189 |
my($lvc) = $ants_[0][$condF]; |
18 | 190 |
|
191 |
#------------------------------------------------ |
|
192 |
# edit pressure outliers outside contiguous range |
|
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
193 |
# - 2dbar resolution |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
194 |
# - histogram shifted by 100dbar to allow for negative values |
18 | 195 |
#------------------------------------------------ |
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
196 |
my($outliers) = my($modeSamp) = 0; my($modeBin,$min,$max); local(@hist); |
18 | 197 |
for (my($s)=0; $s<$nscans; $s++) { |
29 | 198 |
next unless ($ants_[$s][$pressF]>=-100 && $ants_[$s][$pressF]<=6500); |
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
199 |
my($b) = ($ants_[$s][$pressF]+100) / 2; |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
200 |
$hist[$b]++; |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
201 |
next unless ($hist[$b] > $modeSamp); |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
202 |
$modeSamp = $hist[$b]; $modeBin = $b; |
18 | 203 |
} |
33 | 204 |
printf(STDERR "\n\tvalid pressure guess: %d dbar (%d samples)",2*$modeBin-100,$modeSamp) |
205 |
if ($opt_v > 1); |
|
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
206 |
($min,$max) = validRange($modeBin); |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
207 |
$min = 2*$min-100; $max = 2*$max-100; |
18 | 208 |
for (my($s)=0; $s<$nscans; $s++) { |
209 |
if ($ants_[$s][$pressF] > $max) { $outliers++; $ants_[$s][$pressF] = nan; } |
|
210 |
if ($ants_[$s][$pressF] < $min) { $outliers++; $ants_[$s][$pressF] = nan; } |
|
211 |
} |
|
212 |
&antsAddParams("pressure_outliers",sprintf("%d",$outliers)); |
|
213 |
printf(STDERR "\n\tcontinuous pressure range: %d..%d dbar (%d outliers removed)",$min,$max,$outliers) if ($opt_v > 1); |
|
29 | 214 |
croak("$CNVfile: pressure editing removed too many 'outliers'\n") |
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
215 |
unless ($opt_o || $outliers < 100); |
18 | 216 |
|
217 |
#---------------------------------------------------- |
|
218 |
# edit conductivity outliers outside contiguous range |
|
219 |
#---------------------------------------------------- |
|
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
220 |
$outliers = $modeSamp = 0; |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
221 |
undef(@hist); |
18 | 222 |
for (my($s)=0; $s<$nscans; $s++) { |
223 |
next unless ($ants_[$s][$condF] > 0); |
|
224 |
my($b) = $ants_[$s][$condF]*$condHistRes; # 1/10 S/m histogram resolution (1 mS/cm) |
|
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
225 |
$hist[$b]++; |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
226 |
next unless ($hist[$b] > $modeSamp); |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
227 |
$modeSamp = $hist[$b]; $modeBin = $b; |
18 | 228 |
} |
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
229 |
($min,$max) = validRange($modeBin); |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
230 |
$min /= $condHistRes; $max /= $condHistRes; |
18 | 231 |
for (my($s)=0; $s<$nscans; $s++) { |
232 |
if ($ants_[$s][$condF] > $max) { $outliers++; $ants_[$s][$condF] = nan; } |
|
233 |
if ($ants_[$s][$condF] < $min) { $outliers++; $ants_[$s][$condF] = nan; } |
|
234 |
} |
|
235 |
&antsAddParams("conductivity_outliers",sprintf("%d",$outliers)); |
|
236 |
printf(STDERR "\n\tcontinuous conductivity range: %.1f..%.1f S/m (%d outliers removed)",$min,$max,$outliers) if ($opt_v > 1); |
|
29 | 237 |
croak("$CNVfile: conductivity editing removed too many 'outliers'\n") |
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
238 |
unless ($opt_o || $outliers/$nscans < 0.4); |
18 | 239 |
|
240 |
#---------------------------------------------------- |
|
241 |
# edit temperature outliers outside contiguous range |
|
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
242 |
# - Stan's NBP0901 profiles require resolution of 1deg |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
243 |
# - otherwise 0.2deg seems to be fine |
18 | 244 |
#---------------------------------------------------- |
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
245 |
$outliers = $modeSamp = 0; |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
246 |
undef(@hist); |
18 | 247 |
for (my($s)=0; $s<$nscans; $s++) { |
248 |
next unless ($ants_[$s][$tempF] > 0); |
|
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
249 |
my($b) = $ants_[$s][$tempF]*1; # 10th of a degree histogram resolution |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
250 |
$hist[$b]++; |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
251 |
next unless ($hist[$b] > $modeSamp); |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
252 |
$modeSamp = $hist[$b]; $modeBin = $b; |
18 | 253 |
} |
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
254 |
printf(STDERR "\n\ttemperature mode: %.1f degC (%d samples)",$modeBin/10,$modeSamp) |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
255 |
if ($opt_v > 1); |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
256 |
($min,$max) = validRange($modeBin); |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
257 |
$min /= 1; $max /= 1; |
18 | 258 |
for (my($s)=0; $s<$nscans; $s++) { |
259 |
if ($ants_[$s][$tempF] > $max) { $outliers++; $ants_[$s][$tempF] = nan; } |
|
260 |
if ($ants_[$s][$tempF] < $min) { $outliers++; $ants_[$s][$tempF] = nan; } |
|
261 |
} |
|
262 |
&antsAddParams("temperature_outliers",sprintf("%d",$outliers)); |
|
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
263 |
printf(STDERR "\n\tcontinuous temperature range: %.1f..%.1f degC (%d outliers removed)",$min,$max,$outliers) |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
264 |
if ($opt_v > 1); |
29 | 265 |
croak("$CNVfile: temperature editing removed too many 'outliers'\n") |
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
266 |
unless ($opt_o || $outliers/$nscans < 0.4); |
18 | 267 |
|
268 |
#---------------------------------------- |
|
269 |
# edit pressure spikes based on gradients |
|
270 |
#---------------------------------------- |
|
271 |
||
29 | 272 |
for (my($s)=1; $s<$nscans; $s++) { # calculate pressure gradients (across gaps) |
18 | 273 |
if (numberp($ants_[$s][$pressF])) { |
274 |
$dp[$s-1] = $ants_[$s][$pressF] - $lvp; |
|
275 |
$lvp = $ants_[$s][$pressF]; |
|
276 |
} else { |
|
277 |
$dp[$s-1] = nan; |
|
278 |
} |
|
279 |
} |
|
280 |
||
281 |
my($ns1,$ns2) = (0,0); |
|
282 |
for (my($s)=0; $s<$nscans-2; $s++) { # consecutive large pressure gradients of opposite sign |
|
283 |
if (($dp[$s]*$dp[$s+1] < 0) && # tests return false if either of the dps is not defined |
|
284 |
(abs($dp[$s]) > 10) && |
|
285 |
(abs($dp[$s+1]) > 10)) { |
|
286 |
$ants_[$s+1][$pressF] = nan; |
|
287 |
$dp[$s] = $dp[$s+1] = undef; |
|
288 |
$ns1++; |
|
289 |
} |
|
290 |
} |
|
291 |
for (my($s)=0; $s<$nscans-3; $s++) { # 3 consecutive large pressure gradients of opposite sign |
|
292 |
if (($dp[$s]>2 && $dp[$s+1]<-4 && $dp[$s+2]>2) || |
|
293 |
($dp[$s]<-2 && $dp[$s+1]>4 && $dp[$s+2]<-2)) { |
|
294 |
$ants_[$s+1][$pressF] = $ants_[$s+2][$pressF] = nan; |
|
295 |
$dp[$s] = $dp[$s+1] = $dp[$s+2] = undef; |
|
296 |
$ns2+=2; |
|
297 |
} |
|
298 |
} |
|
299 |
&antsAddParams("pressure_spikes_removed",sprintf("%d+%d",$ns1,$ns2)); |
|
300 |
printf(STDERR "\n\t%d+%d pressure spikes removed",$ns1,$ns2) if ($opt_v > 1); |
|
301 |
||
29 | 302 |
#-------------------------------------------------- |
303 |
# edit conductivity spikes based on large gradients |
|
304 |
# - $lvc = $ants_[0] is guaranteed numeric here |
|
305 |
#-------------------------------------------------- |
|
306 |
||
307 |
my($nv) = my($ns) = my($last_dcond) = my($lvcs) = 0; |
|
308 |
my($run_start) = my($run_dcond) = 0; |
|
309 |
for (my($s)=1; $s<$nscans; $s++) { # calculate conductivity gradients (across gaps) |
|
310 |
next unless numberp($ants_[$s][$condF]); |
|
311 |
||
312 |
my($dcond) = $ants_[$s][$condF] - $lvc; # integrate gradient across runs |
|
313 |
# print(STDERR "ldc/dc/rdc: $last_dcond/$dcond/$run_dcond\n"); |
|
314 |
if ($last_dcond*$dcond >= 0) { # run is continuing (gradient does not change sign) |
|
315 |
$run_dcond += $dcond; |
|
316 |
} else { # run has ended |
|
317 |
# print(STDERR "new run at $lvcs (run_dcond = $last_dcond)\n"); |
|
318 |
$run_start = $lvcs; |
|
319 |
$run_dcond = $dcond; |
|
320 |
} |
|
321 |
||
322 |
if (abs($run_dcond) <= 0.02) { # small integrated gradient => okay |
|
323 |
$lvc = $ants_[$s][$condF]; # update stored previous values |
|
324 |
$lvcs = $s; |
|
325 |
$last_dcond = $dcond; |
|
326 |
next; # process next scan |
|
327 |
} |
|
328 |
||
329 |
# print(STDERR "large gradient ($ants_[$run_start][$condF]-$ants_[$s][$condF], $run_dcond) run $run_start-$s\n"); |
|
330 |
my($i); |
|
331 |
my($max_spike_length) = 24; |
|
332 |
for ($i=$s; $i<=$run_start+$max_spike_length && $i<$nscans-1; $i++) { # large gradient => check whether values return within 10s |
|
333 |
next unless (numberp($ants_[$i][$condF]) && |
|
334 |
numberp($ants_[$i+1][$condF])); |
|
335 |
last if ((abs($ants_[$i][$condF]-$ants_[$run_start][$condF]) <= 0.005) && # 2 vals to avoid large gradients |
|
336 |
(abs($ants_[$i+1][$condF]-$ants_[$run_start][$condF]) <= 0.005)); |
|
337 |
} |
|
338 |
||
339 |
if ($i>$run_start+$max_spike_length || $i==$nscans-1) { # values don't return => leave data alone |
|
340 |
# print(STDERR "values don't return => new run at $s\n"); |
|
341 |
$run_start = $s; $run_dcond = 0; |
|
342 |
$lvc = $ants_[$run_start][$condF]; # start new run |
|
343 |
$lvcs = $run_start; |
|
344 |
$last_dcond = 0; |
|
345 |
next; # process next scan |
|
346 |
} |
|
347 |
||
348 |
$ns++; |
|
349 |
for (my($j)=$run_start+1; $j<$i; $j++) { # values return => remove bad spike |
|
350 |
$ants_[$j][$condF] = nan; |
|
351 |
$nv++; |
|
352 |
} |
|
353 |
# print(STDERR "values return at $i ($ants_[$i][$condF]) => deleted $run_start+1-$i; new run\n"); |
|
354 |
$run_start = $i; $run_dcond = 0; |
|
355 |
$lvc = $ants_[$run_start][$condF]; |
|
356 |
$lvcs = $run_start; |
|
357 |
$last_dcond = 0; |
|
358 |
} |
|
359 |
&antsAddParams("conductivity_spikes_removed",$ns); |
|
360 |
printf(STDERR "\n\t%d conductivity values removed from %d spikes",$nv,$ns) if ($opt_v > 1); |
|
361 |
||
18 | 362 |
printf(STDERR "\n") if ($opt_v); |
363 |
||
364 |
} # if $opt_r |
|
365 |
||
366 |
#---------------------------------------------------------------------- |
|
367 |
# Correcting for pressure bias |
|
368 |
#---------------------------------------------------------------------- |
|
369 |
||
370 |
print(STDERR "Correcting for pressure bias...") if ($opt_v); |
|
371 |
||
372 |
my($minP) = 9e99; |
|
373 |
for (my($s)=0; $s<$nscans; $s++) { |
|
374 |
$minP = $ants_[$s][$pressF] |
|
375 |
if numberp($ants_[$s][$pressF]) && ($ants_[$s][$pressF] < $minP); |
|
376 |
} |
|
377 |
croak("$CNVfile: no valid CTD pressure data below 25dbar\n") |
|
378 |
unless ($minP < 9e99); |
|
29 | 379 |
|
380 |
if ($minP < 25) { |
|
381 |
&antsAddParams('pressure_bias',$minP); |
|
382 |
printf(STDERR "\n\tsubtracting %.1f dbar",$minP) if ($opt_v > 1); |
|
383 |
for (my($s)=0; $s<$nscans; $s++) { |
|
384 |
$ants_[$s][$pressF] -= $minP |
|
385 |
if numberp($ants_[$s][$pressF]); |
|
386 |
} |
|
387 |
} else { |
|
388 |
printf(STDERR "\n\tpartial-depth cast below %.1f dbar (no correction applied)",$minP) if ($opt_v > 1); |
|
18 | 389 |
} |
390 |
printf(STDERR "\n") if ($opt_v); |
|
391 |
||
392 |
#---------------------------------------------------------------------- |
|
393 |
# Binning data |
|
394 |
#---------------------------------------------------------------------- |
|
395 |
||
396 |
my($sps) = round(1 / $sampint / $opt_s); |
|
397 |
print(STDERR "Creating ${opt_s}Hz time series ($sps samples per bin)...") if ($opt_v); |
|
398 |
&antsAddParams('sampling_frequency',1/$opt_s); |
|
399 |
&antsAddParams('sampling_rate',$opt_s); |
|
400 |
||
401 |
my(@press,@temp,@cond); |
|
402 |
my($sp,$np,$st,$nt,$sc,$nc); |
|
403 |
||
404 |
$sp = $st = $sc = $np = $nt = $nc = 0; |
|
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
405 |
$slat = $slon = $nGPS = 0; |
18 | 406 |
for (my($rec)=1,my($s)=0; $s<$nscans; $s++) { |
407 |
if ($s*$sampint > $rec/$opt_s) { |
|
408 |
$rec++; |
|
409 |
push(@press,$np>0?$sp/$np:nan); |
|
410 |
push(@temp, $nt>0?$st/$nt:nan); |
|
411 |
push(@cond, $nc>0?$sc/$nc:nan); |
|
412 |
$sp = $st = $sc = $np = $nt = $nc = 0; |
|
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
413 |
if (defined($latF)) { |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
414 |
push(@lat,$nGPS>0?$slat/$nGPS:nan); |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
415 |
push(@lon,$nGPS>0?$slon/$nGPS:nan); |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
416 |
$slat = $slon = $nGPS = 0; |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
417 |
} |
18 | 418 |
} |
419 |
$sp+=$ants_[$s][$pressF],$np++ if numberp($ants_[$s][$pressF]); |
|
420 |
$st+=$ants_[$s][$tempF],$nt++ if numberp($ants_[$s][$tempF]); |
|
421 |
$sc+=$ants_[$s][$condF],$nc++ if numberp($ants_[$s][$condF]); |
|
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
422 |
if (defined($latF) && numberp($ants_[$s][$latF])) { |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
423 |
$slat += $ants_[$s][$latF]; |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
424 |
$slon += $ants_[$s][$lonF]; |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
425 |
$nGPS++; |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
426 |
} |
18 | 427 |
} |
428 |
||
429 |
printf(STDERR "\n") if ($opt_v); |
|
430 |
||
431 |
#---------------------------------------------------------------------- |
|
432 |
# Calculating derived quantities |
|
433 |
#---------------------------------------------------------------------- |
|
434 |
||
435 |
print(STDERR "Calculating vertical package velocity & sound speed...") if ($opt_v); |
|
436 |
||
29 | 437 |
my($maxP) = -9e99; my($atBtm); |
438 |
my($min_sspd) = 9e99; my($max_sspd) = -9e99; |
|
18 | 439 |
for (my($r)=0; $r<@press; $r++) { |
29 | 440 |
$maxP=$press[$r],$atBtm=$r if ($press[$r] > $maxP); |
18 | 441 |
$elapsed[$r] = $r/$opt_s; |
442 |
$depth[$r] = &depth($press[$r],$lat); |
|
29 | 443 |
$salin[$r] = &salin($cond[$r],$temp[$r],$press[$r]); |
444 |
$sspd[$r] = &sVel($salin[$r],$temp[$r],$press[$r]); |
|
445 |
$min_sspd = $sspd[$r] if ($sspd[$r] < $min_sspd); |
|
446 |
$max_sspd = $sspd[$r] if ($sspd[$r] > $max_sspd); |
|
18 | 447 |
} |
448 |
||
449 |
$w[0] = nan; |
|
450 |
for (my($r)=1; $r<@depth-1; $r++) { |
|
451 |
$w[$r] = numbersp($depth[$r-1],$depth[$r+1]) |
|
452 |
? ($depth[$r+1] - $depth[$r-1]) * $opt_s |
|
453 |
: nan; |
|
454 |
} |
|
455 |
push(@w,nan); |
|
456 |
||
457 |
printf(STDERR "\n") if ($opt_v); |
|
458 |
||
459 |
#---------------------------------------------------------------------- |
|
460 |
# Low-pass filter velocity data |
|
461 |
# - interpolate missing vertical velocities first |
|
462 |
#---------------------------------------------------------------------- |
|
463 |
||
464 |
if ($opt_l > 0) { |
|
465 |
print(STDERR "Low-pass filtering vertical package velocity...") if ($opt_v); |
|
466 |
&antsAddParams('w_lowpass_cutoff',$opt_l); |
|
467 |
||
468 |
my($trimmed) = 0; |
|
469 |
shift(@w),shift(@depth),shift(@elapsed),shift(@sspd),$trimmed++ |
|
470 |
until numberp($w[0]); |
|
471 |
my($interpolated) = 0; |
|
472 |
for ($r=1; $r<@w; $r++) { |
|
473 |
next if numberp($w[$r]); |
|
474 |
my($lv) = $r-1; |
|
475 |
for ($nv=$r+1; $nv<@depth && !numberp($w[$nv]); $nv++) {} |
|
476 |
if ($nv < @depth) { |
|
477 |
while ($r < $nv) { |
|
478 |
$w[$r] = $w[$lv] + ($r-$lv)/($nv-$lv) * ($w[$nv]-$w[$lv]); |
|
479 |
$interpolated++; |
|
480 |
$r++; |
|
481 |
} |
|
482 |
||
483 |
} else { |
|
484 |
$trimmed += @w-$r; |
|
485 |
splice(@w,$r); splice(@depth,$r); |
|
486 |
splice(@elapsed,$r); splice(@sspd,$r); |
|
487 |
} |
|
488 |
} |
|
489 |
&antsAddParams('w_interpolated',$interpolated); |
|
490 |
printf(STDERR "\n\t%d/%d vertical velocities trimmed/interpolated",$trimmed,$interpolated) if ($opt_v > 1); |
|
491 |
||
492 |
||
493 |
#-------------------- |
|
494 |
# Zero Pad Data |
|
495 |
#-------------------- |
|
496 |
||
497 |
for ($pot=1; $pot<@w; $pot<<=1) {} # determine power of two |
|
498 |
||
499 |
for ($r=0; $r<@w; $r++) { # copy data |
|
500 |
$fftbuf[2*$r] = $w[$r]; |
|
501 |
$fftbuf[2*$r+1] = 0; |
|
502 |
} |
|
503 |
printf(STDERR "\n\t%d zero records added",$pot-$r) if ($opt_v > 1); |
|
504 |
while ($r < $pot) { # pad with zeroes |
|
505 |
$fftbuf[2*$r] = $fftbuf[2*$r+1] = 0; |
|
506 |
$r++; |
|
507 |
} |
|
508 |
||
509 |
#-------------------- |
|
510 |
# Low-Pass Filter |
|
511 |
#-------------------- |
|
512 |
||
513 |
@fco = &FOUR1(-1,@fftbuf); # forward FFT |
|
514 |
$n = @fco/2; |
|
515 |
for (my($ip)=2; $ip<=$n; $ip+=2) { # +ve freq fco |
|
516 |
my($in) = 2*$n-$ip; # -ve freq fco |
|
517 |
my($f) = $ip/2/$n*$opt_s; # frequency |
|
518 |
$fco[$ip] = $fco[$ip+1] = $fco[$in] = $fco[$in+1] = 0 |
|
519 |
if ($f > 1/$opt_l); # low-pass filter |
|
520 |
} |
|
521 |
@w_lp = &FOUR1(1,@fco); # inverse FFT |
|
522 |
||
523 |
printf(STDERR "\n") if ($opt_v); |
|
524 |
} else { |
|
525 |
@w_lp = @w; |
|
526 |
} |
|
527 |
||
528 |
#---------------------------------------------------------------------- |
|
529 |
||
29 | 530 |
if (defined($opt_p)) { |
531 |
print(STDERR "Plotting data...\n") if ($opt_v); |
|
532 |
my(@pfmt) = split(',',$opt_p); |
|
533 |
croak("$0: cannot decode -p $opt_p\n") |
|
534 |
unless (@pfmt == 2); |
|
535 |
||
536 |
my($xmin) = $elapsed[0]/60; |
|
537 |
my($xmax) = $elapsed[$#elapsed]/60; |
|
538 |
my($ymin) = -3; my($ymax) = 3; |
|
539 |
my($plotsize) = '13c'; |
|
31 | 540 |
GMT_begin(sprintf($pfmt[0],$id),"-JX${plotsize}","-R$xmin/$xmax/$ymin/$ymax",'-X6 -Y4 -P'); |
29 | 541 |
GMT_psxy('-W2,255/127/80'); |
542 |
for ($r=0; $r<@w; $r++) { |
|
543 |
printf(GMT "%f %f\n",$elapsed[$r]/60,$w_lp[2*$r]/@w_lp); |
|
544 |
GMT_psxy('-W2,46/139/87') if ($r == $atBtm); |
|
545 |
} |
|
546 |
GMT_psbasemap('-Bg60a30f5:"Elapsed Time [min]":/g1a1f0.1:"Vertical Package Velocity [ms@+-1@+]":WeSn'); |
|
547 |
GMT_unitcoords(); |
|
31 | 548 |
GMT_pstext('-G255/127/80'); print(GMT "0.98 0.98 14 0 0 TR dc\n"); |
549 |
GMT_pstext('-G46/139/87'); print(GMT "0.98 0.94 14 0 0 TR uc\n"); |
|
550 |
GMT_pstext('-Gblue'); |
|
551 |
if (defined($outfile)) { printf(GMT "0.02 0.02 14 0 0 BL $outfile\n",$id); } |
|
552 |
else { printf(GMT "0.02 0.02 14 0 0 BL %03d\n",$id); } |
|
29 | 553 |
GMT_end(); |
554 |
||
555 |
my($xmin) = round($min_sspd-3,5); |
|
556 |
my($xmax) = round($max_sspd+3,5); |
|
557 |
my($ymin) = 0; my($ymax) = round($depth[$atBtm]+70,100); |
|
558 |
my($plotsize) = '13c'; |
|
31 | 559 |
GMT_begin(sprintf($pfmt[1],$id),"-JX${plotsize}/-${plotsize}","-R$xmin/$xmax/$ymin/$ymax",'-X6 -Y4 -P'); |
29 | 560 |
GMT_psbasemap('-Bg10a10f1:"Speed of Sound [m/s]":/g1000a500f100:"Depth [m]":WeSn'); |
561 |
GMT_psxy('-W8,255/127/80'); |
|
562 |
for ($r=0; $r<@w; $r++) { |
|
563 |
printf(GMT "%f %f\n",$sspd[$r],$depth[$r]); |
|
564 |
GMT_psxy('-W6,46/139/87') if ($r == $atBtm); |
|
565 |
} |
|
566 |
GMT_unitcoords(); |
|
31 | 567 |
GMT_pstext('-G255/127/80'); print(GMT "0.98 0.02 14 0 0 TR dc\n"); |
568 |
GMT_pstext('-G46/139/87'); print(GMT "0.98 0.06 14 0 0 TR uc\n"); |
|
569 |
GMT_pstext('-Gblue'); |
|
570 |
if (defined($outfile)) { printf(GMT "0.02 0.98 14 0 0 BL $outfile\n",$id); } |
|
571 |
else { printf(GMT "0.02 0.98 14 0 0 BL %03d\n",$id); } |
|
29 | 572 |
GMT_end(); |
573 |
} |
|
574 |
||
575 |
#---------------------------------------------------------------------- |
|
576 |
||
18 | 577 |
print(STDERR "Writing output...\n") if ($opt_v); |
578 |
||
29 | 579 |
@antsNewLayout = ('elapsed','press','temp','cond','depth','salin','sspd','w.raw','w'); |
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
580 |
push(@antsNewLayout,'lat','lon') if (defined($latF)); |
18 | 581 |
for ($r=0; $r<@w; $r++) { |
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
582 |
&antsOut($elapsed[$r],$press[$r],$temp[$r],$cond[$r],$depth[$r],$salin[$r], |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
583 |
$sspd[$r],$w[$r],$w_lp[2*$r]/@w_lp,$lat[$r],$lon[$r]); |
18 | 584 |
} |
585 |
||
586 |
exit(0); # don't flush @ants_ |