adcp2ep.m

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

function epDataFile = adcp2ep(adcpFile, epDataFile, ADCPtype, dlgFile)

%function epDataFile = adcp2ep(adcpFile, epDataFile, ADCPtype, dlgFile)
% This function is used to translate RDI ADCP data into variables that are
% in terms of earth coordinates and create an epic compatible data file.
% If the data is in Beam coordinates it will be transformed by runbm2g.m into
% Earth coordinates.  This transformation can be run on workhorse and broad
% band data, but the ADCP type must be specified.
%
% Magnetic Declination
%		If a magnetic declination was provided to the insturment prior to 
%	deployment or to rdi2cdf in post-processing it will be applied at this
%	time for both Earth and Beam coordinat data
%
%Inputs:
%	adcpFile = the ADcp data file in beam coordinates 
%		(Note:if running routines in sequence it should be the trimFile.)
%	epDataFile = the new Epic compatable file that will be created 
%	ADCPtype = WH or BB; will default to WH if not specified
%		WH = workhorse, BB = broad band
%		note: if BB, do not need a dlgFile
%	dlgFile = the dialog file that was created when the ADCP was "deployed"
%
%	Note: If the names of the files are not given, they will be requested.	
%
% Output:
%	epDataFile = same as input


% sub-functions:
%	runbm2g.m
%	ep_time.m
%	gregorian.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.1
%   fixed 30-Apr- 2003 depth calculation for down but no pressure sensor
%   updated 9-Apr-2003  depth calculations to use prssure sensor if available and able to calculate downward looking
%   added note to pressure variable
%   updated so runbm2g history comment will be added to epic file history attributes
%   updated latitude variable so same as longitude variable
%
% updated 19-March-2003 (ALR) added pressure variable
% updated 10-Jan-2003 (ALR) added ';' to stop data stream...
% updated 06-aug-2002 (ALR) Added OCSD to Experiment and LA Shelf ADCP to Description
% updated 01-Jul-2002 (ALR) time conversion won't crash on very short files (line 143)
% updated 29-Apr-2002 (ALR) fixed time stamp error; when moving time stamp to middle of ens, need to add instead of subtract
%       (line 132 corrected)

% version 1.0
% updated 03-Jul-2001 fixed spelling of Buoy (ALR)
% updated 28-Jun-2001 will now work with Matlab 6.0, r.12; problem was with getinfo.m for latitude and
%       longitude values line 428 to 433. (ALR)
% updated 28-Dec-2000 added linefeeds to history/comment attribute (ALR)
% updated 20-Oct-2000 so when computing data in EARTH coordinates the fill values
%		remain the same even after the data is corrected for magnetic declination. ALR
% updated 05-Oct-2000 added additional experiment names and locations.  A. Ramsey
% updated 03-Aug-2000 so creation date will include time.  Fran Hotchkiss
% updated 15-Mar-2000 08:49:03,
%		fixed Pgd and AGC corruption
% updated 02-Feb-2000 09:34:12,
%		fixed temperature weirdness and added rotation for earth data
% updated 16-Dec-1999 12:10:41, put in some clears to help memory
% updated 10-Dec-1999 10:52:22, modified to handel BB
% updated 29-Oct-1999 09:20:26
% updated 19-Oct-1999 13:43:47

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

if running(batch)
   adcpFile = get(batch);
   eval(['adcpFile = ' adcpFile ';'])
   epDataFile = get(batch);
   eval(['epDataFile = ' epDataFile ';'])   
   ADCPtype = get(batch);
   eval(['ADCPtype = ' ADCPtype ';'])
else

if nargin < 1, help(mfilename), adcpFile='';, end
if nargin < 2, epDataFile='';, end
if nargin < 3, ADCPtype = '';, end
if nargin < 4, dlgFile='';, end

if isempty(adcpFile), adcpFile = '*'; end
if isempty(epDataFile), epDataFile = '*';, end
if isempty(ADCPtype), ADCPtype = 'WH';, end
if isempty(dlgFile), dlgFile = '*'; end

% Open ADCP beam file.
if any(adcpFile == '*')
	[theFile, thePath] = uigetfile(adcpFile, 'Select ADCP File:');
	if ~any(theFile), return, end
	if thePath(end) ~= filesep, thePath(end+1) = filesep; end
	adcpFile = [thePath theFile];
end

[path,name,ext,ver]=fileparts(adcpFile);
suggest=[name(1:end-1) '.nc'];
%create ADCP Geographic coordinates File
if any(epDataFile == '*')
	[theFile, thePath] = uiputfile(suggest, 'Save Data in geographic coordinates as:');
	if ~any(theFile), return, end
	if thePath(end) ~= filesep, thePath(end+1) = filesep; end
	epDataFile = [thePath theFile];
end

end %if batch

B = netcdf(adcpFile );
if isempty(B), return, end

%let's deal with the time
ppens=B.pings_per_ensemble(:);
tp=B.time_between_ping_groups(:);
%in order to ping as fast as possible tp may be set to 0,
%but in reality it is approximately 0.288 sec per ping.
if isequal(tp,0);
   tp = 0.288;
end

tens=ppens*tp; 

allTIM=B{'TIM'}(:);
% calculate time step for output file, in seconds (FH 10 May 2000)
timsecs=B{'TIM'}(:)*24*3600;
delta=mean(diff(timsecs))%;
clear timsecs

gtD=gregorian(allTIM);
start_time = datestr(datenum(gtD(1,1),gtD(1,2),gtD(1,3),gtD(1,4),gtD(1,5),gtD(1,6)),0);
stop_time = datestr(datenum(gtD(end,1),gtD(end,2),gtD(end,3),gtD(end,4),gtD(end,5),gtD(end,6)),0);

%Tim is the beginning of the ensemble, make it the middle
tsec=gtD(:,4)*3600 + gtD(:,5)*60 + gtD(:,6);
tmid=tsec+tens/2;
disp(['TIM was corrected by ' num2str(tens/2) ' sec = half the ensemble time']);
disp(' ')
   sec=rem(tmid,60);
   hmt=(tmid-sec)/60;
   minu=rem(hmt,60);
   hr=(hmt-minu)/60;
   Mid_Ens_time=[gtD(:,1) gtD(:,2) gtD(:,3) hr minu sec];
   
disp('Converting "TIM" to "time" and "time2"');
[m,n] = size(Mid_Ens_time);
for ii=1:m;
	alltime(ii,:)=ep_time(Mid_Ens_time(ii,:));
end

Time=alltime(:,1);
Time2=alltime(:,2);

%get some information
theFillValue = fillval(B{'vel1'});
bin = size(B('bin'),1);
ensemble = size(B('ensemble'),1);
xducer_off = B{'D'}.xducer_offset_from_bottom (:);   
wdepth = B{'D'}.water_depth (:); 
serial = B.ADCP_serial_number(:);
coord = B.transform(:);
orientation = B.orientation(:);
%defined these 2 varibales earlier since getting screwed up!
TTX = B{'Tx'};
tempC = TTX(:);


%Depth Calculation: (added 04-Apr-2003)
%   1.  Use pressure sensor if available
%   2.  Use values from Surface.exe calculated using Trimbins.m (only for upward looking)
%   3.  Use user input values calculated using Trimbins.m for upward looking or user input values asked for now for downward looking

switch orientation
case 'UP'
%Is there a pressure variable?
if (B{'Pressure'}) ~=0
    Press = B{'Pressure'}(:);       % pressure at transducer head in pascals
    mnPress = mean(Press);      % average of pressure
    depth_head = mnPress/9806.65;   % depth at transducer head in meters
    center_first_bin = B{'D'}.center_first_bin(:);
    bin_size = B{'D'}.bin_size(:);
    
    depth_head_corrected = depth_head+xducer_off;
    depth = depth_head_corrected - B{'D'}(:);
    dnote = 'Depth values were calculated using the ADCP Pressure Sensor, assuming 9806.65 Pascals per meter.';
    wdepth = depth_head_corrected;
    if ~isempty(B{'height'})
        height = theFillValue*ones(1,ensemble); %Do not need the height variable from surface because the pressure reading is more accurate (ALR 4/4/03)
    end
else
    Press = theFillValue*ones(1,ensemble);
    depth = wdepth - B{'D'}(:);     %wdepth was calculated from surface using trimbins or was a user input value given to trimbins
    %is there a height variable?
    if ~isempty(B{'height'})
        height=B{'height'}(:);
        dnote = 'Depth values were calculated using Surface.exe output';
    else
        height = theFillValue*ones(1,ensemble);
        dnote = 'Depth values were calculated using user input values';
    end
end

case 'DOWN'
   %Is there a pressure variable?
   if (B{'Pressure'}) ~=0
     Press = B{'Pressure'}(:);       % pressure at transducer head in pascals
     mnPress = mean(Press);      % average of pressure
     depth_head = mnPress/9806.65;   % depth at transducer head in meters
     center_first_bin = B{'D'}.center_first_bin(:);
     bin_size = B{'D'}.bin_size(:);
     
    bin1 = depth_head+center_first_bin;
    binEnd = ((bin-1)*bin_size)+bin1; 
    depth = (bin1:bin_size:binEnd)';
    wdepth = depth_head + xducer_off; %water depth is depth at head plus transducer offset.
    dnote = 'Depth values were calculated using the ADCP Pressure Sensor, assuming 9806.65 Pascals per meter.';
    height = theFillValue*ones(1,ensemble);
    
   else
      if running(batch)
          MSL = get(batch);
           eval(['MSL = ' MSL ';'])
           disp(['User input value for depth: ' (num2str(MSL)) ' meters' ])
        else
            disp('User must input the water depth information')
	         Depth_Information.mean_sea_level.value = {0};
   	          Depth_Information.mean_sea_level.units = {'meters'};
	         Depth_Information = uigetinfo(Depth_Information');	
           
   	          %check units
      	       infoD = getinfo(Depth_Information,'mean_sea_level');
                unitD = getinfo(infoD,'units');
              
	       	    if ~isequal(unitD,'meters')
            	    disp('User error!! Depth Information must be in meters')
             	    pause(3)
                  Depth_Information = uigetinfo(Depth_Information)
               end
           
         	  	infoMSL = getinfo(Depth_Information,'mean_sea_level');
                MSL = getinfo(infoMSL,'value');
                 disp(['User input value for depth: ' (num2str(MSL)) ' meters' ])
             end %end if/else batch
             center_first_bin = B{'D'}.center_first_bin(:);
             bin_size = B{'D'}.bin_size(:);
             depth_head = MSL - xducer_off;
             
             bin1 = depth_head+center_first_bin;
             binEnd = ((bin-1)*bin_size)+bin1; 
             depth = (bin1:bin_size:binEnd)';
             dnote = 'Depth values were calculated using user input values';
             Press = theFillValue*ones(1,ensemble);
             height = theFillValue*ones(1,ensemble);
             wdepth = MSL; %water depth equals user input mean sea level
      end %end if pressure/else
end %end switch

disp(['The file has ' num2str(ensemble) ' ensembles and ' num2str(bin) ' bins']);

disp(' ')
disp('Averaging echo intensity')
% Average echo intensity.
i1 = B{'AGC1'};
i1 = autonan(i1,1);
i2 = B{'AGC2'};
i2 = autonan(i2,1);
i3 = B{'AGC3'};
i3 = autonan(i3,1);
i4 = B{'AGC4'};
i4 = autonan(i4,1);
iall(:,:,1) = i1(:,:);
iall(:,:,2) = i2(:,:);
iall(:,:,3) = i3(:,:);
iall(:,:,4) = i4(:,:);
iavg = mean(iall,3);

%free up some memory
clear i1 i2 i3 i4 iall

disp(' ')
disp('Averaging percent good')

switch coord
case 'BEAM'
   
% Average percent good.
p1 = B{'PGd1'};
p1 = autonan(p1,1);
p2 = B{'PGd2'};
p2 = autonan(p2,1);
p3 = B{'PGd3'};
p3 = autonan(p3,1);
p4 = B{'PGd4'};
p4 = autonan(p4,1);
pall(:,:,1) = p1(:,:);
pall(:,:,2) = p2(:,:);
pall(:,:,3) = p3(:,:);
pall(:,:,4) = p4(:,:);
pavg = mean(pall,3);

case 'EARTH'
   p4 = B{'PGd4'};
	p4 = autonan(p4,1);
   pavg = p4(:,:);
end

%free up some memory
clear p1 p2 p3 p4 pall


%check for fill values
	TF = isnan(iavg);
	MTF = max(max(TF));
	if MTF
	   fill_flag = 1;
	end
   
   TF = isnan(pavg);
	MTF = max(max(TF));
	if MTF
   	fill_flag = 1;
	end

%get the velocity data and set the new velocity variable
%dependent on being in Beam or Earth coordinates
switch coord
   
case 'BEAM'
   disp('Data in Beam coordinates is being transformed to Earth')
   %if 0
   if running(batch)
      dlgFile = get(batch);
   	eval(['dlgFile = ' dlgFile ';'])      
   end
   
   cur = runbm2g(adcpFile, ADCPtype, dlgFile);
   %end
   %load cur.mat
   close(B)
   B = netcdf(adcpFile ); %added 09-Apr-2003. Now runbm2g history comment will be added to epic file.
   
   
case 'EARTH'
   disp('Data is in Earth coordinates')
   for ii = 1:4;
      vel{ii} = B{['vel' int2str(ii)]};
   end
   
   %need to adjust for magnetic declination
   u = vel{1}(:);
   v = vel{2}(:);
   magnetic = B{'Hdg'}.heading_bias(:);
   theta = -1*magnetic;
   [ur,vr] = uv_rotate(u,v,theta); 
   %need to reset theFillValue
   urmax = max(max(ur));  vrmax = max(max(vr));
   ur(ur == urmax) = theFillValue;
   vr(vr == vrmax) = theFillValue;
   %fill in rotated velocity values
   [vel{1},vel{2}] = deal(ur,vr);
      
	cur = cell(size(vel));   % Output currents and error.
	p = zeros(4, bin);
   for ii=1:ensemble
		for k = 1:4
  	 		p(k, :) = vel{k}(ii, :);	
			cur{k}(ii, :) = p(k, :);
   	end
      if ~rem(ii,100), 
      disp(sprintf('%d ensembles copied',ii)), 
   	end
   end
end

%scale as cm/sec
q = zeros(ensemble, bin);
for k =1:4;
  	q(:, :) = cur{k}(:, :);	
   q(q == theFillValue) = nan;
   q = q./10;
   q(isnan(q)) = theFillValue;
   cur{k}(:,:) = q;
 end

%Get the mins and maxes for later
for k=1:4
   good_cur=(cur{k}(:,:) ~= theFillValue);
   scur=cur{k}(logical(good_cur));
   minsc(k) = min(min(scur));
   maxsc(k) = max(max(scur));
end
clear p q good_cur scur

% Check for fill values.
   
   for k=1:4;
		TF = isnan(cur{k}(:,:));
		MTF =(max(max(TF)));
		if MTF
   		fill_flag = 1;
		else
   		fill_flag = 0;
	   end
	end