1 #!/usr/bin/perl |
1 #!/usr/bin/perl |
2 #====================================================================== |
2 #====================================================================== |
3 # L A D C P _ W _ C T D |
3 # L A D C P _ W _ C T D |
4 # doc: Mon Nov 3 17:34:19 2014 |
4 # doc: Mon Nov 3 17:34:19 2014 |
5 # dlm: Thu Mar 31 05:52:01 2016 |
5 # dlm: Thu May 26 20:42:15 2016 |
6 # (c) 2014 A.M. Thurnherr |
6 # (c) 2014 A.M. Thurnherr |
7 # uE-Info: 62 42 NIL 0 0 72 2 2 4 NIL ofnI |
7 # uE-Info: 639 69 NIL 0 0 72 2 2 4 NIL ofnI |
8 #====================================================================== |
8 #====================================================================== |
9 |
9 |
10 $antsSummary = 'pre-process SBE 9plus CTD data for LADCP_w'; |
10 $antsSummary = 'pre-process SBE 9plus CTD data for LADCP_w'; |
11 |
11 |
12 # HISTORY: |
12 # HISTORY: |
58 # Mar 16, 2016: - adapted to gmt5 |
58 # Mar 16, 2016: - adapted to gmt5 |
59 # Mar 19, 2016: - added support for $VERB environment var |
59 # Mar 19, 2016: - added support for $VERB environment var |
60 # - added $libSBE_quiet flag |
60 # - added $libSBE_quiet flag |
61 # Mar 30, 2016: - added station number of ASCII format |
61 # Mar 30, 2016: - added station number of ASCII format |
62 # Mar 31, 2016: - changed version %PARAM |
62 # Mar 31, 2016: - changed version %PARAM |
|
63 # May 26, 2016: - renamed w[.raw] to CTD_w[.raw] |
|
64 # - added winch_w |
63 |
65 |
64 ($ANTS) = (`which ANTSlib` =~ m{^(.*)/[^/]*$}); |
66 ($ANTS) = (`which ANTSlib` =~ m{^(.*)/[^/]*$}); |
65 ($WCALC) = ($0 =~ m{^(.*)/[^/]*$}); |
67 ($WCALC) = ($0 =~ m{^(.*)/[^/]*$}); |
66 $WCALC = '.' if ($WCALC eq ''); |
68 $WCALC = '.' if ($WCALC eq ''); |
67 |
69 |
68 require "$WCALC/version.pl"; |
70 require "$WCALC/version.pl"; |
69 require "$ANTS/ants.pl"; |
71 require "$ANTS/ants.pl"; |
70 require "$ANTS/fft.pl"; |
72 require "$ANTS/fft.pl"; |
|
73 require "$ANTS/libstats.pl"; |
71 require "$ANTS/libGMT.pl"; |
74 require "$ANTS/libGMT.pl"; |
72 require "$ANTS/libSBE.pl"; |
75 require "$ANTS/libSBE.pl"; |
73 require "$ANTS/libEOS83.pl"; |
76 require "$ANTS/libEOS83.pl"; |
74 &antsAddParams('LADCP_w_CTD::version',$VERSION); |
77 &antsAddParams('LADCP_w_CTD::version',$VERSION); |
75 |
78 |
76 $antsParseHeader = 0; # usage |
79 $antsParseHeader = 0; # usage |
77 &antsUsage('ai:l:orp:s:v:',1, |
80 &antsUsage('ai:l:orp:s:v:w:',1, |
78 '[-v)erbosity <level[0]>]', |
81 '[-v)erbosity <level[0]>]', |
79 '[use -a)lternate sensor pair]', |
82 '[use -a)lternate sensor pair]', |
80 '[-r)etain all data (no editing)] [allow infinite -o)utliers]', |
83 '[-r)etain all data (no editing)] [allow infinite -o)utliers]', |
81 '[-s)ampling <rate[6Hz]>]', |
84 '[-s)ampling <rate[6Hz]>]', |
82 '[-l)owpass w <cutoff[2s]>]', |
85 '[-l)owpass w_CTD <cutoff[2s]>] [-w)inch-speed <granularity[10s]>]', |
83 '[profile -i)d <id>]', |
86 '[profile -i)d <id>]', |
84 '[-p)lot_basenames <[%03d_w_CTD.ps],[%03d_sspd.ps]>]', |
87 '[-p)lot_basenames <[%03d_w_CTD.ps],[%03d_sspd.ps]>]', |
85 '<SBE CNV file>'); |
88 '<SBE CNV file>'); |
86 |
89 |
87 &antsFloatOpt(\$opt_l,2); # defaults |
90 &antsFloatOpt(\$opt_l,2); # default low-pass cutoff for w_CTD |
88 &antsCardOpt(\$opt_s,6); |
91 &antsCardOpt(\$opt_s,6); # default output samplint rate (Hz) |
|
92 &antsFloatOpt(\$opt_w,10); # winch velocity granularity |
89 &antsCardOpt(\$opt_v,$ENV{VERB}); # support VERB env variable |
93 &antsCardOpt(\$opt_v,$ENV{VERB}); # support VERB env variable |
90 |
94 |
91 $CNVfile = $ARGV[0]; # open CNV file |
95 $CNVfile = $ARGV[0]; # open CNV file |
92 open(F,&antsFileArg()); |
96 open(F,&antsFileArg()); |
93 &antsActivateOut(); # activate ANTS file |
97 &antsActivateOut(); # activate ANTS file |
577 printf(STDERR "\n") if ($opt_v); |
581 printf(STDERR "\n") if ($opt_v); |
578 } else { |
582 } else { |
579 @w_lp = @w; |
583 @w_lp = @w; |
580 } |
584 } |
581 |
585 |
582 #---------------------------------------------------------------------- |
586 #---------------------------------------- |
|
587 # Estimate winch speed |
|
588 #---------------------------------------- |
|
589 |
|
590 print(STDERR "Estimating winch velocity...") if ($opt_v); |
|
591 &antsAddParams('winch_velocity_granularity',$opt_w); |
|
592 |
|
593 my($from_r) = 0; my($to_r); # step 1: bin average in time |
|
594 for (my($from_r)=my($to_r)=0; $from_r<@elapsed; $from_r=$to_r+1) { |
|
595 my($sumw) = $w_lp[2*$from_r]/@w_lp; my($n) = 1; |
|
596 for ($to_r=$from_r+1; $to_r<@elapsed && $elapsed[$to_r]-$elapsed[$from_r]<$opt_w; $to_r++) { |
|
597 $sumw += $w_lp[2*$to_r]/@w_lp; $n++; |
|
598 } |
|
599 $winch[$from_r] = $sumw/$n; |
|
600 } |
|
601 |
|
602 my($pwinch) = $winch[0]; |
|
603 for (my($to_r),my($from_r)=0; $from_r<@elapsed; ) { # step 2: fill after median filtering |
|
604 for ($to_r=$from_r+1; $to_r<@elapsed && !defined($winch[$to_r]); $to_r++) {} |
|
605 my($nwinch) = $to_r<@elapsed ? $winch[$to_r] : $winch[$from_r]; |
|
606 my($winch) = median($pwinch,$winch[$from_r],$nwinch); |
|
607 $pwinch = $winch[$from_r]; $winch[$from_r] = $winch; |
|
608 while (++$from_r < $to_r) { $winch[$from_r] = $winch[$from_r-1]; } |
|
609 } |
|
610 |
|
611 printf(STDERR "\n") if ($opt_v); |
|
612 |
|
613 #---------------------------------------- |
|
614 # Plot Data |
|
615 #---------------------------------------- |
583 |
616 |
584 if (defined($opt_p)) { |
617 if (defined($opt_p)) { |
585 print(STDERR "Plotting data...\n") if ($opt_v); |
618 print(STDERR "Plotting data...\n") if ($opt_v); |
586 my(@pfmt) = split(',',$opt_p); |
619 my(@pfmt) = split(',',$opt_p); |
587 croak("$0: cannot decode -p $opt_p\n") |
620 croak("$0: cannot decode -p $opt_p\n") |
595 GMT_psxy('-W1,coral'); |
628 GMT_psxy('-W1,coral'); |
596 for ($r=0; $r<@w; $r++) { |
629 for ($r=0; $r<@w; $r++) { |
597 printf(GMT "%f %f\n",$elapsed[$r]/60,$w_lp[2*$r]/@w_lp); |
630 printf(GMT "%f %f\n",$elapsed[$r]/60,$w_lp[2*$r]/@w_lp); |
598 GMT_psxy('-W1,SeaGreen') if ($r == $atBtm); |
631 GMT_psxy('-W1,SeaGreen') if ($r == $atBtm); |
599 } |
632 } |
|
633 GMT_psxy('-W1,magenta'); |
|
634 for ($r=0; $r<@w; $r++) { |
|
635 printf(GMT "%f %f\n",$elapsed[$r]/60,$winch[$r]); |
|
636 } |
600 GMT_psbasemap('-Bg60a30f5:"Elapsed Time [min]":/g1a1f0.1:"Vertical Package Velocity [ms@+-1@+]":WeSn'); |
637 GMT_psbasemap('-Bg60a30f5:"Elapsed Time [min]":/g1a1f0.1:"Vertical Package Velocity [ms@+-1@+]":WeSn'); |
601 GMT_unitcoords(); |
638 GMT_unitcoords(); |
602 GMT_pstext('-F+f14,Helvetica,coral+jTR'); print(GMT "0.98 0.98 dc\n"); |
639 GMT_pstext('-F+f14,Helvetica,coral+jBR'); print(GMT "0.98 0.96 dc\n"); |
603 GMT_pstext('-F+f14,Helvetica,SeaGreen+jTR'); print(GMT "0.98 0.94 uc\n"); |
640 GMT_pstext('-F+f14,Helvetica,SeaGreen+jBR'); print(GMT "0.98 0.93 uc\n"); |
|
641 GMT_pstext('-F+f14,Helvetica,magenta+jBR'); print(GMT "0.98 0.89 winch\n"); |
604 GMT_pstext('-F+f14,Helvetica,blue+jTL -N'); |
642 GMT_pstext('-F+f14,Helvetica,blue+jTL -N'); |
605 if (defined($outfile)) { printf(GMT "0.01 1.06 $outfile\n",$id); } |
643 if (defined($outfile)) { printf(GMT "0.01 1.06 $outfile\n",$id); } |
606 else { printf(GMT "0.01 1.06 %03d\n",$id); } |
644 else { printf(GMT "0.01 1.06 %03d\n",$id); } |
607 GMT_end(); |
645 GMT_end(); |
608 |
646 |
616 for ($r=0; $r<@w; $r++) { |
654 for ($r=0; $r<@w; $r++) { |
617 printf(GMT "%f %f\n",$sspd[$r],$depth[$r]); |
655 printf(GMT "%f %f\n",$sspd[$r],$depth[$r]); |
618 GMT_psxy('-W1.5,SeaGreen') if ($r == $atBtm); |
656 GMT_psxy('-W1.5,SeaGreen') if ($r == $atBtm); |
619 } |
657 } |
620 GMT_unitcoords(); |
658 GMT_unitcoords(); |
621 GMT_pstext('-F+f14,Helvetica,coral+jTR'); print(GMT "0.98 0.02 dc\n"); |
659 GMT_pstext('-F+f14,Helvetica,coral+jTR'); print(GMT "0.98 0.02 dc\n"); |
622 GMT_pstext('-F+f14,Helvetica,SeaGreen+jTR'); print(GMT "0.98 0.06 uc\n"); |
660 GMT_pstext('-F+f14,Helvetica,SeaGreen+jTR'); print(GMT "0.98 0.06 uc\n"); |
623 GMT_pstext('-F+f14,Helvetica,blue+jTL -N'); |
661 GMT_pstext('-F+f14,Helvetica,blue+jTL -N'); |
624 if (defined($outfile)) { printf(GMT "0.01 -0.06 $outfile\n",$id); } |
662 if (defined($outfile)) { printf(GMT "0.01 -0.06 $outfile\n",$id); } |
625 else { printf(GMT "0.01 -0.06 %03d\n",$id); } |
663 else { printf(GMT "0.01 -0.06 %03d\n",$id); } |
626 GMT_end(); |
664 GMT_end(); |
627 } |
665 } |
628 |
666 |
629 #---------------------------------------------------------------------- |
667 #---------------------------------------------------------------------- |
630 |
668 |
631 print(STDERR "Writing output...\n") if ($opt_v); |
669 print(STDERR "Writing output...\n") if ($opt_v); |
632 |
670 |
633 @antsNewLayout = ('elapsed','press','temp','cond','depth','salin','sspd','w.raw','w','lat','lon'); |
671 @antsNewLayout = ('elapsed','press','temp','cond','depth','salin','sspd','w_CTD.raw','w_CTD','w_winch','lat','lon'); |
634 for ($r=0; $r<@w; $r++) { |
672 for ($r=0; $r<@w; $r++) { |
635 &antsOut($elapsed[$r],$press[$r],$temp[$r],$cond[$r],$depth[$r],$salin[$r], |
673 &antsOut($elapsed[$r],$press[$r],$temp[$r],$cond[$r],$depth[$r],$salin[$r], |
636 $sspd[$r],$w[$r],$w_lp[2*$r]/@w_lp,$lat[$r],$lon[$r]); |
674 $sspd[$r],$w[$r],$w_lp[2*$r]/@w_lp,$winch[$r], |
|
675 $lat[$r],$lon[$r]); |
637 } |
676 } |
638 |
677 |
639 exit(0); # don't flush @ants_ |
678 exit(0); # don't flush @ants_ |