goodends.m

一个研究声多普勒计程仪很好的工具箱

function [minens, maxens, nens, trimFile] = goodends(theBeamFile,theMaskFile,trimFile)

%function [minens, maxens, nens, trimFile] = goodends(theBeamFile,theMaskFile,trimFile)
%Make an educated guess at the first good ensemble and the last good.
%This function performs a series of checks at the ends of the record:
%		1. If the number of good bins is greater than 75%
%		2. If the pitch and roll vary by less than 2 degrees
%		3. Based on the above criteria there are four possible asnwers, 
%			one for each of the four beams.  AS a last check see if one of those
%			answers also matches the recovery and deployment dates.
%			
%Input:
%	theBeamFile = the masked netcdf data file in beam coordinates
%	theMaskFile = the mask used to screen the data (*.msk),
%						**must be the same dimensions as theBeamFile
%	trimFile = the shortened file to be produced by cutting
%		at the beginning and end of the ensemble record
%
%Outputs:
%	minens = minimum ensemble used to 
%	trimFile = the shortened file to be produced by cutting
%		at the beginning and end of the ensemble record

%Requires additional functions
%	Netcdf toolbox (including nctrim)
%	The following RPS stuff (written by R. Signell)	
%		gregorian.m
%		julian.m

% Written by Jessica M. Cote
% for the U.S. Geological Survey
% Coastal and Marine Geology Program
% Woods Hole, MA
% http://woodshole.er.usgs.gov/
% Please report bugs to [e-mail]
%
%version 1.0
% updated 05-Feb-2003 (ALR) correct comma separator error in versions 6.2 and higher
% updated 21-Sep-2001 - revised trimming of out of water data for Earth and added 
%       more comments to explain the code (ALR)
% updated 28-Dec-2000 08:53:22 - added line feeds to history attribute (ALR)
% updated 27-Dec-2000 14:50:19 - trims recovery time/date more accurately for all data (ALR)
% updated 19-Dec-2000 12:36:57 - trims Earth data more accurately (ALR)
% updated 30-Jun-2000 11:11:50 - better for Earth
% updated 13-Jun-2000 10:53:45
% updated 02-Feb-2000 09:34:12 - runs batch
% updated 12-Jan-2000 09:38:01 - create catch for out of range dates
% updated 15-Oct-1999 09:31:48 

ncquiet

%tell us what function is running
Mname=mfilename;
disp('')
disp([ Mname ' is currently running']);

if nargin < 1, help(mfilename), end
if nargin< 1, theBeamFile = ''; end
if nargin < 2, theMaskFile = ''; end
if nargin < 3, trimFile = ''; end

if isempty(theBeamFile), theBeamFile = '*'; end
if isempty(theMaskFile), theMaskFile = '*'; end
if isempty(trimFile), trimFile = '*'; end

% Get ADCP Beam filename.
if any(theBeamFile == '*')
	[theFile, thePath] = uigetfile(theBeamFile, 'Select ADCP Netcdf File:');
	if ~any(theFile), return, end

	if thePath(end) ~= filesep, thePath(end+1) = filesep; end
	theBeamFile = [thePath theFile];
end
g=netcdf(theBeamFile,'nowrite');

[path,name,ext,ver]=fileparts(theBeamFile);
suggest=[name(1:end-1) 'T.cdf'];


% Get trimmed ADCP filename.
if any(trimFile == '*')
	[theFile, thePath] = uiputfile(suggest, 'Save the Trimmed ADCP File As:');
	if ~any(theFile), return, end
	if thePath(end) ~= filesep, thePath(end+1) = filesep; end
	trimFile = [thePath theFile];
end

%get some information from the data file
dbin=g('bin');
nbins=dbin(:);
ens = g('ensemble');
nens = ens(:);
	%if only 1 ensemble skip to end
   if nens <= 2
      disp('ensemble trimming was skipped due to insufficient ensembles')
      if isunix
      	eval(['!cp  ' theBeamFile ' ' trimFile])
		elseif any(findstr(lower(computer), 'pcwin')) | isVMS
      	eval(['!copy  ' theBeamFile ' ' trimFile])
		elseif any(findstr(lower(computer), 'mac')) & ...
			exist('aduplicate') == 2
			feval('aduplicate', theBeamFile, trimFile)
		else
			fcopy(theBeamFile, trimFile)
      end

      thecomment=sprintf('%s\n','No ensembles were trimmed by goodends.m');
   	history(trimFile,thecomment);
      
      %the required outputs
      minens = nens(1);
      maxens = nens(end);
return
end

coord = g.transform (:);

time=g{'TIM'}(:);
%get recovery and deployment dates and find ensemble where this occurs
ddate=g.Deployment_date(:)
if ~isempty(ddate)
	[y, m, d] = datevec(ddate);
	jr = julian(y,m,d,0);
   ddens = find(time == jr);
   if isempty(ddens)
      %the ADCP may not have been started on the hour
      tim1=time(1);
      start_time = tim1 - floor(tim1);
      newjr = (jr-1) + start_time;
      ddens = find(time == newjr);
      
      if isempty(ddens)	% if still empty
         ddens = 1;
      end
      
   end
else
   ddens = 1;
end

rdate=g.Recovery_date(:)
if ~isempty(rdate)
	[y, m, d] = datevec(rdate);
   jr = julian(y,m,d,0);
   rdens = find(time == jr);
   
   if isempty(rdens)
      %the ADCP may not have been started on the hour
      newjr = (jr+1) - start_time;
     	rdens = find(time == newjr);
   	if isempty(rdens)
      	goodjr = find(floor(time) == jr);
         if ~isempty(goodjr) 
            rdens = goodjr(end);
      	else   
            rdens = nens;
         end
      end 
   end
else
   rdens = nens;
end


% The section below checks the data for bad velocites in order to find the first good ensemble when the ADCP was
%   actually in the water, and the last good ensemble when the ADCP was taken out of the water.
% It checks 'Beam' data in the first portion, and 'Earth' data in the second portion. 
% See each section on how the checking procedure works.

% The Beam data have a mask file which will be used for the check of 'in-the-water-data'
switch coord
case 'BEAM'
   
% Get ADCP Mask filename
if any(theMaskFile == '*')
	[theFile, thePath] = uigetfile(theMaskFile, 'Select ADCP Mask File:');
	if ~any(theFile), return, end

	if thePath(end) ~= filesep, thePath(end+1) = filesep; end
	theMaskFile = [thePath theFile];
end
f=netcdf(theMaskFile,'nowrite');

% Pulls out the masked velocity beam data into 4 cell arrays, one for each beam
for k = 1:4
	vmask{k} = f{['vel' int2str(k)]};
end

CHUNK=100;  % For use below, to break up processing into 'chunks' of 100 for easier/quicker processing
inens=ones(4,1);    % For use below, to store the first good ensemble for each beam
oxens=ones(4,1);    % For use below, to store the last good ensemble for each beam

% Creates an ncvariable 'p' for each cell array of masked beam data
for k = 1:4;
   theVarname = ['vel' int2str(k)];
   disp([' checking ' theVarname ' ...']);
	p = vmask{k};  % new ncvariable for each masked beam data
	m= size(p,1);  % size of each beam (i.e. number of ensembles)
    m=m(1); 
    i = 0;
    t=zeros(m,1);  % temp array of zeros for good/bad data flags for each beam and ensemble to be stored
  
      % Below checks velocity data for bad ensembles by checking the maksed data for ones (bad) and zeros (good).
      % Then sums each column of beam data and if the sum is greater than 25% of the total number of depth bins
      %     the beam for that ensemble is marked with a '1' to designate a bad ensemble for that beam.  
      % For example, the column has a length of 20 (meaning 20 depth bins), and the sum for the column is 10, meaning
      %     that 10 bins were bad velocities, than the Beam data for that ensemble would be marked bad.
      % This check is done in 100 ensemble increments for easier/quicker processing
      while i < m
		 	j = i+1:min(i+CHUNK,nens); % Check data in 'chunks' of 100 ensembles
       	    q=t(j,:);   
            s = sum(p(j,:),2);  % Sums the column of masked data
       	    bad = s > 0.25*nbins;	% if # of bad bins is greater than 25%, it is marked with a '1' for 'bad'
        	q =  q | bad;  
       	    t(j,:)= q;    % Fills in the temp array with flags for 'good' and 'bad' data
       	    i = i + CHUNK;  
      end
      
      id=find(t==0);  % Finds the good velocity data for each beam, for each ensemble
      if isempty(id)
         id = NaN;
      end
      
      % The first/last good ensemble of velocity data (inens/oxens) for each beam is saved here
      inens(k,:)=id(1);
      oxens(k,:)=id(end);
      
end	%for k=1:4 loop
   
% Will give conservative estimates for the min and max ensembles using the first/last good ensembles of velocity