LADCP_w: comparison LADCP_w

equal deleted inserted replaced

-:d9309804b6cf
+:5006e9158207
 #!/usr/bin/perl
 #======================================================================
 #                    L A D C P _ W _ P O S T P R O C
 #                    doc: Fri Apr 24 17:15:59 2015
-#                    dlm: Thu May 26 18:21:11 2016
+#                    dlm: Thu Nov  1 10:55:34 2018
 #                    (c) 2015 A.M. Thurnherr
-#                    uE-Info: 101 35 NIL 0 0 72 2 2 4 NIL ofnI
+#                    uE-Info: 71 77 NIL 0 0 72 2 2 4 NIL ofnI
 #======================================================================
 $antsSummary = 'edit and re-grid LADCP vertical-velocity samples';
 # HISTORY:
 #	Mar 16, 2016: - adapted to gmt5
 #   Mar 31, 2016: - changed version %PARAM
 #	Apr 14, 2016: - added profile id to warning messages
 #	May 24, 2016: - improved plot
 #	May 26, 2016: - added automatic directory creation on -d
+#	Nov 28, 2017: - removed wcorr plot
+#	Dec  9, 2017: - added $antsSuppressCommonOptions = 1;
+#	Oct 31, 2018: - improved label (no longer explained/residual stddev)
+#	Nov  1, 2018: - made layout consistent for dual- and single-head profiles
 ($ANTS)  = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
 ($WCALC) = ($0              =~ m{^(.*)/[^/]*$});
 $WCALC = '.' if ($WCALC eq '');
 require "$ANTS/ants.pl";
 require "$ANTS/libstats.pl";
 require "$ANTS/libGMT.pl";
 &antsAddParams('LADCP_w_postproc::version',$VERSION);
+$antsSuppressCommonOptions = 1;
 &antsUsage('b:d:i:k:l:o:p:v:w:',1,
 	'[profile -i)d <id>]',
 	'[-o)utput bin <resolution>] [-k) require <min> samples]',
 	'[-v)alid bins <DL first>,<DL last>[,<UL first>,<UL_last>]',
 	'[-w) <DL_dc_field>,<DL_uc_field>[,<UL_dc_field>,<UL_uc_field>]',
 #----------------------------------------------------------------------
 # Redirect STDOUT to %.wprof & create plots if STDOUT is a tty
 #----------------------------------------------------------------------
 if (-t STDOUT) {
-	$opt_p = defined($opt_d) ? "$opt_d/%03d_wprof.ps,$opt_d/%03d_wcorr.ps"
+#	$opt_p = defined($opt_d) ? "$opt_d/%03d_wprof.ps,$opt_d/%03d_wcorr.ps"
-							 : '%03d_wprof.ps,%03d_wcorr.ps'
+#							 : '%03d_wprof.ps,%03d_wcorr.ps'
+#		unless defined($opt_p);
+	$opt_p = defined($opt_d) ? "$opt_d/%03d_wprof.ps" : '%03d_wprof.ps'
 		unless defined($opt_p);
 	$outfile = defined($opt_d) ? sprintf('%s/%03d.wprof',$opt_d,$id)
 							   : sprintf('%03d.wprof',$id);
 	open(STDOUT,">$outfile") || die("$outfile: $!\n");
 }
 		$sxx += $xt * $xt;
 		$syy += $yt * $yt;
 		$sxy += $xt * $yt;
 	}
 	my($R) = $sxy/(sqrt($sxx * $syy) + 1e-16);
-	my($esig) = sqrt(($R**2)*$sxx/$n);
+	my($rms) = sqrt($sxx/$n);
-	my($rsig) = sqrt((1-$R**2)*$sxx/$n);
+	return ($R,$rms);
-	return ($R,$esig,$rsig);
 }
 sub uc_corr(@)
 {
 	my($n) = 0;
 		$sxx += $xt * $xt;
 		$syy += $yt * $yt;
 		$sxy += $xt * $yt;
 	}
 	my($R) = $sxy/(sqrt($sxx * $syy) + 1e-16);
-	my($esig) = sqrt(($R**2)*$sxx/$n);
+	my($sig) = sqrt(($sxx+$syy)/(2*$n));			# stddev of average w signal from DL & UL
-	my($rsig) = sqrt((1-$R**2)*$sxx/$n);
+	return ($R,$sig);
-	return ($R,$esig,$rsig);
 }
 #----------------------------------------------------------------------
 $dcwF = &fnr($Ddwf); $ucwF = &fnr($Duwf);
 	unless defined($dayNoP);
 if (defined($opt_p)) {
 	($sumPF,$corrPF) = split(/,/,$opt_p);
 	croak("$0: cannot decode -p $opt_p\n")
-		unless (length($corrPF)>0);
+		unless (length($sumPF)>0);
 }
 my($R,$R2);
 if (defined($opt_p)) {												# begin summary plot
 	$xmin = -0.1; $x2min = -700;
 					  ? $DL_dc_median[$bi] - $UL_dc_median[$bi] : nan;
 		$uc_diff[$bi] = numberp($DL_uc_median[$bi]) && numberp($UL_uc_median[$bi])
 					  ? $DL_uc_median[$bi] - $UL_uc_median[$bi] : nan;
 	}
-	($dc_R,$dc_esig,$dc_rsig) = &dc_corr();										# correlation statistics
+	($dc_R,$dc_sig) = &dc_corr();												# correlation statistics
-	($uc_R,$uc_esig,$uc_rsig) = &uc_corr();
+	($uc_R,$uc_sig) = &uc_corr();
-	&antsAddParams('dc_R',$dc_R,'uc_R',$uc_R,
+	&antsAddParams('dc_R',$dc_R,'uc_R',$uc_R,'dc_w.sig',$dc_sig,'uc_w.sig',$uc_sig);
-				   'dc_explained_stddev',$dc_esig,'uc_explained_stddev',$uc_esig,
-				   'dc_residual_stddev',$dc_rsig,'uc_residual_stddev',$uc_rsig);
 	if (defined($opt_p)) {														# plot 2nd-instrument profiles
 		GMT_psxy('-W1,coral,.');
 		for (my($bi)=0; $bi<=$#dcw1; $bi++) {
 			printf(GMT "%f %f\n",$UL_dc_median[$bi],($bi+0.5)*$opt_o);
 #----------------------------------------------------------------------
 # Average and Output Profiles, Add to Summary Plot
 #----------------------------------------------------------------------
+@antsNewLayout = ('depth','hab',
+				  'dc_elapsed','dc_w','dc_w.mad','dc_w.nsamp',
+'uc_elapsed','uc_w','uc_w.mad','uc_w.nsamp',
+				  'dc_w.diff','uc_w.diff');										# DL-UL differences
 if ($dual_head) {																# dual-head output
-	@antsNewLayout = ('depth','hab',
-					  'dc_elapsed','dc_w','dc_w.mad','dc_w.nsamp',
-	                  'uc_elapsed','uc_w','uc_w.mad','uc_w.nsamp',
-					  'dc_w.diff','uc_w.diff');									# DL-UL differences
 	&antsAddParams('profile_id',$id,'lat',$P{lat},'lon',$P{lon});				# selected %PARAMs
 	&antsAddParams($dayNoP,$dn,'run_label',"$first_label & $P{run_label}");
 	&antsAddParams('outgrid_dz',$opt_o,'outgrid_minsamp',$opt_k);
 	&antsAddParams('min_depth',round($min_depth),'max_depth',round($max_depth));
 	&antsAddParams('water_depth',$P{water_depth},'water_depth.sig',$P{water_depth.sig});
 	&antsAddParams('ADCP_bin_length',$bin_length,'ADCP_pulse_length',$pulse_length);
 	&antsAddParams('ADCP_blanking_distance',$blanking_dist);
 	undef($antsOldHeaders);
-} else {
-	@antsNewLayout = ('depth','hab',											# single-head output
-					  'dc_elapsed','dc_w','dc_w.mad','dc_w.nsamp',
-	                  'uc_elapsed','uc_w','uc_w.mad','uc_w.nsamp');
 }
 #&antsInfo("WARNING: unknown water depth (no height-above-bottom)")
 #	unless numberp($P{water_depth});
 			 (($dcns[$bi]>=$opt_k)?$dcwm[$bi]:nan),(($dcns[$bi]>=$opt_k)?$dcwmad[$bi]:nan),
 			 scalar(@{$dcw[$bi]}),
 			 avg(@{$uce[$bi]}),
 			 (($ucns[$bi]>=$opt_k)?$ucwm[$bi]:nan),(($ucns[$bi]>=$opt_k)?$ucwmad[$bi]:nan),
 			 scalar(@{$ucw[$bi]}));
-	push(@{$out[$bi]},$dc_diff[$bi],$uc_diff[$bi])
+	if ($dual_head) {
-		if ($dual_head);
+		push(@{$out[$bi]},$dc_diff[$bi],$uc_diff[$bi]);
+	} else {
+		push(@{$out[$bi]},nan,nan);
+	}
 	&antsOut(@{$out[$bi]});
 }
 if (defined($opt_p)) {																# complete summary plot
 	GMT_setR($R);
 		GMT_psxy('-W1,100/100/255');												# surface layer limit
 			print(GMT "-0.1 $opt_s\n0.07 $opt_s\n");
 		GMT_unitcoords();
 		if ($dc_R < 0.3 || !numberp($dc_R)) {										# correlation statistics
 			&antsInfo("WARNING: low dc correlation (r = %.1f) between UL and DL data in profile #$id",$dc_R);
-			GMT_pstext('-F+f12,Helvetica,coral+jTL -Gred');
+									GMT_pstext('-F+f12,Helvetica,coral+jTL -Gred');
 		} elsif ($dc_R < 0.5) { 	GMT_pstext('-F+f12,Helvetica,coral+jTL -Gyellow'); }
 		else {						GMT_pstext('-F+f12,Helvetica,coral+jTL -Gwhite'); }
-	        printf(GMT "%f %f r = %.1f\n",0.61,0.01,$dc_R);
+	        printf(GMT "%f %f r = %.1f\n",0.62,0.01,$dc_R);
-	}
+		GMT_pstext('-F+f12,Helvetica,coral+jTR -Gwhite');
-	GMT_pstext('-F+f12,Helvetica,coral+jTR -Gwhite');
+			printf(GMT "%f %f [%.1f cm/s stddev]\n",0.99,0.01,100*$dc_sig);
-		printf(GMT "%f %f [%.1f/%.1f cm/s @~s@~\@-e/r\@-]\n",
-			0.99,0.01,100*$dc_esig,100*$dc_rsig) if ($dual_head);
-	if ($dual_head) {
 		if ($uc_R < 0.3 || !numberp($uc_R)) {
 			&antsInfo("WARNING: low uc correlation (r = %.1f) between UL and DL data in profile #$id",$uc_R);
-			GMT_pstext('-F+f12,Helvetica,SeaGreen+jTL -Gred');
+									GMT_pstext('-F+f12,Helvetica,SeaGreen+jTL -Gred');
 		} elsif ($uc_R < 0.5) { 	GMT_pstext('-F+f12,Helvetica,SeaGreen+jTL -Gyellow'); }
 		else {						GMT_pstext('-F+f12,Helvetica,SeaGreen+jTL -Gwhite'); }
-	        printf(GMT "%f %f r = %.1f\n",0.61,0.05,$uc_R);
+	        printf(GMT "%f %f r = %.1f\n",0.62,0.05,$uc_R);
-	}
+		GMT_pstext('-F+f12,Helvetica,SeaGreen+jTR -Gwhite');
-	GMT_pstext('-F+f12,Helvetica,SeaGreen+jTR -Gwhite');
+			printf(GMT "%f %f [%.1f cm/s stddev]\n",0.99,0.05,100*$uc_sig);
-		printf(GMT "%f %f [%.1f/%.1f cm/s @~s@~\@-e/r\@-]\n",
+	}
-			0.99,0.05,100*$uc_esig,100*$uc_rsig) if ($dual_head);
 	GMT_pstext('-F+f14,Helvetica,blue+jTL -N');
 	if (defined($outfile)) { print(GMT "0.01 -0.06 $outfile [$P{run_label}]\n"); }
 	else { 					printf(GMT "0.01 -0.06 %03d\n [$P{run_label}]",$id); }
 	GMT_pstext('-F+f12,Helvetica+jMR');
 	GMT_pstext('-F+f9,Helvetica,orange+jBR -N -Gwhite');
 		print(GMT "0.99 0.99 V$VERSION\n");
 	GMT_end();
-	if ($dual_head) {																# correlation plot
+	if ($dual_head && length($corrPF)>0) {												# correlation plot
 		my($mwm) = 0.05; # max(|w|) for axes
 		for (my($bi)=0; $bi<@DL_dc_median; $bi++) {
 			next unless numberp($DL_dc_median[$bi]) && numberp($UL_dc_median[$bi]);
 			$mwm = abs($DL_dc_median[$bi]) if abs($DL_dc_median[$bi]) > $mwm;
 			$mwm = abs($UL_dc_median[$bi]) if abs($UL_dc_median[$bi]) > $mwm;
 			if (defined($outfile)) { printf(GMT "%f %f $outfile [$P{run_label}]\n",-$mwm,1.1*$mwm); }
 		    else { 					 printf(GMT "%f %f %03d\n [$P{run_label}]",$id,-$mwm,1.1*$wmw); }
 		GMT_psbasemap('-Bf0.01a0.05:"DL Vertical Velocity [m/s]":/f0.01a0.05:"UL Vertical Velocity [m/s]":WeSn');
 		GMT_end();
-} # if dual_head
+} # if dual_head && length(corrPF) > 0
 }
 &antsExit(0);

changeset 49	5006e9158207
parent 47	2ccb81b7cea5
child 56	8f120b9f795a