1 #====================================================================== |
|
2 # L A D C P P R O C . B E S T L A G |
|
3 # doc: Tue Sep 28 21:58:48 2010 |
|
4 # dlm: Fri Jul 8 02:54:29 2011 |
|
5 # (c) 2010 A.M. Thurnherr |
|
6 # uE-Info: 55 101 NIL 0 0 72 2 2 4 NIL ofnI |
|
7 #====================================================================== |
|
8 |
|
9 # TODO: |
|
10 # - first lag is always(?) nan unless CTD is turned on with LADCP in water |
|
11 |
|
12 # HISTORY: |
|
13 # Sep 28, 2010: - created |
|
14 # Dec 9, 2010: - adapted to %CTD |
|
15 # Dec 10, 2010: - hardened bestlag failure test to require 1/3 agreeing lags |
|
16 # Jan 5, 2011: - changed first guess from 80% down to 10% down |
|
17 # - added LADCP time lag to %PARAMs |
|
18 # - added support of -i |
|
19 # Jul 7, 2011: - added code to remove window-mean of w before lagging to |
|
20 # make it work in regions of crazy ocean w (IWISE 16007) |
|
21 |
|
22 sub interp_LADCP_w($$) |
|
23 { |
|
24 my($elapsed,$ens) = @_; |
|
25 my($sc) = ($elapsed - $LADCP{ENSEMBLE}[$ens-1]->{ELAPSED_TIME}) / |
|
26 ($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} - |
|
27 $LADCP{ENSEMBLE}[$ens-1]->{ELAPSED_TIME}); |
|
28 unless (numberp($LADCP{ENSEMBLE}[$ens-1]->{W})) { |
|
29 $nGaps++; |
|
30 return $LADCP{ENSEMBLE}[$ens]->{W}; |
|
31 } |
|
32 return $LADCP{ENSEMBLE}[$ens-1]->{W} + |
|
33 $sc * ($LADCP{ENSEMBLE}[$ens]->{W} - $LADCP{ENSEMBLE}[$ens-1]->{W}); |
|
34 } |
|
35 |
|
36 sub bestLag($) |
|
37 { |
|
38 my($ws) = @_; # window start index |
|
39 |
|
40 my($best); |
|
41 my($bestmad) = 9e99; # mean absolute deviation |
|
42 for (my($Llag)=-int($opt_w/2); $Llag<int($opt_w/2); $Llag++) { |
|
43 my($mCw,$mLw,$nw) = (0,0,0); # first calc means |
|
44 for (my($Ci)=0; $Ci<$opt_w; $Ci++) { |
|
45 my($Li) = $Ci + $Llag; |
|
46 next if ($Li<0 || $Li>=$opt_w); |
|
47 next unless numberp($CTD{w}[$ws+$Ci]) && numberp($LADCP_w[$ws+$Li]); |
|
48 $mCw += $CTD{w}[$ws+$Ci]; |
|
49 $mLw += $LADCP_w[$ws+$Li]; |
|
50 $nw++; |
|
51 } |
|
52 next unless ($nw > 0); |
|
53 $mCw /= $nw; $mLw /= $nw; |
|
54 |
|
55 my($sad) = my($nad) = 0; # calc mad with means removed |
|
56 for (my($Ci)=0; $Ci<$opt_w; $Ci++) { |
|
57 my($Li) = $Ci + $Llag; |
|
58 next if ($Li<0 || $Li>=$opt_w); |
|
59 next unless numberp($CTD{w}[$ws+$Ci]) && numberp($LADCP_w[$ws+$Li]); |
|
60 $sad += abs($CTD{w}[$ws+$Ci]-$mCw - ($LADCP_w[$ws+$Li]-$mLw)); |
|
61 $nad++; |
|
62 } |
|
63 if ($sad/$nad < $bestmad) { |
|
64 $best = $Llag; |
|
65 $bestmad = $sad/$nad; |
|
66 } |
|
67 } |
|
68 return $best; |
|
69 } |
|
70 |
|
71 sub lagLADCP2CTD() |
|
72 { |
|
73 #------------------------------------------------------------------------ |
|
74 # find 1st rec & ensemble >=10% down to max depth & make 1st guess at lag |
|
75 #------------------------------------------------------------------------ |
|
76 |
|
77 my($first_guess_lag); # in units of CTD records |
|
78 |
|
79 if (defined($opt_i)) { |
|
80 $first_guess_lag = -$opt_i / $CTD{sampint}; |
|
81 } else { |
|
82 my($CTD_10pct_down) = 0; |
|
83 $CTD_10pct_down++ # "until" formulation allows for missing pressures |
|
84 until ($CTD{press}[$CTD_10pct_down]-$CTD{press}[0] >= 0.1*($CTD{maxpress}-$CTD{press}[0])); |
|
85 |
|
86 my($LADCP_10pct_down) = 0; |
|
87 $LADCP_10pct_down++ |
|
88 while ($LADCP{ENSEMBLE}[$LADCP_10pct_down]->{DEPTH} < 0.1*$LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH}); |
|
89 |
|
90 $first_guess_lag = ($LADCP{ENSEMBLE}[$LADCP_10pct_down]->{ELAPSED_TIME} - |
|
91 $CTD_10pct_down*$CTD{sampint}) / $CTD{sampint}; |
|
92 |
|
93 printf(STDERR "\n\t1st guess offset [CTD pressure, LADCP estimated depth] = %ds [%ddbar, %dm]\n", |
|
94 $first_guess_lag*$CTD{sampint},$CTD{press}[$CTD_10pct_down],$LADCP{ENSEMBLE}[$LADCP_10pct_down]->{DEPTH}) |
|
95 if ($opt_d); |
|
96 croak("\n$0: cannot determine valid 1st guess offset ($LADCP_10pct_down,$CTD_10pct_down,$CTD{press}[$CTD_10pct_down])\n") |
|
97 if ($first_guess_lag*$CTD{sampint} > 1200); |
|
98 } |
|
99 |
|
100 #------------------------------------------------------------------------------------ |
|
101 # Linearly interpolate LADCP time series onto a new grid with $CTD{sampint} resolution |
|
102 # ALSO: apply first_guess_lag to make lags small, which keeps the bestlag data |
|
103 # chunks large |
|
104 #------------------------------------------------------------------------------------ |
|
105 |
|
106 $nGaps = 0; |
|
107 |
|
108 for (my($ens)=$LADCP_start,my($r)=0; $ens<=$LADCP_end; $ens++) { |
|
109 while ($r*$CTD{sampint} < $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}) { |
|
110 $LADCP_w[$r-$first_guess_lag] = interp_LADCP_w($r*$CTD{sampint},$ens) |
|
111 unless ($first_guess_lag > $r); |
|
112 $r++; |
|
113 } |
|
114 } |
|
115 |
|
116 print(STDERR "\t$nGaps gaps in w timeseries") |
|
117 if ($opt_d); |
|
118 |
|
119 print(STDERR "\n"); |
|
120 |
|
121 #---------------------------------------------------------------------- |
|
122 # Calculate lags |
|
123 #---------------------------------------------------------------------- |
|
124 |
|
125 printf(STDERR "\tcalculating $opt_n lags from %ds-long windows [s]:",$opt_w); |
|
126 $opt_w = int($opt_w / $CTD{sampint}); |
|
127 |
|
128 #--------------------------------------------------------------- |
|
129 # carry out opt_n lag correlations and keep tally of the results |
|
130 #--------------------------------------------------------------- |
|
131 my(%nBest); |
|
132 my($nLags) = 0; |
|
133 my($lags) = ''; |
|
134 my($lastLag) = 9e99; my($nSame) = 1; |
|
135 for (my($window)=0; $window<$opt_n; $window++) { |
|
136 my($ws) = $window * int(@LADCP_w/$opt_n); # window start |
|
137 $ws = @LADCP_w-$opt_w if ($ws+$opt_w >= @LADCP_w); |
|
138 $bestLag = bestLag($ws); |
|
139 if (defined($bestLag)) { |
|
140 if (defined($lastLag) && $bestLag == $lastLag) { |
|
141 $nSame++; |
|
142 } else { |
|
143 printf(STDERR "(x%d)",$nSame) |
|
144 if ($nSame > 1); |
|
145 printf(STDERR " %d",$bestLag*$CTD{sampint}); |
|
146 $nSame = 1; |
|
147 $lastLag = $bestLag; |
|
148 } |
|
149 $lags .= sprintf(" %s",$bestLag*$CTD{sampint}); |
|
150 $nBest{$bestLag}++; |
|
151 $nLags++; |
|
152 } else { |
|
153 if (!defined($lastLag)) { |
|
154 $nSame++; |
|
155 } else { |
|
156 printf(STDERR "(x%d)",$nSame) |
|
157 if ($nSame > 1); |
|
158 printf(STDERR " nan"); |
|
159 $nSame = 1; |
|
160 $lastLag = $bestLag; |
|
161 } |
|
162 } |
|
163 } |
|
164 printf(STDERR "(x%d)",$nSame) if ($nSame > 1); |
|
165 &antsAddParams('time_lags',$lags); |
|
166 |
|
167 #---------------------- |
|
168 # find most popular lag |
|
169 #---------------------- |
|
170 my($best_lag); |
|
171 foreach my $i (keys(%nBest)) { |
|
172 $best_lag = $i if ($nBest{$i} > $nBest{$best_lag}); |
|
173 } |
|
174 croak("\n$0: cannot determine a valid lag\n") |
|
175 unless ($nBest{$best_lag} > $opt_n/3); |
|
176 print(STDERR "\n\n\t\tWARNING: only $nBest{$best_lag} of the lag estimates agree!\n") |
|
177 if ($nBest{$best_lag} < $opt_n/2); |
|
178 |
|
179 if ($nBest{$best_lag} == $nLags) { |
|
180 printf(STDERR "\n\t\tunanimous lag = %ds\n",$best_lag*$CTD{sampint}); |
|
181 } else { |
|
182 printf(STDERR "\n\t\tmost popular lag = %ds\n",$best_lag*$CTD{sampint}); |
|
183 } |
|
184 |
|
185 &antsAddParams('LADCP_time_lag',($first_guess_lag + $best_lag) * $CTD{sampint}); |
|
186 return ($first_guess_lag + $best_lag) * $CTD{sampint}; |
|
187 } |
|
188 |
|
189 1; |
|