.
authorA.M. Thurnherr <athurnherr@yahoo.com>
Tue, 20 May 2014 09:08:43 +0000
changeset 15 dfcb6bef9d42
parent 14 fea65697bc7b
child 16 29e867b3e070
.
LADCP_w
acoustic_backscatter.pl
bottom_tracking.pl
edit_data.pl
svel_corrections.pl
time_lag.pl
--- a/LADCP_w	Thu Nov 21 10:24:23 2013 +0000
+++ b/LADCP_w	Tue May 20 09:08:43 2014 +0000
@@ -2,9 +2,9 @@
 #======================================================================
 #                    L A D C P _ W 
 #                    doc: Fri Dec 17 18:11:13 2010
-#                    dlm: Thu Sep  5 16:23:16 2013
+#                    dlm: Mon May 19 22:21:35 2014
 #                    (c) 2010 A.M. Thurnherr
-#                    uE-Info: 495 0 NIL 0 0 72 2 2 4 NIL ofnI
+#                    uE-Info: 932 27 NIL 0 0 72 2 2 4 NIL ofnI
 #======================================================================
 
 # TODO:
@@ -136,6 +136,8 @@
 #				  - BUG: $bad_beam did not discard BT_VELOCITY data
 #	Jun  6, 2013: - BUG: error message had -a instead of -d
 #	Sep  5, 2013: - BUG: w12/w34 do not work for earth-coordinate data, of course
+#	Apr 17, 2014: - BUG: edit_tilt was never called when all recorded bins are valid
+#	Apr 21, 2014: - updated comments
 
 # CTD REQUIREMENTS
 #	- elapsed		elapsed seconds; see note below
@@ -589,17 +591,15 @@
 	last;
 }
 
-if (defined($first_bad_bin)) {
-	$fprm = $pte = 0;
-	for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
-		$pte += editTilt($ens,$opt_t);
-		$fprm += editFarBins($ens,$first_bad_bin) if defined($first_bad_bin);
-	}
-	progress("\tattitude threshold (max_tilt = %d deg): %d velocites removed (%d%% of total)\n",
-		$opt_t,$pte,round(100*$pte/$nvv));
-	progress("\tvelocities beyond bin $first_bad_bin (<%d%% valid values): %d velocites removed (%d%% of total in bins $LADCP_firstBin-$LADCP_lastBin)\n",
-		round(100*$per_bin_valid_frac_lim),$fprm,round(100*$fprm/$nvw));
+$fprm = $pte = 0;
+for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
+	$pte += editTilt($ens,$opt_t);
+	$fprm += editFarBins($ens,$first_bad_bin) if defined($first_bad_bin);
 }
+progress("\tattitude threshold (max_tilt = %d deg): %d velocites removed (%d%% of total)\n",
+	$opt_t,$pte,round(100*$pte/$nvv));
+progress("\tvelocities beyond bin $first_bad_bin (<%d%% valid values): %d velocites removed (%d%% of total in bins $LADCP_firstBin-$LADCP_lastBin)\n",
+	round(100*$per_bin_valid_frac_lim),$fprm,round(100*$fprm/$nvw));
 
 #--------------
 # Read CTD data
@@ -658,15 +658,23 @@
 # Construct sound-speed correction profile from CTD 1Hz downcast data
 #	very simple algorithm that stores the last value found
 #	in each 1m bin
+# For PPI filtering, a sound speed profile to the surface is required.
+# 	This is ensured by extrapolating the first value up to zero
 #--------------------------------------------------------------------
 
 progress("Constructing sound-speed correction profile\n");
 
 my($scans_per_sec) = int(1/$CTD{DT}+0.5);
+my($min_depth) = 9e99;
 for (my($s)=0; $s<=$CTD_atbottom; $s+=$scans_per_sec) {
 	next unless ($CTD{DEPTH}[$s] >= 0 && numberp($CTD{SVEL}[$s]));
+	$min_depth = $s if ($s < $min_depth);
 	$sVelProf[int($CTD{DEPTH}[$s])] = $CTD{SVEL}[$s];
 }
+while ($min_depth > 0) {
+	$sVelProf[$min_depth-1] = $sVelProf[$min_depth];
+	$min_depth--;
+}
 
 #-------------------
 # Determine time lag
@@ -855,7 +863,9 @@
 calc_backscatter_profs($firstGoodEns,$lastGoodEns);
 
 #----------------------------------------------------------------------------
-# Remove data contaminated by sidelobe reflection from seabed and sea surface
+# Edit data
+#	1) contaminated by sidelobe reflection from seabed and sea surface
+#	2) PPI 
 #----------------------------------------------------------------------------
 
 if ($LADCP{ENSEMBLE}[$LADCP_atbottom]->{XDUCER_FACING_DOWN}) {
@@ -894,15 +904,37 @@
 		($nvrm,$nerm) = editSideLobes($firstGoodEns,$lastGoodEns,$water_depth);
 		progress("\t$nvrm velocities from $nerm ensembles removed\n");
 
+		progress("Editing data to remove PPI from seabed...\n");
+		  progress("\tConstructing travel-time profile...\n");
+		  my($tt) = ($water_depth - $#sVelProf) / $sVelProf[$#sVelProf];  # $#sVelProf = max_depth(profile) in meters
+		  $ttProf[$#sVelProf] = $tt;
+		  for (my($d)=$#sVelProf-1; $d>=0; $d--) {
+			  $tt += 1 / $sVelProf[$d];
+			  $ttProf[$d] = $tt;
+          }
+		($nvrm,$nerm) = editPPI($firstGoodEns,$lastGoodEns,$water_depth);
+		progress("\t$nvrm velocities from $nerm ensembles removed\n");
 	} else {
-		info("no seabed found in backscatter profiles --- no sidelobe editing done\n");
+		info("no seabed found in backscatter profiles --- cannot edit sidelobe or PPI\n");
 	}
 	
 } else {
 	&antsAddParams('ADCP_orientation','uplooker');
+
 	progress("Editing data to remove sidelobe interference from sea surface...\n");
 	($nvrm,$nerm) = editSideLobes($firstGoodEns,$lastGoodEns,undef);
 	progress("\t$nvrm velocities from $nerm ensembles removed\n");
+
+	progress("Editing data to remove PPI from sea surface...\n");
+	  progress("\tConstructing travel-time profile...\n");
+	  my($tt) = 0;
+	  $ttProf[0] = $tt;
+	  for (my($d)=1; $d<=$#sVelProf; $d++) {
+		  $tt += 1 / $sVelProf[$d];
+		  $ttProf[$d] = $tt;
+      }
+	($nvrm,$nerm) = editPPI($firstGoodEns,$lastGoodEns,undef);
+	progress("\t$nvrm velocities from $nerm ensembles removed\n");
 }
 
 #----------------------------------------------------------------------
--- a/acoustic_backscatter.pl	Thu Nov 21 10:24:23 2013 +0000
+++ b/acoustic_backscatter.pl	Tue May 20 09:08:43 2014 +0000
@@ -1,9 +1,9 @@
 #======================================================================
 #                    A C O U S T I C _ B A C K S C A T T E R . P L 
 #                    doc: Wed Oct 20 13:02:27 2010
-#                    dlm: Fri Oct 21 11:29:13 2011
+#                    dlm: Thu Apr 17 08:49:21 2014
 #                    (c) 2010 A.M. Thurnherr
-#                    uE-Info: 17 0 NIL 0 0 72 2 2 4 NIL ofnI
+#                    uE-Info: 19 34 NIL 0 0 72 0 2 4 NIL ofnI
 #======================================================================
 
 # HISTORY:
@@ -15,6 +15,8 @@
 #				  - BUG: acoustic-backscatter assumed 0 deg C
 #				  - SV now saved in ensemble
 #	Oct 21, 2011: - BUG: made code work for uplooker again
+#	Mar  4, 2014: - added support for missing PITCH/ROLL (TILT)
+#	Apr 17, 2017: - BUG: missing ;
 
 #----------------------------------------------------------------------
 # Volume Scattering Coefficient, following Deines (IEEE 1999)
@@ -76,7 +78,7 @@
 		my(@bd) = calc_binDepths($ens);
 		for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
 			my($depth) = int($bd[$bin]);
-			next if ($depth < 0);
+			next if ($depth<0 || !defined($LADCP{ENSEMBLE}[$ens]->{TILT}));
 			my($range_to_bin) = abs($bd[$bin] - $LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH})
 									/ cos(rad($LADCP{ENSEMBLE}[$ens]->{TILT}))
 									/ $cosBeamAngle;
--- a/bottom_tracking.pl	Thu Nov 21 10:24:23 2013 +0000
+++ b/bottom_tracking.pl	Tue May 20 09:08:43 2014 +0000
@@ -1,9 +1,9 @@
 #======================================================================
 #                    B O T T O M _ T R A C K I N G . P L 
 #                    doc: Wed Oct 20 21:05:37 2010
-#                    dlm: Thu May 16 22:03:15 2013
+#                    dlm: Tue Mar  4 13:54:15 2014
 #                    (c) 2010 A.M. Thurnherr
-#                    uE-Info: 79 111 NIL 0 0 72 2 2 4 NIL ofnI
+#                    uE-Info: 15 43 NIL 0 0 72 2 2 4 NIL ofnI
 #======================================================================
 
 # HISTORY:
@@ -12,6 +12,7 @@
 #	Oct 11, 2011: - moved defaults to [defaults.pl]
 #	Oct 24, 2011: - disabled not-very-useful %BT-params
 #	Apr 22, 2013: - replace output_bin_size by opt_o
+#	Mar  4, 2014: - removed old unused code
 
 # This code is essentially identical to the one used in LADCPproc. Differences:
 #	1) velocity editing is simpler: no wake editing, no PPI editing, no shear
@@ -24,7 +25,6 @@
 #	input:	ensemble number, water depth (with uncertainty)
 #	output:	@BTu,@BTv,@BTw				main result
 #			editflags 					for information
-#			@BTbtmu, @BTbtmv, @BTtilt 	stats to find reasons for quality differences
 #----------------------------------------------------------------------
 
 my($nBTfound,$nBTdepthFlag,$nBTvalidVelFlag,$nBTwFlag) = (0,0,0,0);
@@ -98,11 +98,6 @@
 	$nBTwFlag++,return										# velocity error is too great
 		if (abs($seafloor_w-$LADCP{ENSEMBLE}[$ens]->{REFLR_W}) > $BT_max_w_error);
 
-	push(@BTbtmu,$seafloor_u);								# calc stats to try to find reasons for quality
-	push(@BTbtmv,$seafloor_v);
-	push(@BTbtmw,$seafloor_w);
-	push(@BTtilt,$LADCP{ENSEMBLE}[$ens]->{TILT});
-
 	for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
 		next unless defined($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
 		my($gi) = int($bd[$bin]) / $opt_o;
@@ -135,14 +130,6 @@
 		$BT{MEDIAN_W}[$gi] = median(@{$BTw[$gi]});
         $BT{MAD_W}[$gi] = mad2($BT{W}[$gi],@{$BTw[$gi]});
 	}
-
-#	&antsAddParams('BT_rms_seafloor_u',round(rms(@BTbtmu),0.01),
-#				   'BT_rms_seafloor_v',round(rms(@BTbtmv),0.01),
-#				   'BT_rms_seafloor_w',round(rms(@BTbtmw),0.01),
-#				   'BT_avg_seafloor_u',round(avg(@BTbtmu),0.01),
-#				   'BT_avg_seafloor_v',round(avg(@BTbtmv),0.01),
-#				   'BT_avg_seafloor_w',round(avg(@BTbtmw),0.01),
-#				   'BT_rms_tilt',round(rms(@BTtilt),0.1));
 }
 
 1;
--- a/edit_data.pl	Thu Nov 21 10:24:23 2013 +0000
+++ b/edit_data.pl	Tue May 20 09:08:43 2014 +0000
@@ -1,9 +1,9 @@
 #======================================================================
 #                    E D I T _ D A T A . P L 
 #                    doc: Sat May 22 21:35:55 2010
-#                    dlm: Tue Nov 12 03:09:49 2013
+#                    dlm: Mon May 19 22:24:40 2014
 #                    (c) 2010 A.M. Thurnherr
-#                    uE-Info: 33 29 NIL 0 0 72 2 2 4 NIL ofnI
+#                    uE-Info: 285 28 NIL 0 0 72 2 2 4 NIL ofnI
 #======================================================================
 
 # HISTORY:
@@ -27,6 +27,7 @@
 #				  - added correctAttitude()
 #	Oct 15, 2012: - BUG: editSurfLayer() counted also ensembles without CTD depth
 #	Nov 12, 2013: - added comments on editCorr_Earthcoords()
+#	Mar  4, 2013: - added support for missing PITCH/ROLL (TILT) & HEADING
 
 # NOTES:
 #	- editCorr_Earthcoords() is overly conservative and removed most
@@ -44,9 +45,9 @@
 sub correctAttitude($$$$)
 {
 	my($ens,$pitch_bias,$roll_bias,$heading_bias) = @_;
-	$LADCP{ENSEMBLE}[$ens]->{PITCH}   -= $pitch_bias;
-	$LADCP{ENSEMBLE}[$ens]->{ROLL}    -= $roll_bias;
-	$LADCP{ENSEMBLE}[$ens]->{HEADING} -= $heading_bias;
+	$LADCP{ENSEMBLE}[$ens]->{PITCH}   -= $pitch_bias 	if defined($LADCP{ENSEMBLE}[$ens]->{PITCH});
+	$LADCP{ENSEMBLE}[$ens]->{ROLL}    -= $roll_bias		if defined($LADCP{ENSEMBLE}[$ens]->{ROLL});
+	$LADCP{ENSEMBLE}[$ens]->{HEADING} -= $heading_bias	if defined($LADCP{ENSEMBLE}[$ens]->{HEADING});
 }
 
 #======================================================================
@@ -138,18 +139,13 @@
 	$LADCP{ENSEMBLE}[$ens]->{TILT} =
 		&angle_from_vertical($LADCP{ENSEMBLE}[$ens]->{PITCH},$LADCP{ENSEMBLE}[$ens]->{ROLL});
 
-	return 0 if ($LADCP{ENSEMBLE}[$ens]->{TILT} <= $lim);
+	return 0 unless ($LADCP{ENSEMBLE}[$ens]->{TILT} > $lim);
 
 	my($nrm) = 0;
 	for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
 		next unless defined($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
 		undef($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
 		$nrm++;
-#		for (my($beam)=0; $beam<4; $beam++) {
-#			next unless defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$beam]);
-#			undef($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$beam]);
-#			$nrm++;
-#		}
 	}
 	return $nrm;
 }
@@ -282,6 +278,43 @@
 
 
 #======================================================================
+# ($nvrm,$nerm) = editPPI($fromEns,$toEns,$range)
+#
+# NOTES:
+#	1) When this code is executed the travel-time profile (@ttProf at 1m resolution)
+#	   has been constructed.
+#======================================================================
+
+sub editPPI($$$)
+{
+	my($fe,$te,$wd) = @_;	# first & last ens to process, water depth for downlooker
+	my($nvrm) = 0;			# of velocities removed
+	my($nerm) = 0;			# of ensembles affected
+	for (my($e)=$fe; $e<=$te; $e++) {
+		next unless numberp($LADCP{ENSEMBLE}[$e]->{CTD_DEPTH});
+		my($range) = $LADCP{ENSEMBLE}[$e]->{XDUCER_FACING_UP}
+				   ? $LADCP{ENSEMBLE}[$e]->{CTD_DEPTH}
+				   : $wd - $LADCP{ENSEMBLE}[$e]->{CTD_DEPTH};
+		my($sscorr) = $CTD{SVEL}[$LADCP{ENSEMBLE}[$e]->{CTD_SCAN}] / 1500;
+		my($goodBins) =   ($range - $sscorr*$LADCP{DISTANCE_TO_BIN1_CENTER}) * cos(rad($LADCP{BEAM_ANGLE}))
+						/ ($sscorr*$LADCP{BIN_LENGTH})
+						- 1.5;
+
+		my($dirty) = 0;
+		for (my($bin)=int($goodBins); $bin<$LADCP{N_BINS}; $bin++) { 	# NB: 2 good bins implies that bin 2 is bad
+			next unless ($bin>=0 && defined($LADCP{ENSEMBLE}[$e]->{W}[$bin]));
+			$dirty = 1;
+			$nvrm++;
+			undef($LADCP{ENSEMBLE}[$e]->{W}[$bin]);
+		}
+
+		$nerm += $dirty;
+	}
+	return ($nvrm,$nerm);
+}
+
+
+#======================================================================
 # $nerm = editSurfLayer($fromEns,$toEns,$surface_layer_depth)
 #
 # NOTES:
--- a/svel_corrections.pl	Thu Nov 21 10:24:23 2013 +0000
+++ b/svel_corrections.pl	Tue May 20 09:08:43 2014 +0000
@@ -1,9 +1,9 @@
 #======================================================================
 #                    S V E L _ C O R R E C T I O N S . P L 
 #                    doc: Thu Dec 30 01:35:18 2010
-#                    dlm: Wed Oct 26 21:54:12 2011
+#                    dlm: Thu Apr 17 09:02:29 2014
 #                    (c) 2010 A.M. Thurnherr
-#                    uE-Info: 15 69 NIL 0 0 72 2 2 4 NIL ofnI
+#                    uE-Info: 49 24 NIL 0 0 72 0 2 4 NIL ofnI
 #======================================================================
 
 # HISTORY:
@@ -13,6 +13,7 @@
 #	Oct 26, 2011: - BUG? in calc_binDepth() on very shallow station 38 in
 #					2010 Gom Spill data set the uplooker code did not stop
 #				 	at the surface, requiring additon of another test
+#	Mar  4, 2014: - added support for missing TILT (PITCH/ROLL)
 
 # NOTES:
 #	In an effort to track down the scale bias, NBP0901 stn 160 was reprocessed with various
@@ -43,6 +44,10 @@
 	my($ens) = @_;
 	my(@bindz);
 
+	# if the following assertion fails, the entire code needs to be searched for
+	# each call of calc_binDepths() needs to be protected by a test
+	die("ensemble $ens") unless defined($LADCP{ENSEMBLE}[$ens]->{TILT});
+
 	my($tanSqBeamAngle) = tan(rad($LADCP{BEAM_ANGLE}))**2;
 	my($curdz) = 0;												# calc avg sndspeed btw transducer & 1st bin
 	$curdz-- until numberp($sVelProf[int($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH}+$curdz)]);
--- a/time_lag.pl	Thu Nov 21 10:24:23 2013 +0000
+++ b/time_lag.pl	Tue May 20 09:08:43 2014 +0000
@@ -1,9 +1,9 @@
 #======================================================================
 #                    T I M E _ L A G . P L 
 #                    doc: Fri Dec 17 21:59:07 2010
-#                    dlm: Sat May 18 11:35:50 2013
+#                    dlm: Thu Apr 17 08:54:30 2014
 #                    (c) 2010 A.M. Thurnherr
-#                    uE-Info: 60 15 NIL 0 0 72 0 2 4 NIL ofnI
+#                    uE-Info: 74 0 NIL 0 0 72 0 2 4 NIL ofnI
 #======================================================================
 
 # HISTORY:
@@ -50,6 +50,7 @@
 #	Mar 23, 2012: - adapted to piece-wise time lagging
 #	Apr 22, 2013: - replaced $max_allowed_w by $opt_m, $TL_required_top_three_fraction by $opt_3
 #	May 14, 2013: - opt_m => w_max_lim
+#	Mar  3, 2014: - BUG: var-name typo
 
 # DIFFICULT STATIONS:
 #	NBP0901#131		this requires the search-radius doubling heuristic
@@ -226,7 +227,7 @@
 		} else {
 			warning(0,"lag too close to edge of search --- trying again after doubling the search radius\n");
 			$search_radius *= 2;
-			$search_raidus =- $w_size if ($search_radius > $w_size);
+			$search_radius =- $w_size if ($search_radius > $w_size);
 		}
 		undef(%nBest); undef(%madBest); undef(@best_lag);
 		goto RETRY;
@@ -238,7 +239,7 @@
 		} else {
 			warning(0,"lag too close to edge of search --- trying again after doubling the search radius\n");
 			$search_radius *= 2;
-			$search_raidus =- $w_size if ($search_radius > $w_size);
+			$search_radius =- $w_size if ($search_radius > $w_size);
 		}
 		undef(%nBest); undef(%madBest); undef(@best_lag);
 		goto RETRY;