trimbins.m

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

function [MSL, Dstd] = trimbins(numRawFile,rawdata1,rawdata2,trimFile,MSL,Dstd,ADCP_offset,progPath)

%function [MSL, Dstd] = trimbins(numRawFile,rawdata1,rawdata2,trimFile,MSL,Dstd,ADCP_offset)
%Modifies the data file so the bins all fall below mean sea level(MSL) plus 
%the tidal range using a MSL provided by the user or from RDsurface.m,
%a Matlab function that runs the Dos version of the RDI surface program.
%
%INPUTS:
%	rawdata =  RDI ADCP binary output file (*.000)
%	trimFile = use the file create by trimming the bad ensembles in goodends.m if 
%					going through runADCP prcessing steps OR any ADCP file that is in netcdf	
%	ADCP_offset = is the an attribute of the depth variable, 
%		called xducer_offset_from_bottom in the Netcdf files
%		If not given, the program requests a Netcdf file to find it.
%	progPath = the full path for the RDI surface program on your computer
%
%
%If the following variables are defined than RDsurface.m will not be run
%	MSL = same as output variable
%	Dstd = same as output
%
%OUTPUTS:
%	MSL = mean sea level based on RDI surface output
%	Dstd = standard of deviation to give an approximate tidal variation

%sub-functions:
%	Netcdf toolbox
%	RDsurface.m
%	DOS version of surface.exe

% 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 03-Apr-2003 (ALR) clarified note for Height variable
% updated 16-Aug-2001 If using surface.exe, uses 94% of MSL + 1/2 tidal range 
%		to trim surface bins (ALR) 
% updated 10-Jul-2001 added ability to use one OR two binary raw ADCP files (ALR)
% updated 28-Dec-2000 09:28:32 - added linefeeds to comment/history attribute (ALR)
% updated 12-Jan-2000 11:07:35
%	-runs with batch and has all revisions
% updated 22-Oct-1999 14:40:41  
% updated 15-Oct-1999 16:49:45

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


if nargin < 1, numRawFile = ''; end
if nargin < 2, rawdata1 = '';, end
if nargin < 3, rawdata2 = '';, end
if nargin < 4, trimFile = '';, end
if nargin < 7, ADCP_offset = '';, end
if nargin < 8, progPath = '';, end

if isempty(numRawFile)
   numRawFile = menu('How many bianry files used?','1','2');
end
if isempty(rawdata1), rawdata1 = '*'; end
if isempty(rawdata2), rawdata2 = '*'; end
if isempty(trimFile), trimFile = '*'; end
%ADCP_offset will be taken care of in the RDsurface program if isempty.

%Get the files in case the inputs are empty
% Get ADCP raw data filename.
switch numRawFile
case 1
   if any(rawdata1 == '*')
      [theFile1, thePath1] = uigetfile(rawdata1, 'Select Binary ADCP File:');
      if ~any(theFile1), return, end
      if thePath1(end) ~= filesep, thePath1(end+1) = filesep; end
      rawdata1 = [thePath1 theFile1];
   end
end

switch numRawFile
case 2
   if any(rawdata1 == '*')
      [theFile1, thePath1] = uigetfile(rawdata1, 'Select 1st Binary ADCP File:');
      if ~any(theFile1), return, end
      if thePath1(end) ~= filesep, thePath1(end+1) = filesep; end
      rawdata1 = [thePath1 theFile1];
   end
   
   if any(rawdata2 == '*')
      [theFile2, thePath2] = uigetfile(rawdata2, 'Select 2nd Binary ADCP File:');
      if ~any(theFile2), return, end
      if thePath2(end) ~= filesep, thePath2(end+1) = filesep; end
      rawdata2 = [thePath2 theFile2];
   end
end


% Get Ensemble trimmed ADCP data file.
if any(trimFile == '*')
	[theFile, thePath] = uigetfile(trimFile, 'Select Ensemble trimmed ADCP File:');
	if ~any(theFile), return, end
	if thePath(end) ~= filesep, thePath(end+1) = filesep; end
	trimFile = [thePath theFile];
end

%open the netcdf file to trim and get some info
f=netcdf(trimFile,'write');
if isempty(f), return, end;
ADCP_offset = f{'D'}.xducer_offset_from_bottom(:);
depth=f{'D'};
theVars=var(f);
bad_num = fillval(depth);

B = f('bin'); 
theBinNum = B(:);
ensembles = f{'Rec'}(:);
EnsNum = length(ensembles);

if isempty(theBinNum)
    disp(' ## the number of bins not found.')
    close(f)
    return
end


%get the mean sea level and tidal range for the given ensemble numbers
%[MSL,Dstd,Dout] = RDsurface(numRawFile,rawdata1,rawdata2,ADCP_offset,ensembles,progPath,DepthFile);
if nargin < 5
   if running(batch)
      progPath = get(batch);     
      eval(['progPath = ' progPath ';'])
	end
	[MSL,Dstd,Dsurf] = rdsurface(numRawFile,rawdata1,rawdata2,ADCP_offset,ensembles,progPath);
		%to append history later
      thecomment=sprintf('%s\n','Bins were trimmed by trimBins.m using 94% of the RDI surface output.');
else
   disp('User input water_depth and tidal variation')
   Dsurf=[ ];
   %to append history later
   thecomment=sprintf('%s\n','Bins were trimmed by trimBins.m based on user input depth information.');
end


%for files with _FillValue
if ~isempty(bad_num)
   D = autonan(depth,1);
   depth=D(:);
   
else
bad_num=f.VAR_FILL
   
   if ischar(bad_num)==1;
   		bad_num=str2num(bad_num);
   	end
      
   depth=depth(:);
   idgood=find(depth < bad_num & depth > 0);
   depth=depth(idgood);
   clear idgood  
   
end


%find the bins1 that fall below the given sea level including tidal fluctuation   
disp(' ')
disp('Finding the good bins');
%goodBins = find(depth <= MSL+Dstd); Only trims bins out of water
% Because data within 6% of the surface is contaminated, we are trimming those bins as well.
% Reference: Principles of Operation: A Practical Primer (for ADCPs) pg. 38
goodBins=find(depth <= (0.94*(MSL+(Dstd))));

if isempty(goodBins)  
   disp('## No bins were found within the specified depth range');
   close(f)
   return
elseif isequal(goodBins,length(depth));
   disp('No bad bins were found, file is unmodified');
   f{'D'}.water_depth(:) =MSL;
   close(f)
   return
 
end

%Check to make sure that everything is on order
%The goodbins should run from 1(closets to the ADCP) to the depth.
goodBins=sort(goodBins);
%make sure indices are unique
Bd = find(diff(goodBins) == 0);
if any(Bd) 
   goodBins(Bd) = []; 
	disp('Warning: Empty bins were found')   
end  


% compute bin locations
offset=f{'D'}.xducer_offset_from_bottom(:);
bin1=f{'D'}.center_first_bin(:);
binsize=f{'D'}.bin_size(:);
bincnt=length(goodBins);
EndD = (((bincnt-1)*binsize)+bin1) + offset;

%resize the dimension bins if all ok
if depth(goodBins(end)) - EndD > binsize
   disp(' ')
   disp('**Depths are not incremental.')
   disp('**Bad bins were not trimmed.')
   close(f);
   return
else 
   disp('Redefining the "Bin" dimension')
   disp('May take a few minutes...')
	B=resize(B,length(goodBins));
end

%Record some information 
f{'D'}.water_depth(:) = MSL;

disp(['# in Dsurf ' num2str(length(Dsurf))])
disp(['# in ensemble ' num2str(EnsNum)])

%and define a new variables
if ~isempty(Dsurf)
	f{'height'} = ncfloat('ensemble') 
	f{'height'}.long_name = ncchar('height of sea surface from transducer head');
	f{'height'}.units = ncchar('m');
	f{'height'}.FillValue_ = 1.0e35;
	endef(f)

	f{'height'}(1:EnsNum) = Dsurf;
end

close(f)
    
%Done
disp(' ')
disp(['File ' trimFile ' has been modified'])
disp(['## ' num2str(theBinNum-bincnt) ' bins were removed from the top of the water column'])

history(trimFile,thecomment);

if ~running(batch)
	starbeam(trimFile);
end