ADCP_tools: comparison mkProfile

equal deleted inserted replaced

-:9c3b147b4372
+:0f89b1523648
 #!/usr/bin/perl
 #======================================================================
 #                    M K P R O F I L E
 #                    doc: Sun Jan 19 18:55:26 2003
-#                    dlm: Fri Apr 12 09:24:14 2013
+#                    dlm: Tue May 14 11:37:04 2013
 #                    (c) 2003 A.M. Thurnherr
-#                    uE-Info: 269 72 NIL 0 0 72 2 2 4 NIL ofnI
+#                    uE-Info: 169 15 NIL 0 0 72 2 2 4 NIL ofnI
 #======================================================================
 # Make an LADCP Profile by Integrating W (similar to Firing's scan*).
 # HISTORY:
 #				  - BUG: post-recovery rotations were always zero
 #	Sep  9, 2011: - BUG: range calculation for Earth coordinate data included bins without
 #						 valid velocities
 #	Sep 21, 2010: - added %rms_heave_acceleration
 #	Apr 12, 2013: - added -p
+#	May 10, 2013: - BUG: mkProfile bombed when ADCP file is truncated at deepest location
+#	May 14, 2013: - added heading to output
+#				  - added err_vel to output
+#				  - finally removed -d/-g
 # NOTES:
 #	- the battery values are based on transmission voltages (different
 #	  from battery voltages) and reported without units (raw 8-bit a2d
 #	  values)
 $USAGE = "$0 @ARGV";
 die("Usage: $0 " .
 	"[-Q)uiet] [-F)ilter <script>] " .
 	"[-s)uppress checkensemble()] " .
-	"[require -4)-beam solutions] " .
+	"[require -4)-beam solutions] [-d)iscard <beam#>] [apply beamvel-m)ask <file>] " .
 	"[-r)ef-layer <bin|1,bin|6>] [-n) vels <min|2>] " .
 	"[-e)rr-vel <max[0.1]] [-c)orrelation <min>] [-p)ct-good <min[100]>] " .
-	"[-m)ax <gap>] " .
+	"[max -g)ap <len>] " .
-	"[-d)rift <dx,dy>] [-g)ps <start lat,lon/end lat,lon>] " .
 	"[output -f)ields <field[,...]> " .
 	"[-M)agnetic <declination>] [profile -B)ottom <depth>] " .
 	"<RDI file>\n")
-		unless (&Getopts("4AB:F:M:Qd:r:n:e:c:g:f:m:sp:") && @ARGV == 1);
+		unless (&Getopts("4AB:F:M:Qd:g:r:n:e:c:f:m:sp:") && @ARGV == 1);
 $RDI_Coords::minValidVels = 4 if ($opt_4);			# no 3-beam solutions
 require $opt_F if defined($opt_F);					# load filter
 $opt_r = "1,6" 	unless defined($opt_r);				# defaults
 $opt_n = 2	   	unless defined($opt_n);
 $opt_e = 0.1   	unless defined($opt_e);
 $opt_c = 70	   	unless defined($opt_c);
-$opt_m = 120	unless defined($opt_m);
+$opt_g = 120	unless defined($opt_g);
 $opt_p = 100	unless defined($opt_p);
 ($minb,$maxb) = split(',',$opt_r);					# reference layer
 die("$0: can't decode -r $opt_r\n") unless defined($maxb);
-if ($opt_g) {										# GPS info
-	($s_lat,$s_lon,$e_lat,$e_lon) = gps_to_deg($opt_g);
-	$lat = $s_lat/2 + $e_lat/2;
-	$lon = $s_lon/2 + $e_lon/2;
-	$ddx = dist($lat,$s_lon,$lat,$e_lon);
-	$ddy = dist($s_lat,$lon,$e_lat,$lon);
-}
-if ($opt_d) {										# ship drift
-	($ddx,$ddy) = split(',',$opt_d);
-	die("$0: can't decode -d $opt_d\n") unless defined($ddy);
-}
 print(STDERR "Reading $ARGV[0]...");				# read data
 readData($ARGV[0],\%dta);
 print(STDERR "done\n");
 die("$ARGV[0]: not enough bins for choice of -r\n")	# enough bins?
 	unless ($dta{N_BINS} >= $maxb);
 if ($dta{BEAM_COORDINATES}) {						# coords used
 	$beamCoords = 1;
 } elsif (!$dta{EARTH_COORDINATES}) {
 	die("$ARGV[0]: only beam and earth coordinates implemented so far\n");
+}
+if (defined($opt_m) && -r $opt_m) {
+	die("$ARGV[0]: -m only implemented for data collected in beam coordinates\n")
+		unless ($beamCoords);
+	print(STDERR "Masking beam velocities as prescribed in $opt_m...");
+	open(BVM,$opt_m) || die("$opt_m: $!\n");
+	while (<BVM>) {
+		s/#.*//;
+		s/^\s*$//;
+		next if ($_ eq '');
+		my($fe,$te,$db) = split;
+		die("$opt_m: cannot decode $_\n")
+			unless (numberp($fe) && numberp($te) && $te>=$fe && $db>=1 && $db<=4);
+		die("$0: assertion failed")
+			unless ($dta{ENSEMBLE}[$fe-1]->{NUMBER} == $fe &&
+					$dta{ENSEMBLE}[$te-1]->{NUMBER} == $te);
+		for (my($ens)=$fe-1; $ens<=$te-1; $ens++) {
+			$nens++;
+			for (my($bin)=0; $bin<$dta{N_BINS}; $bin++) {
+				undef($dta{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$db-1]);
+			}
+	    }
+}
+	close(BVM);
+	print(STDERR " $nens ensembles edited\n");
+}
+if (defined($opt_d)) {								# discard entire beam
+	die("$ARGV[0]: -d only implemented for data collected in beam coordinates\n")
+		unless ($beamCoords);
+	print(STDERR "Discarding beam-$opt_d velocities...");
+for (my($ens)=0; $ens<=$#{$dta{ENSEMBLE}}; $ens++) {
+for (my($bin)=0; $bin<$dta{N_BINS}; $bin++) {
+undef($dta{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$opt_d-1]);
+}
+}
+	print(STDERR "done\n");
 }
 if (defined($opt_M)) {								# magnetic declination
 	$dta{HEADING_BIAS} = -1*$opt_M;
 } else {
 #	print(STDERR "addFields = $addFields\n");
 #	print(STDERR "\@f = @f\n");
 }
 #======================================================================
-# Misc funs used to decode options
-#======================================================================
-sub dist($$$$)										# distance
-{
-	my($lat1,$lon1,$lat2,$lon2) = @_;
-	my($a) = 6378139; 								# Earth's radius
-	$lat1 = rad($lat1); $lon1 = rad($lon1);
-	$lat2 = rad($lat2); $lon2 = rad($lon2);
-	$ct1 = cos($lat1); $st1 = sin($lat1); $cp1 = cos($lon1);
-	$sp1 = sin($lon1); $ct2 = cos($lat2); $st2 = sin($lat2);
-	$cp2 = cos($lon2); $sp2 = sin($lon2);
-	$cos = $ct1*$cp1*$ct2*$cp2 + $ct1*$sp1*$ct2*$sp2 + $st1*$st2;
-	$cos =  1 if ($cos >  1);
-	$cos = -1 if ($cos < -1);
-	return $a * acos($cos);
-}
-sub deg_to_dec($)									# parse degrees
-{
-	my($deg,$min) = split(':',$_[0]);
-	return $deg + $min/60;
-}
-sub gps_to_deg($)									# decode lat/lon
-{
-	my($start,$end) = split('/',$_[0]);
-	my($sa,$so,$ea,$eo);
-	my($lat,$lon) = split(',',$start);
-	if ($lat =~ m{N$}) 		{ $sa =  deg_to_dec($`); }
-	elsif ($lat =~ m{S$}) 	{ $sa = -deg_to_dec($`); }
-	else 					{ $sa = $lat; }
-	if ($lon =~ m{E$}) 		{ $so =  deg_to_dec($`); }
-	elsif ($lon =~ m{W$}) 	{ $so = -deg_to_dec($`); }
-	else 					{ $so = $lon; }
-	my($lat,$lon) = split(',',$end);
-	if ($lat =~ m{N$}) 		{ $ea =  deg_to_dec($`); }
-	elsif ($lat =~ m{S$}) 	{ $ea = -deg_to_dec($`); }
-	else 					{ $ea = $lat; }
-	if ($lon =~ m{E$}) 		{ $eo =  deg_to_dec($`); }
-	elsif ($lon =~ m{W$}) 	{ $eo = -deg_to_dec($`); }
-	else 					{ $eo = $lon; }
-	return ($sa,$so,$ea,$eo);
-}
-#======================================================================
 # Step 0: Check data & Calculate Ping Rates
 #======================================================================
 unless ($dta{NARROW_BANDWIDTH}) {
 	print(STDERR "WARNING: $0 WIDE BANDWIDTH!\n");
 #======================================================================
 # Step 1: Integrate w & determine water depth
 #======================================================================
 ($firstgood,$lastgood,$atbottom,$w_gap_time,$zErr,$maxz,$rms_heave_accel) =
-	mk_prof(\%dta,!$opt_s,$opt_F,$minb,$maxb,$opt_c,$opt_e,$opt_m,$opt_p);
+	mk_prof(\%dta,!$opt_s,$opt_F,$minb,$maxb,$opt_c,$opt_e,$opt_g,$opt_p);
 unless (($atbottom > $firstgood) && ($lastgood > $atbottom)) {
 	if ($opt_Q) {
 		print(STDERR "$ARGV[0]: no valid cast data found\n");
 		exit(0);
 }
 }
 $gb = ($gb[0]+$gb[1]+$gb[2]+$gb[3]) / 4;
 #======================================================================
-# Step 5: Remove Ship Drift (probably not useful)
+# Step 5: Remove Ship Drift (probably not useful => removed)
 #======================================================================
-if (defined($opt_M) && defined($ddx)) {						# remove barotropic
-	$du = $ddx / $dta{ENSEMBLE}[$lastgood]->{ELAPSED_TIME};# mean drift vel
-	$dv = $ddy / $dta{ENSEMBLE}[$lastgood]->{ELAPSED_TIME};
-	$iu = $dta{ENSEMBLE}[$lastgood]->{X} /				# mean obs vel
-			$dta{ENSEMBLE}[$lastgood]->{ELAPSED_TIME};
-	$iv = $dta{ENSEMBLE}[$lastgood]->{Y} /
-			$dta{ENSEMBLE}[$lastgood]->{ELAPSED_TIME};
-	for ($e=$firstgood; $e<=$lastgood; $e++) {
-		next unless (defined($dta{ENSEMBLE}[$e]->{X}) &&
-					 defined($dta{ENSEMBLE}[$e]->{Y}));
-		$dta{ENSEMBLE}[$e]->{U} -= $du;
-		$dta{ENSEMBLE}[$e]->{V} -= $dv;
-		$dta{ENSEMBLE}[$e]->{X} += $dta{ENSEMBLE}[$e]->{ELAPSED_TIME} * ($du-$iu);
-		$dta{ENSEMBLE}[$e]->{Y} += $dta{ENSEMBLE}[$e]->{ELAPSED_TIME} * ($dv-$iv);
-	}
-}
 #======================================================================
 # Step 6: Pitch, Roll, Rotation
 #======================================================================
 $uprot   /= 360;
 $postrot /= 360;
 $prrms 	 = sqrt($prrms/($lastgood-$firstgood));
 $dnprrms = sqrt($dnprrms/($atbottom-$firstgood));
-$upprrms = sqrt($upprrms/($lastgood-$atbottom));
+if ($lastgood == $atbottom) {
+	print(STDERR "WARNING: $0 NO UPCAST DATA\n");
+	$upprrms = nan;
+} else {
+	$upprrms = sqrt($upprrms/($lastgood-$atbottom));
+}
 $rotrms = sqrt($rotrms/($lastgood-$firstgood));
 #======================================================================
 # PRODUCE OUTPUT
 print("#ANTS# [] $USAGE\n");
 $uFields = "{u} {u_err} {v} {v_err} {x} {x_err} {y} {y_err}"
 	if defined($opt_M);
 print("#ANTS#FIELDS# {ens} {time} {elapsed} {secno} {downcast} " .
-					"{w} {w_err} {depth} {depth_err} {depth_BT} " .
+					"{w} {w_err} {err_vel} {depth} {depth_err} {depth_BT} " .
-					"{pitchroll} {rotation} " .
+					"{pitchroll} {heading} {rotation} " .
 					"$uFields $addFields\n");
 printf("#ANTS#PARAMS# date{$dta{ENSEMBLE}[$firstgood]->{DATE}} " .
 			   "start_time{$dta{ENSEMBLE}[$firstgood]->{TIME}} " .
 			  "bottom_time{$dta{ENSEMBLE}[$atbottom]->{TIME}} " .
 	$uErr, $dta{ENSEMBLE}[$lastgood]->{X}, $x_err,
 	$dta{ENSEMBLE}[$lastgood]->{Y} /
 		$dta{ENSEMBLE}[$lastgood]->{ELAPSED_TIME},
 	$vErr, $dta{ENSEMBLE}[$lastgood]->{Y}, $y_err,
 ) if defined ($opt_M);
-print("#ANTS#PARAMS# start_lat{$s_lat} start_lon{$s_lon} " .
-					  "end_lat{$e_lat} end_lon{$e_lon} " .
-						  "lat{$lat} lon{$lon}\n")
-	if defined($lat);
 if ($dta{TIME_BETWEEN_PINGS} == 0) {
 	 print("#ANTS#PARAMS# pinging_rate{staggered}\n");
 } else {
 	 printf("#ANTS#PARAMS# pinging_rate{%.2f}\n",
 		1/$dta{TIME_BETWEEN_PINGS});
 }
-printf("#ANTS#PARAMS# drift_x{%d} drift_y{%d} " .
-					 "drift_u{%.3f} drift_v{%.3f} " .
-	   "\n",$ddx,$ddy,$du,$dv) if defined($ddx);
 if (defined($water_depth)) {
 	printf("#ANTS#PARAMS# water_depth{%d} sig-water_depth{%d}\n",
 				$water_depth,$sig_wd);
 } else {
 	print("#ANTS#PARAMS# water_depth{nan} sig-water_depth{nan}\n");
 	p($dta{ENSEMBLE}[$e]->{ELAPSED_TIME});
 	p($dta{ENSEMBLE}[$e]->{SECNO});
 	pb($dta{ENSEMBLE}[$e]->{UNIX_TIME} < $dta{ENSEMBLE}[$atbottom]->{UNIX_TIME});
 	p($dta{ENSEMBLE}[$e]->{W});
 	p($dta{ENSEMBLE}[$e]->{W_ERR});
+	p($dta{ENSEMBLE}[$e]->{ERR_VEL});
 	p($dta{ENSEMBLE}[$e]->{DEPTH});
 	p($dta{ENSEMBLE}[$e]->{DEPTH_ERR});
 	p($dta{ENSEMBLE}[$e]->{DEPTH_BT});
 	p($dta{ENSEMBLE}[$e]->{PITCHROLL});
+	p($dta{ENSEMBLE}[$e]->{HEADING});
 	p($dta{ENSEMBLE}[$e]->{ROTATION});
 	if (defined($opt_M)) {
 		p($dta{ENSEMBLE}[$e]->{U}); p($dta{ENSEMBLE}[$e]->{U_ERR});
 		p($dta{ENSEMBLE}[$e]->{V}); p($dta{ENSEMBLE}[$e]->{V_ERR});
 		p($dta{ENSEMBLE}[$e]->{X}); p($dta{ENSEMBLE}[$e]->{X_ERR});

changeset 12	0f89b1523648
parent 11	9c3b147b4372
child 14	8c79b38a7086