1 #!/usr/bin/perl |
1 #!/usr/bin/perl |
2 #====================================================================== |
2 #====================================================================== |
3 # M E A N P R O F |
3 # M E A N P R O F |
4 # doc: Fri Feb 22 08:40:18 2008 |
4 # doc: Fri Feb 22 08:40:18 2008 |
5 # dlm: Wed Mar 16 07:02:38 2016 |
5 # dlm: Tue Apr 10 09:01:50 2018 |
6 # (c) 2008 A.M. Thurnherr |
6 # (c) 2008 A.M. Thurnherr |
7 # uE-Info: 14 49 NIL 0 0 72 0 2 4 NIL ofnI |
7 # uE-Info: 19 35 NIL 0 0 72 0 2 4 NIL ofnI |
8 #====================================================================== |
8 #====================================================================== |
9 |
9 |
10 # extract time-averaged mean profile from ADCP data |
10 # extract time-averaged mean profile from ADCP data |
11 |
11 |
12 # HISTORY: |
12 # HISTORY: |
13 # Feb 22, 2008: - created from [listBins] |
13 # Feb 22, 2008: - created from [listBins] |
14 # Mar 16, 2016: - adapted to new Getopt library |
14 # Mar 16, 2016: - adapted to new Getopt library |
|
15 # Apr 2, 2018: - BUG: velBeamToInstrument() used old usage |
|
16 # Apr 9, 2018: - adapted to "new" readData() ensemble limits |
|
17 # - added -l to set final bin |
|
18 # - BUG: division by zero in empty bins |
|
19 # Apr 10, 2018: - activate output |
15 |
20 |
16 # Soundspeed Correction: |
21 # Soundspeed Correction: |
|
22 # - based on first ensemble only |
17 # - applied as described in the RDI coord-trans manual |
23 # - applied as described in the RDI coord-trans manual |
18 # - sound-speed variation over range is ignored (valid for small gradients) |
24 # - sound-speed variation over range is ignored (valid for small gradients) |
19 # => - same simple correction for all velocity components |
25 # => - same simple correction for all velocity components |
20 # - simple correction for cell depths |
26 # - simple correction for cell depths |
21 |
27 |
22 use Getopt::Std; |
28 use Getopt::Std; |
23 |
29 |
24 $0 =~ m{(.*/)[^/]+}; |
30 $ADCP_tools_minVersion = 2.2; |
25 require "$1RDI_BB_Read.pl"; |
31 ($ADCP_TOOLS) = ($0 =~ m{(.*/)[^/]+}); |
26 require "$1RDI_Coords.pl"; |
32 require "$ADCP_TOOLS/ADCP_tools_lib.pl"; |
27 require "$1RDI_Utils.pl"; |
33 |
28 |
34 die("Usage: $0 [-r)ange <first_ens,last_ens>] [-l)ast <bin>] " . |
29 die("Usage: $0 [-r)ange <first_ens,last_ens>] " . |
|
30 "[-Q)uiet (stats only)] " . |
35 "[-Q)uiet (stats only)] " . |
31 "[-S)oundspeed correction <salin|*,temp|*,depth|*> " . |
36 "[-S)oundspeed correction <salin|*,temp|*,depth|*> " . |
32 "[require -4)-beam solutions] [-d)iscard <beam#>] " . |
37 "[require -4)-beam solutions] [-d)iscard <beam#>] " . |
33 "[-%)good <min>] " . |
38 "[-%)good <min>] " . |
34 "[output -b)eam coordinates] " . |
39 "[output -b)eam coordinates] " . |
35 "[-M)agnetic <declination>] " . |
40 "[-M)agnetic <declination>] " . |
36 "[-D)epth <depth>] " . |
41 "[-D)epth <depth>] " . |
37 "<RDI file>\n") |
42 "<RDI file>\n") |
38 unless (&getopts("4bd:D:M:p:r:QS:") && @ARGV == 1); |
43 unless (&getopts("4bd:D:l:M:p:r:QS:") && @ARGV == 1); |
39 |
44 |
40 die("$0: -4 and -d are mutually exclusive\n") |
45 die("$0: -4 and -d are mutually exclusive\n") |
41 if ($opt_4 && defined($opt_d)); |
46 if ($opt_4 && defined($opt_d)); |
42 |
47 |
43 die("$0: -p and -b are mutually exclusive\n") |
48 die("$0: -p and -b are mutually exclusive\n") |
49 |
54 |
50 print(STDERR "WARNING: magnetic declination not set!\n") |
55 print(STDERR "WARNING: magnetic declination not set!\n") |
51 unless defined($opt_M) || defined($opt_b); |
56 unless defined($opt_M) || defined($opt_b); |
52 |
57 |
53 $ifn = $ARGV[0]; |
58 $ifn = $ARGV[0]; |
54 |
|
55 ($first_ens,$last_ens) = split(',',$opt_r) |
|
56 if defined($opt_r); |
|
57 |
59 |
58 ($SS_salin,$SS_temp,$SS_depth) = split(',',$opt_S) |
60 ($SS_salin,$SS_temp,$SS_depth) = split(',',$opt_S) |
59 if defined($opt_S); |
61 if defined($opt_S); |
60 die("$0: Cannot do variable soundspeed correction (implementation restriction)\n") |
62 die("$0: Cannot do variable soundspeed correction (implementation restriction)\n") |
61 if ($SS_salin eq '*' || $SS_temp eq '*' || $SS_depth eq '*'); |
63 if ($SS_salin eq '*' || $SS_temp eq '*' || $SS_depth eq '*'); |
66 |
68 |
67 $P{RDI_file} = $ifn; |
69 $P{RDI_file} = $ifn; |
68 $P{mag_decl} = $opt_M if defined($opt_M); |
70 $P{mag_decl} = $opt_M if defined($opt_M); |
69 |
71 |
70 print(STDERR "reading $ifn: "); |
72 print(STDERR "reading $ifn: "); |
71 readData($ifn,\%dta); |
73 if (defined($opt_r)) { # read selected range |
|
74 my($fe,$le) = split(',',$opt_r); |
|
75 readData($ifn,\%dta,$fe,$le,$opt_l); |
|
76 } else { # read entire file |
|
77 readData($ifn,\%dta,undef,undef,$opt_l); |
|
78 } |
72 printf(STDERR "%d complete ensembles.\n",scalar(@{$dta{ENSEMBLE}})); |
79 printf(STDERR "%d complete ensembles.\n",scalar(@{$dta{ENSEMBLE}})); |
73 $dta{HEADING_BIAS} = -$opt_M; # magnetic declination |
80 $dta{HEADING_BIAS} = -$opt_M; # magnetic declination |
74 |
81 |
75 if ($dta{BEAM_COORDINATES}) { # coords |
82 if ($dta{BEAM_COORDINATES}) { # coords |
76 $beamCoords = 1; |
83 $beamCoords = 1; |
85 $dz[$b] = $dta{DISTANCE_TO_BIN1_CENTER} + $b*$dta{BIN_LENGTH}; |
92 $dz[$b] = $dta{DISTANCE_TO_BIN1_CENTER} + $b*$dta{BIN_LENGTH}; |
86 } |
93 } |
87 |
94 |
88 $lastGoodBin = 0; |
95 $lastGoodBin = 0; |
89 for ($e=0; $e<=$#{$dta{ENSEMBLE}}; $e++) { # check/transform velocities |
96 for ($e=0; $e<=$#{$dta{ENSEMBLE}}; $e++) { # check/transform velocities |
90 next if (defined($first_ens) && |
|
91 $dta{ENSEMBLE}[$e]->{NUMBER} < $first_ens); |
|
92 $P{first_ens} = $dta{ENSEMBLE}[$e]->{NUMBER},$fe = $e |
97 $P{first_ens} = $dta{ENSEMBLE}[$e]->{NUMBER},$fe = $e |
93 unless defined($P{first_ens}); |
98 unless defined($P{first_ens}); |
94 last if (defined($last_ens) && |
|
95 $dta{ENSEMBLE}[$e]->{NUMBER} > $last_ens); |
|
96 $P{last_ens} = $dta{ENSEMBLE}[$e]->{NUMBER}; |
99 $P{last_ens} = $dta{ENSEMBLE}[$e]->{NUMBER}; |
97 $le = $e; |
100 $le = $e; |
98 |
101 |
99 die("3-beams used in ensemble #$dta{ENSEMBLE}[$e]->{NUMBER}\n") |
102 die("3-beams used in ensemble #$dta{ENSEMBLE}[$e]->{NUMBER}\n") |
100 if ($dta{ENSEMBLE}[$e]->{N_BEAMS_USED} < 4); |
103 if ($dta{ENSEMBLE}[$e]->{N_BEAMS_USED} < 4); |
114 undef($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][$i]); |
117 undef($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][$i]); |
115 } |
118 } |
116 } |
119 } |
117 @{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]} = $beamCoords |
120 @{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]} = $beamCoords |
118 ? velInstrumentToEarth(\%dta,$e, |
121 ? velInstrumentToEarth(\%dta,$e, |
119 velBeamToInstrument(\%dta,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]}) |
122 velBeamToInstrument(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]}) |
120 ) |
123 ) |
121 : velApplyHdgBias(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]}) |
124 : velApplyHdgBias(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]}) |
122 unless ($opt_b); |
125 unless ($opt_b); |
123 |
126 |
124 $sum_corr1[$b] += $dta{ENSEMBLE}[$e]->{CORRELATION}[$b][0]; |
127 $sum_corr1[$b] += $dta{ENSEMBLE}[$e]->{CORRELATION}[$b][0]; |
188 for ($b=0; $b<=$lastGoodBin; $b++) { |
191 for ($b=0; $b<=$lastGoodBin; $b++) { |
189 $mean_corr1[$b] = $sum_corr1[$b] / $nEns; $mean_corr2[$b] = $sum_corr2[$b] / $nEns; |
192 $mean_corr1[$b] = $sum_corr1[$b] / $nEns; $mean_corr2[$b] = $sum_corr2[$b] / $nEns; |
190 $mean_corr3[$b] = $sum_corr3[$b] / $nEns; $mean_corr4[$b] = $sum_corr4[$b] / $nEns; |
193 $mean_corr3[$b] = $sum_corr3[$b] / $nEns; $mean_corr4[$b] = $sum_corr4[$b] / $nEns; |
191 $mean_amp1[$b] = $sum_amp1[$b] / $nEns; $mean_amp2[$b] = $sum_amp2[$b] / $nEns; |
194 $mean_amp1[$b] = $sum_amp1[$b] / $nEns; $mean_amp2[$b] = $sum_amp2[$b] / $nEns; |
192 $mean_amp3[$b] = $sum_amp3[$b] / $nEns; $mean_amp4[$b] = $sum_amp4[$b] / $nEns; |
195 $mean_amp3[$b] = $sum_amp3[$b] / $nEns; $mean_amp4[$b] = $sum_amp4[$b] / $nEns; |
193 $mean_pcg1[$b] = $sum_pcg1[$b] / $n_pcg1[$b]; $mean_pcg2[$b] = $sum_pcg2[$b] / $n_pcg2[$b]; |
196 |
194 $mean_pcg3[$b] = $sum_pcg3[$b] / $n_pcg3[$b]; $mean_pcg4[$b] = $sum_pcg4[$b] / $n_pcg4[$b]; |
197 $mean_pcg1[$b] = $sum_pcg1[$b] / $n_pcg1[$b] if ($n_pcg1[$b] > 0); |
195 $mean_u[$b] = $sum_u[$b] / $good_vels[$b]; $mean_v[$b] = $sum_v[$b] / $good_vels[$b]; |
198 $mean_pcg2[$b] = $sum_pcg2[$b] / $n_pcg2[$b] if ($n_pcg2[$b] > 0); |
|
199 $mean_pcg3[$b] = $sum_pcg3[$b] / $n_pcg3[$b] if ($n_pcg3[$b] > 0); |
|
200 $mean_pcg4[$b] = $sum_pcg4[$b] / $n_pcg4[$b] if ($n_pcg4[$b] > 0); |
|
201 |
|
202 next unless ($good_vels[$b] > 0); |
|
203 |
|
204 $mean_u[$b] = $sum_u[$b] / $good_vels[$b]; |
|
205 $mean_v[$b] = $sum_v[$b] / $good_vels[$b]; |
196 $mean_w[$b] = $sum_w[$b] / $good_vels[$b]; |
206 $mean_w[$b] = $sum_w[$b] / $good_vels[$b]; |
197 $mean_e[$b] = $sum_e[$b] / ($good_vels[$b] - $three_beam[$b]); |
207 $mean_e[$b] = ($good_vels[$b] - $three_beam[$b] > 0) ? $sum_e[$b] / ($good_vels[$b] - $three_beam[$b]) : undef; |
198 } |
208 } |
199 |
209 |
200 for ($e=$fe; $e<=$le; $e++) { |
210 for ($e=$fe; $e<=$le; $e++) { |
201 for ($b=0; $b<=$lastGoodBin; $b++) { |
211 for ($b=0; $b<=$lastGoodBin; $b++) { |
202 $sumsq_corr1[$b] += ($mean_corr1[$b] - $dta{ENSEMBLE}[$e]->{CORRELATION}[$b][0])**2; |
212 $sumsq_corr1[$b] += ($mean_corr1[$b] - $dta{ENSEMBLE}[$e]->{CORRELATION}[$b][0])**2; |
208 $sumsq_amp2[$b] += ($mean_amp2[$b] - $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b][1])**2; |
218 $sumsq_amp2[$b] += ($mean_amp2[$b] - $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b][1])**2; |
209 $sumsq_amp3[$b] += ($mean_amp3[$b] - $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b][2])**2; |
219 $sumsq_amp3[$b] += ($mean_amp3[$b] - $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b][2])**2; |
210 $sumsq_amp4[$b] += ($mean_amp4[$b] - $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b][3])**2; |
220 $sumsq_amp4[$b] += ($mean_amp4[$b] - $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b][3])**2; |
211 |
221 |
212 $sumsq_pcg1[$b] += ($mean_pcg1[$b] - $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][0])**2 |
222 $sumsq_pcg1[$b] += ($mean_pcg1[$b] - $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][0])**2 |
213 if defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][0]); |
223 if defined($mean_pcg1[$b]) && defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][0]); |
214 $sumsq_pcg2[$b] += ($mean_pcg2[$b] - $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][1])**2 |
224 $sumsq_pcg2[$b] += ($mean_pcg2[$b] - $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][1])**2 |
215 if defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][1]); |
225 if defined($mean_pcg2[$b]) && defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][1]); |
216 $sumsq_pcg3[$b] += ($mean_pcg3[$b] - $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][2])**2 |
226 $sumsq_pcg3[$b] += ($mean_pcg3[$b] - $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][2])**2 |
217 if defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][2]); |
227 if defined($mean_pcg3[$b]) && defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][2]); |
218 $sumsq_pcg4[$b] += ($mean_pcg4[$b] - $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][3])**2 |
228 $sumsq_pcg4[$b] += ($mean_pcg4[$b] - $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][3])**2 |
219 if defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][3]); |
229 if defined($mean_pcg4[$b]) && defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][3]); |
220 |
230 |
221 next unless ($dta{ENSEMBLE}[$e]->{GOOD_VEL}[$b]); |
231 next unless ($dta{ENSEMBLE}[$e]->{GOOD_VEL}[$b]); |
222 |
232 |
223 $sumsq_u[$b] += ($mean_u[$b] - $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0])**2; |
233 $sumsq_u[$b] += ($mean_u[$b] - $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0])**2; |
224 $sumsq_v[$b] += ($mean_v[$b] - $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][1])**2; |
234 $sumsq_v[$b] += ($mean_v[$b] - $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][1])**2; |
225 $sumsq_w[$b] += ($mean_w[$b] - $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][2])**2; |
235 $sumsq_w[$b] += ($mean_w[$b] - $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][2])**2; |
226 |
236 |
227 next if ($dta{ENSEMBLE}[$e]->{THREE_BEAM}[$b]); |
237 next if ($dta{ENSEMBLE}[$e]->{THREE_BEAM}[$b]); |
228 |
238 |
229 $sumsq_e[$b] += ($mean_e[$b] - $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3])**2; |
239 $sumsq_e[$b] += ($mean_e[$b] - $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3])**2 |
|
240 if defined($mean_e[$b]); |
230 } |
241 } |
231 } |
242 } |
232 |
243 |
233 for ($b=0; $b<=$lastGoodBin; $b++) { |
244 for ($b=0; $b<=$lastGoodBin; $b++) { |
234 $var_corr1[$b] = $sumsq_corr1[$b] / ($nEns-1); $var_corr2[$b] = $sumsq_corr2[$b] / ($nEns-1); |
245 $var_corr1[$b] = $sumsq_corr1[$b] / ($nEns-1); $var_corr2[$b] = $sumsq_corr2[$b] / ($nEns-1); |
249 |
260 |
250 #---------------------------------------------------------------------- |
261 #---------------------------------------------------------------------- |
251 # Calculate Beam Statistics |
262 # Calculate Beam Statistics |
252 #---------------------------------------------------------------------- |
263 #---------------------------------------------------------------------- |
253 |
264 |
|
265 # not implemented yet |
|
266 |
254 #---------------------------------------------------------------------- |
267 #---------------------------------------------------------------------- |
255 # Produce Output |
268 # Produce Output |
256 #---------------------------------------------------------------------- |
269 #---------------------------------------------------------------------- |
257 |
270 |
258 unless ($opt_Q) { |
271 unless ($opt_Q) { |
259 my($ssCorr) = defined($opt_S) |
272 my($ssCorr) = defined($opt_S) |
260 ? ssCorr($dta{ENSEMBLE}[$fe],$SS_salin,$SS_temp,$SS_depth) |
273 ? ssCorr($dta{ENSEMBLE}[$fe],$SS_salin,$SS_temp,$SS_depth) |
261 : 1; |
274 : 1; |
262 |
275 |
|
276 print("#!/usr/bin/perl -S list\n"); |
|
277 chmod(0777&~umask,*STDOUT); |
263 print("#ANTS#PARAMS# "); |
278 print("#ANTS#PARAMS# "); |
264 foreach my $k (keys(%P)) { |
279 foreach my $k (keys(%P)) { |
265 print("$k\{$P{$k}\} "); |
280 print("$k\{$P{$k}\} "); |
266 } |
281 } |
267 printf("soundspeed_correction{%s}",defined($opt_S) ? $opt_S : 'NONE!'); |
282 printf("soundspeed_correction{%s}",defined($opt_S) ? $opt_S : 'NONE!'); |