# HG changeset patch # User A.M. Thurnherr # Date 1308331984 14400 # Node ID 16726a31a399a2c79ddf07d937ebd45ddd64a86e # Parent 54222c82435fdf157e23f566ba3cf66b6c611d27 pre IWISE diff --git a/LADCPintsh b/LADCPintsh --- 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: Fri Dec 10 04:57:50 2010 +# dlm: Wed Feb 16 15:13:46 2011 # (c) 2010 A.M. Thurnherr & E. Firing -# uE-Info: 207 1 NIL 0 0 72 2 2 4 NIL ofnI +# uE-Info: 401 0 NIL 0 0 72 2 2 4 NIL ofnI #====================================================================== $antsSummary = 'integrate LADCP shear'; @@ -28,6 +28,7 @@ # Oct 24, 2010: - fix spuriously small variances that can occur for BT velocities based on very small # samples (i.e. primarily when chosing a small -r) # Dec 9, 2010: - allowed for empty BT file +# Feb 16, 2011: - BUG: gaps in shear data were not handled correctly in baroclinic solution ($ANTS) = (`which list` =~ m{^(.*)/[^/]*$}); require "$ANTS/ants.pl"; @@ -378,16 +379,26 @@ unless (defined($refU)) { my($sumU,$sumV,$sumW,$dc_sumU,$dc_sumV,$dc_sumW,$uc_sumU,$uc_sumV,$uc_sumW); + my($nSumVel,$dc_nSumVel,$uc_nSumVel); for (my($r)=0; $r<@depth; $r++) { - $sumU += $u[$r]; $sumV += $v[$r]; $sumW += $w[$r]; - $dc_sumU += $dc_u[$r]; $dc_sumV += $dc_v[$r]; $dc_sumW += $dc_w[$r]; - $uc_sumU += $uc_u[$r]; $uc_sumV += $uc_v[$r]; $uc_sumW += $uc_w[$r]; + if (numberp($u[$r])) { + $nSumVel++; + $sumU += $u[$r]; $sumV += $v[$r]; $sumW += $w[$r]; + } + if (numberp($dc_u[$r])) { + $dc_nSumVel++; + $dc_sumU += $dc_u[$r]; $dc_sumV += $dc_v[$r]; $dc_sumW += $dc_w[$r]; + } + if (numberp($uc_u[$r])) { + $uc_nSumVel++; + $uc_sumU += $uc_u[$r]; $uc_sumV += $uc_v[$r]; $uc_sumW += $uc_w[$r]; + } } - $refU = $sumU / @depth; $refV = $sumV / @depth; $refW = $sumW / @depth; - $dc_refU = $dc_sumU / @depth; $dc_refV = $dc_sumV / @depth; $dc_refW = $dc_sumW / @depth; - $uc_refU = $uc_sumU / @depth; $uc_refV = $uc_sumV / @depth; $uc_refW = $uc_sumW / @depth; + $refU = $sumU / $nSumVel; $refV = $sumV / $nSumVel; $refW = $sumW / $nSumVel; + $dc_refU = $dc_sumU / $dc_nSumVel; $dc_refV = $dc_sumV / $dc_nSumVel; $dc_refW = $dc_sumW / $dc_nSumVel; + $uc_refU = $uc_sumU / $uc_nSumVel; $uc_refV = $uc_sumV / $uc_nSumVel; $uc_refW = $uc_sumW / $uc_nSumVel; } for (my($r)=0; $r<@depth; $r++) { # reference velocities diff --git a/LADCPproc b/LADCPproc --- 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: Mon Jan 10 13:53:54 2011 +# dlm: Wed Jun 15 15:47:25 2011 # (c) 2010 A.M. Thurnherr & E. Firing -# uE-Info: 382 0 NIL 0 0 72 2 2 4 NIL ofnI +# uE-Info: 479 55 NIL 0 0 72 2 2 4 NIL ofnI #====================================================================== $antsSummary = 'process LADCP data to get shear, time series'; @@ -32,8 +32,11 @@ # Dec 10, 2010: - change -w) default to 120s # - changed nshear output to 0 from nan when there are no samples # Dec 27, 2010: - changed sign of -l to accept lag output from [LADCP_w] -# Jan 10, 2010: - -o => -k added new -o +# Jan 10, 2011: - -o => -k added new -o # - added code to fill CTD sound vel gaps +# Jan 22, 2011: - added -g) lat,lon +# - added -c)ompass corr +# Jun 15, 2011: - added mean backscatter profile to default output ($ANTS) = (`which list` =~ m{^(.*)/[^/]*$}); ($PERL_TOOLS) = (`which mkProfile` =~ m{^(.*)/[^/]*$}); @@ -52,17 +55,18 @@ require "$PERL_TOOLS/RDI_Utils.pl"; $antsParseHeader = 0; -&antsUsage('24ab:df:i:kl:n:o:ps:t:w:',2, +&antsUsage('24ab:c:df:g:i:kl:n:o:ps:t:w:',2, '[use -2)dary CTD sensor pair]', '[require -4)-beam LADCP solutions]', - '[-s)etup ]', + '[-s)etup ] [-g)ps ]', + '[-c)ompass-corr ]', '[enable -p)PI editing]', '[-o)utput grid ]', '[-i)nitial LADCP time lag ]', '[-l)ag LADCP ] [auto-lag -w)indow ] [-n) ] [-f)lag ] [-b)ottom-track ]', - ' [-a)coustic backscatter profiles] [bottom-trac-k) profs]', + ' [per-bin -a)coustic backscatter profiles] [bottom-trac-k) profs]', ' '); $RDI_Coords::minValidVels = 4 if ($opt_4); @@ -78,6 +82,19 @@ &antsFloatOpt($opt_i); &antsCardOpt($opt_o); +if (defined($opt_g)) { + ($CTD{lat},$CTD{lon}) = split(',',$opt_g); + croak("$0: cannot decode -g $opt_g\n") + unless numberp($CTD{lat}) && numberp($CTD{lon}); +} + +if (defined($opt_c)) { + ($CC_offset,$CC_cos_fac,$CC_sin_fac) = split(',',$opt_c); + croak("$0: cannot decode -c $opt_c\n") + unless numberp($CC_offset) && numberp($CC_cos_fac) && numberp($CC_sin_fac); +} + + $LADCP_file = &antsFileArg(); $CTD_file = &antsFileArg(); @@ -209,10 +226,18 @@ &antsAddParams('3_beam_solutions',"$RDI_Coords::threeBeam_1 $RDI_Coords::threeBeam_2 $RDI_Coords::threeBeam_3 $RDI_Coords::threeBeam_4"); } } elsif ($LADCP{EARTH_COORDINATES}) { - printf(STDERR "\n\t\tcorrecting for magnetic declination of %.1f deg...\n",$magdec) - if ($opt_d); + if ($opt_d) { + if ($opt_c) { + printf(STDERR "\n\t\tcalculating tilt and correcting for compass error and magnetic declination of %.1f deg...\n",$magdec); + } else { + printf(STDERR "\n\t\tcalculating tilt and correcting for magnetic declination of %.1f deg...\n",$magdec); + } + } for (my($ens)=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) { $LADCP{ENSEMBLE}[$ens]->{TILT} = &angle_from_vertical($LADCP{ENSEMBLE}[$ens]->{PITCH},$LADCP{ENSEMBLE}[$ens]->{ROLL}); + my($hdg) = rad($LADCP{ENSEMBLE}[$ens]->{HEADING}); + $LADCP{HEADING_BIAS} = ($CC_offset + $CC_cos_fac*cos($hdg) + $CC_sin_fac*sin($hdg)) - $magdec + if ($opt_c); for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) { @{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]} = velApplyHdgBias(\%LADCP,$ens,@{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]}); @@ -421,7 +446,7 @@ @antsNewLayout = ('depth','dc_nshear','dc_u_z','dc_u_z.sig','dc_v_z','dc_v_z.sig','dc_w_z','dc_w_z.sig', 'uc_nshear','uc_u_z','uc_u_z.sig','uc_v_z','uc_v_z.sig','uc_w_z','uc_w_z.sig', - 'nshear','u_z','u_z.sig','v_z','v_z.sig','w_z','w_z.sig'); + 'nshear','u_z','u_z.sig','v_z','v_z.sig','w_z','w_z.sig','Sv','Sv.n'); $commonParams = $antsCurParams; &antsAddParams('ubin_start',$ubin_start,'ubin_end',$ubin_end, # record processing params @@ -434,7 +459,8 @@ 'max_shdev',$max_shdev,'max_shdev_sum',$max_shdev_sum, 'water_depth',round($water_depth),'water_depth.sig',round($sig_water_depth), 'min_hab',round($min_hab), - 'clip_margin',$clip_margin,'first_clip_bin',$first_clip_bin); + 'clip_margin',$clip_margin,'first_clip_bin',$first_clip_bin, + 'Svbin_start',$Svbin_start,'Svbin_end',$Svbin_end); for (my($gi)=0; $gi<@ush_mu; $gi++) { &antsOut(depthOfGI($gi), # depth in center of bin @@ -446,10 +472,12 @@ $uc_ush_mu[$gi],$uc_ush_sig[$gi], $uc_vsh_mu[$gi],$uc_vsh_sig[$gi], $uc_wsh_mu[$gi],$uc_wsh_sig[$gi], - $sh_n[$gi], + $sh_n[$gi], # combined $ush_mu[$gi],$ush_sig[$gi], $vsh_mu[$gi],$vsh_sig[$gi], $wsh_mu[$gi],$wsh_sig[$gi], + $nSv_prof[$gi]?$sSv_prof[$gi]/$nSv_prof[$gi]:nan, + $nSv_prof[$gi], ); } @@ -458,7 +486,7 @@ #---------------------------------------------------------------------- if (defined($opt_a)) { - print(STDERR "Writing acoustic backscatter profiles..."); + print(STDERR "Writing per-bin acoustic backscatter profiles..."); for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) { my($fn) = sprintf("bin%02d.Sv",$bin); diff --git a/LADCPproc.backscatter b/LADCPproc.backscatter --- a/LADCPproc.backscatter +++ b/LADCPproc.backscatter @@ -1,15 +1,17 @@ #====================================================================== # L A D C P P R O C . B A C K S C A T T E R # doc: Wed Oct 20 13:02:27 2010 -# dlm: Fri Dec 10 00:35:45 2010 +# dlm: Wed Jun 15 15:31:43 2011 # (c) 2010 A.M. Thurnherr -# uE-Info: 12 52 NIL 0 0 72 2 2 4 NIL ofnI +# uE-Info: 67 0 NIL 0 0 72 2 2 4 NIL ofnI #====================================================================== # HISTORY: # Oct 20, 2010: - created # Dec 10, 2010: - BUG: backscatter above sea surface made code bomb # when run with uplooker data +# Jun 15, 2011: - added calculation of backscatter profiles from +# subset of bins #---------------------------------------------------------------------- # Volume Scattering Coefficient, following Deines (IEEE 1999) @@ -67,23 +69,30 @@ my($gi) = int(&depthOfBin($ens,$bin) / $GRID_DZ); next if ($gi < 0); my($range_to_bin) = &dzToBin($ens,$bin) / cos(rad($LADCP{BEAM_ANGLE})); - $sSv[$gi][$bin] += Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP}, - $LADCP{TRANSMITTED_PULSE_LENGTH}, - $Er[0],$range_to_bin, - $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][0])/4 + - Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP}, - $LADCP{TRANSMITTED_PULSE_LENGTH}, - $Er[1],$range_to_bin, - $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][1])/4 + - Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP}, - $LADCP{TRANSMITTED_PULSE_LENGTH}, - $Er[2],$range_to_bin, - $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][2])/4 + - Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP}, - $LADCP{TRANSMITTED_PULSE_LENGTH}, - $Er[3],$range_to_bin, - $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][3])/4; + my($Sv) = Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP}, + $LADCP{TRANSMITTED_PULSE_LENGTH}, + $Er[0],$range_to_bin, + $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][0])/4 + + Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP}, + $LADCP{TRANSMITTED_PULSE_LENGTH}, + $Er[1],$range_to_bin, + $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][1])/4 + + Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP}, + $LADCP{TRANSMITTED_PULSE_LENGTH}, + $Er[2],$range_to_bin, + $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][2])/4 + + Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP}, + $LADCP{TRANSMITTED_PULSE_LENGTH}, + $Er[3],$range_to_bin, + $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][3])/4; + + $sSv[$gi][$bin] += $Sv; $nSv[$gi][$bin]++; + + if ($bin>=$Svbin_start && $bin<=$Svbin_end) { + $sSv_prof[$gi] += $Sv; + $nSv_prof[$gi]++; + } } } } diff --git a/LADCPproc.defaults b/LADCPproc.defaults --- 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 Dec 9 19:40:05 2010 +# dlm: Wed Jun 15 15:15:41 2011 # (c) 2010 A.M. Thurnherr -# uE-Info: 40 31 NIL 0 0 72 0 2 4 NIL ofnI +# uE-Info: 23 48 NIL 0 0 72 0 2 4 NIL ofnI #====================================================================== # default parameters for [mkShearProf] @@ -20,6 +20,7 @@ # HISTORY: # Sep 17, 2010: - created # Dec 9, 2010: - added doc for ASCII CTD file support +# Jun 15, 2011: - added Svbin_start, Svbin_end #---------------------------------------------------------------------- # ASCII CTD file support @@ -44,6 +45,23 @@ # $CTD_ASCII_badval = 999; #---------------------------------------------------------------------- +# Backscatter Profile Parameters +#---------------------------------------------------------------------- + +# The output profile of volume-scattering coefficient is calculated +# from a subset of the bins only. This is based on the observation, +# from on a single P403 profile, that the volume-scattering coefficient +# correction of Deines (IEEE 1999) yield similar results only for +# bins 3-9. + +# The seabed-finding routine, on the other hand, uses acoustic +# backscatter from all bins, as it is possible that the seabed is only +# seen in the farthest bins. + +$Svbin_start = 3; +$Svbin_end = 9; + +#---------------------------------------------------------------------- # Shear Processing Parameters #---------------------------------------------------------------------- diff --git a/LADCPproc.loadCTD b/LADCPproc.loadCTD --- a/LADCPproc.loadCTD +++ b/LADCPproc.loadCTD @@ -1,16 +1,17 @@ #====================================================================== # L A D C P P R O C . L O A D C T D # doc: Thu Dec 9 18:39:01 2010 -# dlm: Mon Jan 10 12:13:54 2011 +# dlm: Sat Jan 22 21:40:12 2011 # (c) 2010 A.M. Thurnherr -# uE-Info: 31 30 NIL 0 0 72 2 2 4 NIL ofnI +# uE-Info: 174 0 NIL 0 0 72 2 2 4 NIL ofnI #====================================================================== # HISTORY: # Dec 9, 2010: - exported from LADCPproc # - added support for ASCII files # Dec 16, 2010: - BUG cnv read did not work any more -# Jan 10, 2010: - added code to skip ANTS header +# Jan 10, 2011: - added code to skip ANTS header +# Jan 22, 2011: - adapted to new -g sub readCTD_ASCII($$) { @@ -19,8 +20,10 @@ croak("$fn: unknown pressure field\n") unless defined($CTD_ASCII_press_field); croak("$fn: unknown temperature field\n") unless defined($CTD_ASCII_temp_field); croak("$fn: unknown salinity field\n") unless defined($CTD_ASCII_salin_field); - croak("$fn: unknown latitude field\n") unless defined($CTD_ASCII_lat_field); - croak("$fn: unknown longitude field\n") unless defined($CTD_ASCII_lon_field); + unless (numberp($dtaR->{lat})) { + croak("$fn: unknown latitude field\n") unless defined($CTD_ASCII_lat_field); + croak("$fn: unknown longitude field\n") unless defined($CTD_ASCII_lon_field); + } $CTD_ASCII_badval = 9e99 unless defined($CTD_ASCII_badval); @@ -33,7 +36,7 @@ push(@{$dtaR->{press}},($rec[$CTD_ASCII_press_field-1] == $CTD_ASCII_badval) ? nan : $rec[$CTD_ASCII_press_field-1]); push(@{$dtaR->{temp}}, ($rec[$CTD_ASCII_temp_field-1] == $CTD_ASCII_badval) ? nan : $rec[$CTD_ASCII_temp_field-1]); push(@{$dtaR->{salin}},($rec[$CTD_ASCII_salin_field-1] == $CTD_ASCII_badval) ? nan : $rec[$CTD_ASCII_salin_field-1]); - unless ($rec[$CTD_ASCII_lat_field-1] == $CTD_ASCII_badval) { + unless (!defined($CTD_ASCII_lat_field) || $rec[$CTD_ASCII_lat_field-1] == $CTD_ASCII_badval) { $nPos++; $sumLat += $rec[$CTD_ASCII_lat_field-1]; $sumLon += $rec[$CTD_ASCII_lon_field-1]; diff --git a/TODO b/TODO --- a/TODO +++ b/TODO @@ -1,9 +1,9 @@ ====================================================================== T O D O doc: Thu Oct 14 09:15:02 2010 - dlm: Mon Jan 10 13:29:13 2011 + dlm: Wed Jun 15 15:16:20 2011 (c) 2010 A.M. Thurnherr - uE-Info: 17 51 NIL 0 0 72 3 2 4 NIL ofnI + uE-Info: 40 0 NIL 0 0 72 3 2 4 NIL ofnI ====================================================================== =LADCPintsh= @@ -36,6 +36,7 @@ - implement TILT_BIT (others?) - calculate real acoustic backscatter profile + - currently, profiles are arbitarily calculated from bins 3-9 - clean up LADCPproc.UHcode: - remove weirdnesses