before beginning to use instrument tilt in distance calculations
authorA.M. Thurnherr <athurnherr@yahoo.com>
Fri, 21 Mar 2014 15:13:50 +0000
changeset 26 d778b73f2a43
parent 25 91bd907db97f
child 27 ec8873454890
before beginning to use instrument tilt in distance calculations
LADCP2CTDmatch
LADCPintsh
LADCPproc
LADCPproc.defaults
deleted file mode 100755
--- a/LADCP2CTDmatch
+++ /dev/null
@@ -1,333 +0,0 @@
-#!/usr/bin/perl
-#======================================================================
-#                    L A D C P 2 C T D M A T C H 
-#                    doc: Thu Sep 23 22:01:19 2010
-#                    dlm: Tue Sep 28 21:43:29 2010
-#                    (c) 2010 A.M. Thurnherr
-#                    uE-Info: 312 33 NIL 0 0 70 10 2 4 NIL ofnI
-#======================================================================
-
-$antsSummary = 'match LADCP to CTD time series';
-
-# HISTORY:
-#	Sep 23, 2010: - incepted
-
-($ANTS) = ($0 =~ m{^(.*)/[^/]*$});
-$PERL_TOOLS = '/Data/LADCP/Software/Thurnherr/perl-tools';
-
-require "$ANTS/ants.pl";
-require "$PERL_TOOLS/RDI_BB_Read.pl";
-require "$PERL_TOOLS/RDI_Coords.pl";
-require "$PERL_TOOLS/RDI_Utils.pl";
-
-$antsParseHeader = 0;
-&antsUsage('dw:n:',2,
-	'[-d)ebug] ',
-	'[-w)indow <size[30s]>] [-n) <windows[20]]',
-	'<RDI file> <SeaBird file>');
-
-&antsCardOpt(\$opt_w,30);
-&antsCardOpt(\$opt_n,20);
-
-$LADCP_file = &antsFileArg();
-$CTD_file 	= &antsFileArg();
-
-&antsAddParams('LADCP_file',$LADCP_file,'CTD_file',$CTD_file);
-@antsNewLayout = ('elapsed','CTD_w','LADCP_w');
-&antsActivateOut();
-
-#----------------------------------------------------------------------
-# Step 1: Read Data
-#----------------------------------------------------------------------
-
-print(STDERR "Reading LADCP data ($LADCP_file)...");
-readData($LADCP_file,\%LADCP);
-printf(STDERR "\n\t%d ensembles\n",scalar(@{$LADCP{ENSEMBLE}}));
-
-print(STDERR "Reading CTD data ($CTD_file)...");
-open(F,$CTD_file) || croak("$CTD_file: $!\n");
-while (1) {											# parse header
-	chomp($hdr = <F>);
-	$hdr =~ s/\r*$//;
-	croak("$0: unexpected EOF (format error)\n") unless defined($hdr);
-	last if ($hdr eq '*END*');
-	
-	$CTD_nfields = $',next if ($hdr =~ /nquan = /);		# Layout
-	$CTD_nrecs = $',next if ($hdr =~ /nvalues = /);
-	$pressF = $1,next if ($hdr =~ /name (\d+) = prDM:/);
-
-	&antsAddParams('start_time',$1),next			# selected metadata
-		if ($hdr =~ /start_time = (.*)/);
-
-	&antsAddParams('station',$1),next
-		if ($hdr =~ /Station\s*:\s*(.*)/);
-	&antsAddParams('ship',$1),next
-		if ($hdr =~ /Ship\s*:\s*(.*)\s*$/);
-	&antsAddParams('cruise',$1),next
-		if ($hdr =~ /Cruise\s*:\s*(.*)\s*$/);
-	&antsAddParams('time',$1),next
-		if ($hdr =~ /Time\s*:\s*(.*)/);
-	&antsAddParams('date',$1),next
-		if ($hdr =~ /Date\s*:\s*(.*)/);
-
-	if ($hdr =~ /Latitude\s*[=:]\s*/) {
-		($deg,$min,$NS) = split(/ /,$');
-		$lat = $deg + $min/60;
-		$lat *= -1 if ($NS eq 'S');
-		&antsAddParams('lat',$lat);
-		next;
-	}
-	if ($hdr =~ /Longitude\s*[=:]\s*/) {
-		($deg,$min,$EW) = split(/ /,$');
-		$lon = $deg + $min/60;
-		$lon *= -1 if ($EW eq 'W');
-		&antsAddParams('lon',$lon);
-		next;
-	}
-	
-	if ($hdr =~ /interval = seconds: /) {
-		$CTD_sampint = 1*$';
-		&antsAddParams('CTD_interval',1/$CTD_sampint);
-		next;
-	}
-	
-	$CTD_badval = $',next
-		if ($hdr =~ /bad_flag = /);	
-	$CTD_file_type = $',next
-		if ($hdr =~ /file_type = /);	
-}
-
-croak("$CTD_file: cannot determine CTD file layout\n")
-	unless ($CTD_nfields && $CTD_nrecs);
-croak("$CTD_file: cannot determine missing value\n")
-	unless defined($CTD_badval);
-croak("$CTD_file: not a CTD time series file\n")
-	unless ($CTD_sampint);
-croak("$CTD_file: no pressure field\n")
-	unless defined($pressF);
-
-croak("$0: unknown latitude\n") unless defined($lat);
-
-if ($CTD_file_type eq 'ascii') {
-	while (1) {
-		last unless (@rec = &antsFileIn(F));
-		push(@CTD_press,($rec[$pressF] == $CTD_badval) ? nan : $rec[$pressF]);
-	}
-} elsif ($CTD_file_type eq 'binary') {
-
-	my($fbits) = 8 * length(pack('f',0));
-	croak(sprintf("$0: incompatible native CPU float representation (%d instead of 32bits)\n",fbits))
-		unless ($fbits == 32);	
-		
-	croak("$CTD_file: can't read binary data\n")
-		unless (read(F,$dta,4*$CTD_nfields*$CTD_nrecs) == 4*$CTD_nfields*$CTD_nrecs);
-	print(STDERR "$CTD_file: WARNING: extraneous data at EOF\n") unless eof(F);
-
-	$dta = pack('V*',unpack('N*',$dta))				# big-endian CPU
-		if (unpack('h*', pack('s', 1)) =~ /01/);	# c.f. perlport(1)
-	
-	@dta = unpack("f*",$dta);
-
-	for ($r=0; $r<$CTD_nrecs; $r++) {
-		push(@CTD_press,($dta[$r*$CTD_nfields+$pressF] == $CTD_badval) ? nan : $dta[$r*$CTD_nfields+$pressF]);
-	}
-} else {
-	croak("$CTD_file: unknown CTD file type $CTD_file_type\n");
-}
-
-printf(STDERR "\n\t%d scans\n",scalar(@CTD_press));
-
-#----------------------------------------------------------------------
-# Step 2a: Pre-Process CTD Data
-#----------------------------------------------------------------------
-
-printf(STDERR "Pre-processing data...");
-printf(STDERR "\n\tCTD...");
-
-#------------------------------------
-# calculate w and find deepest record
-#------------------------------------
-$CTD_maxpress = -9e99;
-for (my($r)=1; $r<$CTD_nrecs-1; $r++) {
-	print(STDERR '.') if ($r%500==0);
-	$CTD_w[$r] = 0.99*($CTD_press[$r+1] - $CTD_press[$r-1]) / (2*$CTD_sampint);
-	if ($CTD_press[$r] > $CTD_maxpress) {
-		$CTD_maxpress = $CTD_press[$r];
-		$CTD_bottom = $r;
-    }										
-}
-
-#----------------------------------------------------------------------
-# Step 2b: Pre-Process LADCP Data
-#----------------------------------------------------------------------
-
-print(STDERR "\n\tLADCP...");
-
-#------------------------------------------------------
-# construct a depth-vs-time "profile" from the raw data
-#------------------------------------------------------
-
-my($LADCP_start,$LADCP_end,$LADCP_bottom,$w_gap_time,$zErr,$maxz) =
-    mk_prof(\%LADCP,0,undef,1,6,70,0.1,120);
-croak("\n$LADCP_file: no good ensembles found\n")
-    unless defined($LADCP_start);
-    
-if ($opt_d) {
-	printf(STDERR "\n\t\tStart of cast    : %s (#%5d) at %6.1fm\n",
-	                    $LADCP{ENSEMBLE}[$LADCP_start]->{TIME},
-	                    $LADCP{ENSEMBLE}[$LADCP_start]->{NUMBER},
-	                    $LADCP{ENSEMBLE}[$LADCP_start]->{DEPTH});
-	printf(STDERR "\t\tBottom of cast   : %s (#%5d) at %6.1fm\n",
-	                    $LADCP{ENSEMBLE}[$LADCP_bottom]->{TIME},
-	                    $LADCP{ENSEMBLE}[$LADCP_bottom]->{NUMBER},
-	                    $LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH});
-	if (defined($water_depth)) {
-	    printf(STDERR "\t\tSeabed           :                      at %6.1fm (+-%dm)\n",$water_depth,$sig_wd);
-	} else {
-	    print(STDERR "\t\tSeabed           : not found\n");
-	}
-	printf(STDERR "\t\tEnd of cast      : %s (#%5d) at %6.1fm\n",
-	                    $LADCP{ENSEMBLE}[$LADCP_end]->{TIME},
-	                    $LADCP{ENSEMBLE}[$LADCP_end]->{NUMBER},
-	                    $LADCP{ENSEMBLE}[$LADCP_end]->{DEPTH});
-}
-
-#------------------------------------------------------------------------
-# find 1st rec & ensemble >=80% down to max depth & make 1st guess at lag
-#------------------------------------------------------------------------
-
-my($CTD_80pct_down) = 0;
-$CTD_80pct_down++
-	while ($CTD_press[$CTD_80pct_down] < 0.8*$CTD_maxpress);
-
-my($LADCP_80pct_down) = 0;
-$LADCP_80pct_down++
-	while ($LADCP{ENSEMBLE}[$LADCP_80pct_down]->{DEPTH} < 0.8*$LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH});
-
-
-my($first_guess_lag) = $LADCP{ENSEMBLE}[$LADCP_80pct_down]->{ELAPSED_TIME} -
-					   $CTD_80pct_down*$CTD_sampint;
-
-printf(STDERR "\t\t1st guess offset = %ds\n",$first_guess_lag*$CTD_sampint)
-	if ($opt_d);
-
-#------------------------------------------------------------------------------------
-# Linearly interpolate LADCP time series onto a new grid with $CTD_sampint resolution
-#	ALSO: apply first_guess_lag to make lags small, which keeps the bestlag data
-#		  chunks large
-#------------------------------------------------------------------------------------
-
-my($nGaps) = 0;
-
-sub interp_LADCP_w($$)
-{
-	my($elapsed,$ens) = @_;
-	my($sc) = ($elapsed - $LADCP{ENSEMBLE}[$ens-1]->{ELAPSED_TIME}) /
-			  ($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} -
-					$LADCP{ENSEMBLE}[$ens-1]->{ELAPSED_TIME});
-	unless (numberp($LADCP{ENSEMBLE}[$ens-1]->{W})) {
-		$nGaps++;
-		croak("$0: cannot (yet) handle double gaps --- request feature upgrade\n")
-			unless numberp($LADCP{ENSEMBLE}[$ens]->{W});
-		return $LADCP{ENSEMBLE}[$ens]->{W};
-	}
-	return $LADCP{ENSEMBLE}[$ens-1]->{W} +
-				$sc * ($LADCP{ENSEMBLE}[$ens]->{W} - $LADCP{ENSEMBLE}[$ens-1]->{W});
-}
-
-for (my($ens)=$LADCP_start,my($r)=0; $ens<=$LADCP_end; $ens++) {
-	while ($r*$CTD_sampint < $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}) {
-		$LADCP_w[$r-$first_guess_lag] = interp_LADCP_w($r*$CTD_sampint,$ens)
-			unless ($first_guess_lag > $r);
-		$r++;
-    }
-}
-
-print(STDERR "\t\t$nGaps gaps in w timeseries")
-	if ($opt_d);
-
-print(STDERR "\n");
-
-#----------------------------------------------------------------------
-# Step 3: Calculate Lags
-#----------------------------------------------------------------------
-
-printf(STDERR "Calculating $opt_n lags from %ds-long windows [s]: ",$opt_w);
-$opt_w = int($opt_w / $CTD_sampint);
-
-sub bestLag($)
-{
-	my($ws) = @_;														# window start index
-
-	my($best);
-	my($bestmad) = 9e99;												# mean absolute deviation
-	for (my($Llag)=-int($opt_w/2); $Llag<int($opt_w/2); $Llag++) {
-		my($sad) = my($nad) = 0;
-		for (my($Ci)=0; $Ci<$opt_w; $Ci++) {
-			my($Li) = $Ci + $Llag;
-			next if ($Li<0 || $Li>=$opt_w);
-			next unless numberp($CTD_w[$ws+$Ci]) && numberp($LADCP_w[$ws+$Li]);
-			$sad += abs($CTD_w[$ws+$Ci] - $LADCP_w[$ws+$Li]);
-			$nad++;
-		}
-		next unless ($nad > 0);
-		if ($sad/$nad < $bestmad) {
-			$best = $Llag;
-			$bestmad = $sad/$nad;
-		}
-	}
-	return $best;	
-}
-
-#---------------------------------------------------------------
-# carry out opt_n lag correlations and keep tally of the results
-#---------------------------------------------------------------
-my(@nBest);
-for (my($window)=0; $window<$opt_n; $window++) {
-	my($ws) = $window * int(@LADCP_w/$opt_n);			# window start
-	$ws = @LADCP_w-$opt_w if ($ws+$opt_w >= @LADCP_w);
-
-	$bestLag = bestLag($ws);
-	if (defined($bestLag)) {
-		printf(STDERR "%d ",$bestLag*$CTD_sampint);
-		$nBest{$bestLag}++;
-	} else {
-		printf(STDERR "nan ");
-	}
-}
-
-#----------------------
-# find most popular lag
-#----------------------
-my($best_lag);
-foreach my $i (keys(%nBest)) {
-	$best_lag = $i if ($nBest{$i} > $nBest{$best_lag});
-}
-croak("\n$0: cannot determine a valid lag\n")
-	unless ($nBest{$best_lag} > 1);
-print(STDERR "\n\n\tWARNING: only $nBest{$best_lag} of the lag estimates agree!\n")
-	if ($nBest{$best_lag} < $opt_n/2);
-
-if ($nBest{$best_lag} == $opt_n) {
-	printf(STDERR "\n\tunanimous lag = %ds\n",$best_lag*$CTD_sampint);
-} else {
-	printf(STDERR "\n\tmost popular lag = %ds\n",$best_lag*$CTD_sampint);
-}
-
-#----------------------------------------------------------------------
-# Step 4: associate CTD data with each LADCP ensemble & write time series
-#----------------------------------------------------------------------
-
-print(STDERR "Associating CTD data with LADCP ensembles...\n");
-
-my($overall_dt) = $CTD_sampint * ($best_lag + $first_guess_lag);
-for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
-	my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} - $overall_dt) / $CTD_sampint);
-	$LADCP{ENSEMBLE}[$ens]->{CTD_W} = $CTD_w[$r];
-	&antsOut($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME},
-			 $LADCP{ENSEMBLE}[$ens]->{CTD_W},
-			 $LADCP{ENSEMBLE}[$ens]->{W});
-}
-
-&antsExit();
-
--- a/LADCPintsh
+++ b/LADCPintsh
@@ -2,9 +2,9 @@
 #======================================================================
 #                    L A D C P I N T S H 
 #                    doc: Thu Oct 14 21:22:50 2010
-#                    dlm: Thu Mar 20 11:57:54 2014
+#                    dlm: Thu Mar 20 12:14:53 2014
 #                    (c) 2010 A.M. Thurnherr & E. Firing
-#                    uE-Info: 58 55 NIL 0 0 72 2 2 4 NIL ofnI
+#                    uE-Info: 266 0 NIL 0 0 72 2 2 4 NIL ofnI
 #======================================================================
 
 $antsSummary = 'integrate LADCP shear';
@@ -263,7 +263,6 @@
 			$wz[$r] = $uc_wz[$r];
 			$elapsed[$r] = $uc_elapsed[$r];
         }
-		print(STDERR "uz[$r] := $uz[$r] [$dc_uz[$r]/$uc_uz[$r]] (dcf=$dcf ucf=$ucf)\n");
 	} else {
 		$uz[$r] = $vz[$r] = $wz[$r] = $elapsed[$r] = nan;
 	}
--- a/LADCPproc
+++ b/LADCPproc
@@ -2,9 +2,9 @@
 #======================================================================
 #                    L A D C P P R O C 
 #                    doc: Thu Sep 16 20:36:10 2010
-#                    dlm: Wed Mar 19 23:41:55 2014
+#                    dlm: Thu Mar 20 12:34:04 2014
 #                    (c) 2010 A.M. Thurnherr
-#                    uE-Info: 91 55 NIL 0 0 72 10 2 4 NIL ofnI
+#                    uE-Info: 94 52 NIL 0 0 72 10 2 4 NIL ofnI
 #======================================================================
 
 # NOTES:
@@ -89,6 +89,9 @@
 #	Mar 19, 2014: - moved code to set LADCP_time_lag %PARAM into main prog so it is
 #				    set, even when -l is used
 #				  - added pitch, roll, hdg to -t output
+#	Mar 20, 2014: - BUG: wrong number of samples were recorded when upcast had no
+#						 valid data whatsoever
+#				  - added support for $LADCP_max_gap
 
 ($ANTS)    = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
 ($PERL_TOOLS) = (`which mkProfile` =~ m{^(.*)/[^/]*$});
@@ -357,7 +360,7 @@
 	if ($opt_d);
 	
 ($LADCP_start,$LADCP_end,$LADCP_bottom,$w_gap_time,$zErr,$maxz) =
-    mk_prof(\%LADCP,0,undef,1,6,70,0.1,9999);
+    mk_prof(\%LADCP,0,undef,1,6,70,0.1,$LADCP_max_gap);
 croak("\n$LADCP_file: no good ensembles found\n")
     unless defined($LADCP_start);
 
@@ -628,6 +631,8 @@
 @dc_lash_mu = @lash_mu; @dc_losh_mu = @losh_mu;
 
 print(STDERR "\n\tupcast...") if ($opt_d);
+@sh_n=@ush_mu=@ush_sig=@vsh_mu=@vsh_sig=@wsh_mu=@wsh_sig=@esh_mu=@lash_mu=@losh_mu=undef;
+
 edit_velocity($LADCP_end,$LADCP_bottom);							# upcast
 calc_shear($LADCP_end,$LADCP_bottom,$SHEAR_PREGRID_DZ,0);
 calc_shear($LADCP_end,$LADCP_bottom,$GRID_DZ,1);
--- a/LADCPproc.defaults
+++ b/LADCPproc.defaults
@@ -1,9 +1,9 @@
 #======================================================================
 #                    L A D C P P R O C . D E F A U L T S 
 #                    doc: Fri Sep 17 09:44:21 2010
-#                    dlm: Thu Sep 19 13:34:34 2013
+#                    dlm: Thu Mar 20 12:41:41 2014
 #                    (c) 2010 A.M. Thurnherr
-#                    uE-Info: 142 75 NIL 0 0 72 0 2 4 NIL ofnI
+#                    uE-Info: 81 23 NIL 0 0 72 0 2 4 NIL ofnI
 #======================================================================
 
 # default parameters for [LADCPproc]
@@ -32,6 +32,7 @@
 #	Sep  6, 2013: - BUG: BT_begin_search_above value of 300m was correct; 
 #						 the original bug was in the documentation
 #	Sep 19, 2013: - added support for $BT_range_method
+#	Mar 20, 2014: - added support for $LADCP_max_gap
 
 #----------------------------------------------------------------------
 # Data editing
@@ -45,6 +46,7 @@
 
 # $bad_beam = 1;
 
+
 #----------------------------------------------------------------------
 # ASCII CTD file support
 #----------------------------------------------------------------------
@@ -69,6 +71,18 @@
 # $CTD_ASCII_badval 		= 999;
 
 #----------------------------------------------------------------------
+# Time Lagging
+#----------------------------------------------------------------------
+
+# Calculation of the LADCP time series of depth is re-started
+# after a gap exceeding $LADCP_max_gap seconds. This is useful
+# for dealing with double dips, data from instruments that
+# return occasional velocities while on deck, etc. By default,
+# all gaps are ignored.
+
+$LADCP_max_gap = 9999;
+
+#----------------------------------------------------------------------
 # Backscatter Profile Parameters
 #----------------------------------------------------------------------