edit_data.pl
changeset 57 69e39fcb7f41
parent 56 8f120b9f795a
equal deleted inserted replaced
56:8f120b9f795a 57:69e39fcb7f41
     1 #======================================================================
     1 #======================================================================
     2 #                    E D I T _ D A T A . P L 
     2 #                    E D I T _ D A T A . P L 
     3 #                    doc: Sat May 22 21:35:55 2010
     3 #                    doc: Sat May 22 21:35:55 2010
     4 #                    dlm: Fri Jul  9 13:41:36 2021
     4 #                    dlm: Mon Oct 18 21:39:56 2021
     5 #                    (c) 2010 A.M. Thurnherr
     5 #                    (c) 2010 A.M. Thurnherr
     6 #                    uE-Info: 597 63 NIL 0 0 72 0 2 4 NIL ofnI
     6 #                    uE-Info: 55 60 NIL 0 0 72 72 2 4 NIL ofnI
     7 #======================================================================
     7 #======================================================================
     8 
     8 
     9 # HISTORY:
     9 # HISTORY:
    10 #	May 22, 2010: - created
    10 #	May 22, 2010: - created
    11 #	May 24, 2010: - added editSideLobesFromSeabed()
    11 #	May 24, 2010: - added editSideLobesFromSeabed()
    14 #						 interpreted as along-beam, rather than vertical
    14 #						 interpreted as along-beam, rather than vertical
    15 #				  - replaced editPitchRoll by editTilt
    15 #				  - replaced editPitchRoll by editTilt
    16 #	Dec 25, 2010: - adapted to changes in [LADCP_w]
    16 #	Dec 25, 2010: - adapted to changes in [LADCP_w]
    17 #	Aug  3, 2011: - added editTruncRange()
    17 #	Aug  3, 2011: - added editTruncRange()
    18 #	Oct 10, 2011: - added editFalsePositives()
    18 #	Oct 10, 2011: - added editFalsePositives()
    19 #				  - BUG: when earth velocities were edited, all were
    19 #				  - BUG: when Earth velocities were edited, all were
    20 #						 counted, not just those between first and lastBin
    20 #						 counted, not just those between first and lastBin
    21 #	Oct 11, 2011: - moved defaults to [defaults.pl]
    21 #	Oct 11, 2011: - moved defaults to [defaults.pl]
    22 #	Oct 12, 2011: - added &editSurfLayer()
    22 #	Oct 12, 2011: - added &editSurfLayer()
    23 #				  - BUG: editSideLobes() was slightly loose
    23 #				  - BUG: editSideLobes() was slightly loose
    24 #	Oct 15, 2011: - added editWOutliers()
    24 #	Oct 15, 2011: - added editWOutliers()
    44 #	Oct 13, 2017: - BUG: editPPI() only allowed for nominal transducer frequencies
    44 #	Oct 13, 2017: - BUG: editPPI() only allowed for nominal transducer frequencies
    45 #	May  1, 2018: - added editLargeHSpeeds()
    45 #	May  1, 2018: - added editLargeHSpeeds()
    46 #	Nov 17, 2018: - BUG: spurious letter "z" had crept in at some stage
    46 #	Nov 17, 2018: - BUG: spurious letter "z" had crept in at some stage
    47 #	Mar 23, 2021: - updated PPI doc
    47 #	Mar 23, 2021: - updated PPI doc
    48 #	Jul  9, 2021: - added editHighResidualLayers()
    48 #	Jul  9, 2021: - added editHighResidualLayers()
       
    49 #	Sep  1, 2021: - added Sv editing to editHighResidualLayers()
       
    50 #				  - modified sidelobe editing to include instrument tilt
       
    51 #	Oct 15, 2021: - BUG: new sidelobe editing was stupid, because sidelobe
       
    52 #					contamination works in the time domain and is, therefore,
       
    53 #					independent of tilt
       
    54 #	Oct 18, 2021: - BUG: seabed contamination was missing abs() and did not
       
    55 #						 work correctly with missing Sv data
    49 # END OF HISTORY
    56 # END OF HISTORY
    50 
    57 
    51 # NOTES:
    58 # NOTES:
    52 #	- all bins must be edited (not just the ones between $LADCP_firstBin
    59 #	- all bins must be edited (not just the ones between $LADCP_firstBin
    53 #	  and $LADCP_lastBin to allow reflr calculations to use bins outside
    60 #	  and $LADCP_lastBin to allow reflr calculations to use bins outside
   266 		$nrm++;
   273 		$nrm++;
   267 	}
   274 	}
   268 	return $nrm;
   275 	return $nrm;
   269 }
   276 }
   270 
   277 
   271 #======================================================================
   278 #===========================================================================================
   272 # ($nvrm,$nerm) = editSideLobes($fromEns,$toEns,$water_depth)
   279 # ($nvrm,$nerm) = editSideLobes($fromEns,$toEns,$water_depth)
   273 #
   280 #
   274 # NOTES:
   281 # NOTES:
   275 #	1) When this code is executed the sound speed is known. No attempt is made to correct for
   282 #	- When this code is executed the sound speed is known. No attempt is made to correct for
   276 #	   along-beam soundspeed variation, but the soundspeed at the transducer is accounted for.
   283 #	  along-beam soundspeed variation, but the soundspeed at the transducer is accounted for.
   277 #	2) for surface sidelobes, water_depth == undef; surface sidelobes include the
   284 #	- for surface sidelobes, water_depth == undef; surface sidelobes include the
   278 #	   vessel draft
   285 #	  vessel draft
   279 #	- all velocities are counted, even those outside valid bin range,
   286 #	- all velocities are counted, even those outside valid bin range,
   280 #	  because the %age is not reported
   287 #	  because the %age is not reported
   281 #======================================================================
   288 #	- while this filter removes the sidelobe contamination in most profiles
       
   289 #	  there are still profiles with Sv.diff anomalies near the seabed
       
   290 #     (based on SR1b/2004 data); for these editSeabedContamination has been implemented
       
   291 #==========================================================================================
   282 
   292 
   283 sub editSideLobes($$$)
   293 sub editSideLobes($$$)
   284 {
   294 {
   285 	my($fe,$te,$wd) = @_;	# first & last ens to process, water depth for sidelobes near seabed
   295 	my($fe,$te,$wd) = @_;	# first & last ens to process, water depth for sidelobes near seabed
   286 	my($nvrm) = 0;			# of velocities removed
   296 	my($nvrm) = 0;			# of velocities removed
   288 	for (my($e)=$fe; $e<=$te; $e++) {
   298 	for (my($e)=$fe; $e<=$te; $e++) {
   289 		next unless numberp($LADCP{ENSEMBLE}[$e]->{CTD_DEPTH});
   299 		next unless numberp($LADCP{ENSEMBLE}[$e]->{CTD_DEPTH});
   290 		my($range) = defined($wd) ? $wd - $LADCP{ENSEMBLE}[$e]->{CTD_DEPTH} 
   300 		my($range) = defined($wd) ? $wd - $LADCP{ENSEMBLE}[$e]->{CTD_DEPTH} 
   291 								  : $LADCP{ENSEMBLE}[$e]->{CTD_DEPTH} - $vessel_draft;
   301 								  : $LADCP{ENSEMBLE}[$e]->{CTD_DEPTH} - $vessel_draft;
   292 		$range = 0 if ($range < 0);								  
   302 		$range = 0 if ($range < 0);								  
       
   303 		
       
   304 #		from UH code 
   293 		my($sscorr) = $CTD{SVEL}[$LADCP{ENSEMBLE}[$e]->{CTD_SCAN}] / $LADCP{ENSEMBLE}[$e]->{SPEED_OF_SOUND};
   305 		my($sscorr) = $CTD{SVEL}[$LADCP{ENSEMBLE}[$e]->{CTD_SCAN}] / $LADCP{ENSEMBLE}[$e]->{SPEED_OF_SOUND};
   294 		my($goodBins) =   ($range - $sscorr*$LADCP{DISTANCE_TO_BIN1_CENTER}) * cos(rad($LADCP{BEAM_ANGLE}))
   306 		my($firstBadBin) =   ($range - $sscorr*$LADCP{DISTANCE_TO_BIN1_CENTER}) * cos(rad($LADCP{BEAM_ANGLE}))
   295 						/ ($sscorr*$LADCP{BIN_LENGTH})
   307 						/ ($sscorr*$LADCP{BIN_LENGTH})
   296 						- 1.5;
   308 						- 1.5;
   297 
   309 
   298 		my($dirty) = 0;
   310 		my($dirty) = 0;
   299 		for (my($bin)=int($goodBins); $bin<$LADCP{N_BINS}; $bin++) { 	# NB: 2 good bins implies that bin 2 is bad
   311 		for (my($bin)=int($firstBadBin); $bin<$LADCP{N_BINS}; $bin++) { 	
   300 			next unless ($bin>=0 && defined($LADCP{ENSEMBLE}[$e]->{W}[$bin]));
   312 			next unless ($bin>=0 && defined($LADCP{ENSEMBLE}[$e]->{W}[$bin]));
   301 			$dirty = 1;
   313 			$dirty = 1;
   302 			$nvrm++;
   314 			$nvrm++;
   303 			undef($LADCP{ENSEMBLE}[$e]->{W}[$bin]);
   315 			undef($LADCP{ENSEMBLE}[$e]->{W}[$bin]);
       
   316 			debugmsg("sidelobe at range=$range firstBadBin=$firstBadBin ens=$e bin=$bin at CTD depth = $LADCP{ENSEMBLE}[$e]->{CTD_DEPTH}\n");
   304 		}
   317 		}
   305 
   318 
   306 		$nerm += $dirty;
   319 		$nerm += $dirty;
   307 	}
   320 	}
   308 	return ($nvrm,$nerm);
   321 	return ($nvrm,$nerm);
   571     return $nerm;
   584     return $nerm;
   572 }
   585 }
   573 
   586 
   574 #======================================================================
   587 #======================================================================
   575 # $nbrm = editHighResidualLayers($max_lr_res)
   588 # $nbrm = editHighResidualLayers($max_lr_res)
   576 #	- filter applied after gridding
   589 #	- filter applied after depth binning
       
   590 #	- while filter only removes values from profiles, the corresponding
       
   591 #	  samples are not output to .wsamp, because MEDIAN_W is not defined
   577 #	- filter based on observation that profiles of beam-pair residuals
   592 #	- filter based on observation that profiles of beam-pair residuals
   578 #	  are good indicators of bad data, but very noisy
   593 #	  are good indicators of bad data, but very noisy
   579 #	- current version uses estimates in 5-bin-thick layers (200m by
   594 #	- current version uses estimates in 5-bin-thick layers (200m by
   580 #	  default) 
   595 #	  default) 
   581 #	- filter cutoff based on 2021 A20 cruise which crossed region with
   596 #	- filter cutoff based on 2021 A20 cruise which crossed region with
   588 
   603 
   589 	my($nbrm) = 0;
   604 	my($nbrm) = 0;
   590 	for (my($bi)=0; $bi<=$#{$DNCAST{LR_RMS_BP_RESIDUAL}}; $bi++) {
   605 	for (my($bi)=0; $bi<=$#{$DNCAST{LR_RMS_BP_RESIDUAL}}; $bi++) {
   591 		next unless ($DNCAST{LR_RMS_BP_RESIDUAL}[$bi] > $limit);
   606 		next unless ($DNCAST{LR_RMS_BP_RESIDUAL}[$bi] > $limit);
   592 		$DNCAST{MEDIAN_W}[$bi] = $DNCAST{MEDIAN_W12}[$bi] = $DNCAST{MEDIAN_W34}[$bi] = nan;
   607 		$DNCAST{MEDIAN_W}[$bi] = $DNCAST{MEDIAN_W12}[$bi] = $DNCAST{MEDIAN_W34}[$bi] = nan;
       
   608 #		$DNCAST{SV}[$bi] = nan;
   593 		$nbrm++;
   609 		$nbrm++;
   594 	}
   610 	}
   595 	for (my($bi)=0; $bi<=$#{$UPCAST{LR_RMS_BP_RESIDUAL}}; $bi++) {
   611 	for (my($bi)=0; $bi<=$#{$UPCAST{LR_RMS_BP_RESIDUAL}}; $bi++) {
   596 		next unless ($UPCAST{LR_RMS_BP_RESIDUAL}[$bi] > $limit);
   612 		next unless ($UPCAST{LR_RMS_BP_RESIDUAL}[$bi] > $limit);
   597 		$UPCAST{MEDIAN_W}[$bi] = $UPCAST{MEDIAN_W12}[$bi] = $UPCAST{MEDIAN_W34}[$bi] = nan;
   613 		$UPCAST{MEDIAN_W}[$bi] = $UPCAST{MEDIAN_W12}[$bi] = $UPCAST{MEDIAN_W34}[$bi] = nan;
       
   614 #		$UPCAST{SV}[$bi] = nan;
   598 		$nbrm++;
   615 		$nbrm++;
   599 	}
   616 	}
   600     return $nbrm;
   617     return $nbrm;
   601 }
   618 }
   602 
   619 
   603 #======================================================================
   620 #======================================================================
       
   621 # $nbrm = editSeabedContamination($max_lr_res)
       
   622 #	- filter applied after depth binning
       
   623 #	- while filter only removes values from profiles, the corresponding
       
   624 #	  samples are not output to .wsamp, because MEDIAN_W is not defined
       
   625 #	- filter based on SR1b/2004 observation that some profiles of Sv
       
   626 #	  show clear anomalies near the seabed
       
   627 #	- anomalies are easily detected in dSv/dz plots, with
       
   628 #	  $seabed_contamination_Sv_grad_limit = 0.1 being a suitable limiting 
       
   629 #	  value for WH300 ADCPs
       
   630 #======================================================================
       
   631 
       
   632 sub editSeabedContamination($)
       
   633 {
       
   634 	my($limit) = @_;
       
   635 	my($nbrm) = 0;
       
   636 
       
   637 	for (my($bi)=$#{$DNCAST{SV}}; $bi>0; $bi--) {
       
   638 		next unless numberp($DNCAST{SV}[$bi]) && numberp($DNCAST{SV}[$bi-1]);
       
   639 		last if (abs($DNCAST{SV}[$bi]-$DNCAST{SV}[$bi-1])/$opt_o < $limit);
       
   640 		$DNCAST{MEDIAN_W}[$bi] = $DNCAST{MEDIAN_W12}[$bi] = $DNCAST{MEDIAN_W34}[$bi] = nan;
       
   641 		$DNCAST{SV}[$bi] = nan;
       
   642 		$nbrm++;
       
   643     }
       
   644 
       
   645 	for (my($bi)=$#{$UPCAST{SV}}; $bi>0; $bi--) {
       
   646 		next unless numberp($UPCAST{SV}[$bi]) && numberp($UPCAST{SV}[$bi-1]);
       
   647 		last if (abs($UPCAST{SV}[$bi]-$UPCAST{SV}[$bi-1])/$opt_o < $limit);
       
   648 		$UPCAST{MEDIAN_W}[$bi] = $UPCAST{MEDIAN_W12}[$bi] = $UPCAST{MEDIAN_W34}[$bi] = nan;
       
   649 		$UPCAST{SV}[$bi] = nan;
       
   650 		$nbrm++;
       
   651     }
       
   652 
       
   653     return $nbrm;
       
   654 }
       
   655 
       
   656 #======================================================================
   604 
   657 
   605 1;
   658 1;