rdi2cdf.m

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

function rdi2cdf(infile, outfile, minens, maxens);

% rdi2cdf.m converts RDI ADCP files to netCDF format
%
% function rdi2cdf(infile, outfile, minens, maxens);
%
% where:
%		infile	= the input file name in RDI broadband data format
%			  this includes the Workhorse family of ADCP's
%		outfile = the netCDF file name
%		minens = ensumble number at which to start converting
%		maxens = ensemble number at which to stop converting
%
% If no outfile is given, the header info for the file is displayed
% Version 3.0	10-Jan-2003


% Written by Marinna Martini
% for the U.S. Geological Survey
% Marine and Coastal Program
% Woods Hole Center, Woods Hole, MA
% http://woodshole.er.usgs.gov/
% Please report bugs to [e-mail]
%
% dependents:
%	rdhead.m	reads RDI header data format
%	rdflead.m	reads RDI fixed leader format
%	rdvlead.m	reads RDI variable leader format
%	mcnote.m	annotate netCDF data
%	mexcdf.mex	netCDF file I/O
%
% netCDF attributes and variable names follow EPIC conventions
%
% information on netCDF may be obtained at
% http://www.unidata.ucar.edu
% it's free!!
%
%
prog_ver = '3.0 10-Jan-2003';
disp(['rdi2cdf.m version ',prog_ver]);

% Updated 03-APr-2003 (ALR) Note numbers for depths
% Updated 27-Feb-2003 (ALR) added ADCP Transducer pressure variable to netcdf file
% Updated 05-Feb-2003 (ALR) adjusted heading_bias attribute so would accept decimals 
% Version 3.0 accomodates downward orientation 10-Jan-2003
% updated 10-Jan-2003 (ALR)  Adjusted NOTE in D so will reflect if bins are relative to seabed for upward facing
%   or relative to transducer for downward facing. And computed downward depths are negative.
% updated 02-Oct-2002 (ALR) Will now read downward facing orientation
% updated 29-Jun-2000 15:18:36
% updated 29-Oct-1999 09:20:26
% version 2.2c  by JMC
% updated 22-Oct-1999 15:00:37
% 7/22/99	Version 2.2b  edited by JMC
%				- request for deployment dates in the netcdf file for future checks
% 7/14/99	Version 2.2  edited by JMC
%				- space added in for some attributes that will be placed in the final
%				processed data file	
% 4/16/99	Version 2.1
%				- I am, in fact, reverting back to mexcdf because
%				  that access code is more than 50% faster.
% 3/19/99	Version 2.0
%				- convert netCDF file access from mexcdf to netcdf
%				  the old code is still contained within
% 8/14/98	Version 1.2b
%				- add the instrument serial number to global attributes
% 11/15/97  Version 1.2a
%				- fix some minor details for USGS data processing needs
% 10/28/97  Version 1.2
%           - fix the incorrect estimation of number of ensembles in the file
% 10/12/97  Version 1.1
%           - fix bin2dec & dec2bin to use MATLAB versions
%           - Add a dialog box to get metadata from user
% 5/2/97    Version 1.0
%           - fix time translation so that year is 1996 instead of '96
%           - include provision for turn of the century, assuming RDI's
%             dates will continue to be two digit year numbers
%
% 10/29/96  Version 0.0

mexcdf('SETOPTS',0);
% no necdf equivalent
if nargin == 0,
	error('I need an input file name');
end
if nargin < 2,	% assume the user only wants header info output
	select=zeros(size(32,1));
	fid = fopen(infile,'r');
	disp(['Header information from RDI adcp file ',infile])
   [nbytes, nt, offsets] = rdhead(fid, 1);
   rdflead(fid,1,select);
	return;
end
tic

% offset indicators
% if RD changes the order of their data types
% be sure to update these types and the other
% info defined below
VELOCITY = 3;
CORRELATION = 4;
INTENSITY = 5;
GOOD = 6;
% the first two read types have no meaning, they are holders
% so that the offset indicator can be used as a universal
% index for code readability
read_types=['int16';'int16';'int16';'uchar';'uchar';'uchar'];
% likewise, the first two rec ids are for the fixed and
% variable leader records
rec_ids=[0 128 256 512 768 1024];
longnames = str2mat(' ',' ','velocity','correlation','intensity','percent good');
VARNAMES=['   ';'   ';'vel';'cor';'AGC';'PGd'];

%
% ------ Open the files & write global attributes -----
%
% create the file
cdf = mexcdf('CREATE',outfile,'NC_CLOBBER');
if cdf == -1,
	error(['Problem opening netCDF output file: ',outfile]);
end
%cdf = netcdf(outfile,'clobber');	%clobber = create or overwrite a new file
%if isempty(cdf), error(['Problem opening netCDF output file: ',outfile]); end

if running(batch)
   Mooring_number = get(batch);
   eval(['Mooring_number = ' Mooring_number ';'])
   Deployment_date = get(batch);
   eval(['Deployment_date = ' Deployment_date ';'])
   Recovery_date = get(batch);
   eval(['Recovery_date = ' Recovery_date ';'])
else
	% get the deployment dates
	prompt  = {'Enter 4-digit mooring number:',...
      'Enter the deployment date (dd-mmm-yyyy):', 'Enter the recovery date (dd-mmm-yyyy):'};
	def     = {'0','01-jan-1999','01-jan-1999'};
	title   = 'Input mooring data for ADCP';
	lineNo  = 1;
	dlgresult  = inputdlg(prompt,title,lineNo,def);
	Mooring_number = dlgresult{1};
	Deployment_date = dlgresult{2};
	Recovery_date = dlgresult{3};
end %batch

% enter the global attributes   
mexcdf('ATTPUT',cdf,'GLOBAL','CREATION_DATE','CHAR',11,date);
%cdf.CREATION_DATE = date;

mexcdf('ATTPUT',cdf,'GLOBAL','Mooring_number','CHAR',4,Mooring_number);
mexcdf('ATTPUT',cdf,'GLOBAL','Deployment_date','CHAR',11,Deployment_date);
mexcdf('ATTPUT',cdf,'GLOBAL','Recovery_date','CHAR',11,Recovery_date);

mexcdf('ATTPUT',cdf,'GLOBAL','INST_TYPE','CHAR',19,'RD Instruments ADCP');
%cdf.INST_TYPE = 'RD Instruments Broadband ADCP'
junk = ['Converted to netCDF via MATLAB by rdi2cdf.m ', prog_ver];
mexcdf('ATTPUT',cdf,'GLOBAL','History','CHAR',length(junk),junk);
%cdf.PROG_CMNT1 = junk;

fid = fopen(infile,'r');
%function [nb, nt, off] = rdhead(fid, verbose);
%	Read the header data from a raw ADCP
%	data file opened for binary reading.
%	fid = file handle returned by fopen
%	nb = number of bytes in the ensemble
%	nt = number of data types
%	off = offset to the data for each type
%	Set verbose = 1 in rdhead.m for a text output.
disp(['Header information from adcp file ',infile])
[nbytes, nt, offsets] = rdhead(fid, 1)
notenum = 1;	% for general notes to file

% this is essentially the code from rdflead.m
% adapted for writing attributes to the netcdf file
%
% ------ Process Fixed Leader Data --------
%
%	Read the fixed leader data from a raw ADCP
%	data file opened for binary reading 
%	Returns the contents of the fixed leader
%	as 31 elements of the vector 'data' or an
%	empty matrix if the fixed leader ID is not
%	identified (error condition)
%	Set verbose=1 for a text output.
%	The names and meanings of these attibutes are taken
%	directly from Appendix B in the RDI Broadband 
%	Phase III technical manual.

if exist('verbose') ~= 1,
	verbose = 1;	% verbose must be = 1 to get all header info
end
NFIELDS = 32;
data=zeros(1,NFIELDS);
fld=1;  

if exist('select') ~= 1,
	select = [];
end

if running(batch)
   ADCP_serial_number = get(batch);
   eval(['ADCP_serial_number = ' ADCP_serial_number ';'])
   xducer_offset = get(batch);
   eval(['xducer_offset = ' xducer_offset ';'])
   pred_accuracy = get(batch);
   eval(['pred_accuracy = ' pred_accuracy ';'])
   slow_by = get(batch);
   eval(['slow_by = ' slow_by ';'])
   magnetic = get(batch);
   eval(['magnetic = ' magnetic ';'])

else
	% get the distance of the ADCP head from the bottom
	% get the instrument serial number
	prompt  = {'Enter the ADCP''s serial number:',...
	'Enter the distance between the ADCP transducers and the sea bed in meters:',...
   'Enter the predicted accuracy given by PLAN in cm/s:',...
	'Enter the amount of time the ADCP clock was slow by in seconds:',...
	'Enter the magnetic variation at the mooring location in degrees'};
      
	def     = {'0','0','0','0','0'};
	title   = 'Input metadata from the mooring log';
	lineNo  = 1;
	dlgresult  = inputdlg(prompt,title,lineNo,def);
	ADCP_serial_number = str2num(dlgresult{1});	
	xducer_offset = str2num(dlgresult{2});
   pred_accuracy = str2num(dlgresult{3});
   slow_by = str2num(dlgresult{4});
	magnetic = str2num(dlgresult{5});

	%let's check to make sure these values make sense
	if xducer_offset > 10 | ADCP_serial_number < 100; 
   	prompt  = {'RE-Enter the ADCP''s serial number:',...
         'RE-Enter the distance between the ADCP and sea bed in METERS:'};
		def     = {dlgresult{1},dlgresult{2}};
		title   = 'Check metadata, INVALID values';
		lineNo  = 1;
		dlgresult  = inputdlg(prompt,title,lineNo,def);
		ADCP_serial_number = str2num(dlgresult{1});
		xducer_offset = str2num(dlgresult{2});
   end
   
end %batch

% make sure we're looking at the beginning of
% the fixed leader record by testing for it's ID
data(fld)=fread(fid,1,'ushort');
if(data(fld)~=0),
	disp('Fixed Leader ID not found');
	data=[];
	return;
end
fld=fld+1;
% version number of CPU firmware
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
% revision number of CPU firmware
data(fld)=fread(fid,1,'uchar');
if verbose, disp(sprintf('CPU Version %d.%d',data(fld-1),data(fld))); end;
mexcdf('ATTPUT',cdf,'GLOBAL','firmware_version','FLOAT',1,data(fld-1)+data(fld)/100);
%cdf.firmware_version = ncfloat(data(fld-1)+data(fld)/100);
fld=fld+1;
% configuration, uninterpreted
data(fld)=fread(fid,1,'uchar');
if verbose, 
	disp(sprintf('Hardware Configuration for LSB %d',data(fld))); 
   b=zeros(1,8); %dec2bin(data(fld));
   b(9-length(dec2bin(data(fld))):8)=dec2bin(data(fld));
   b=char(b);
   freqs=[75 150 300 600 1200 2400];
	junk=bin2dec(b(6:8));
	disp(sprintf('	System Frequency = %d kHz',freqs(junk+1))); 
   mexcdf('ATTPUT',cdf,'GLOBAL','frequency','LONG',1,freqs(junk+1));
   %cdf.frequency = nclong(freqs(junk+1));
	if b(5) == '0', 
		disp('	Concave Beam'); 
      mexcdf('ATTPUT',cdf,'GLOBAL','beam_pattern','CHAR',7,'concave');
      %cdf.beam_pattern = 'concave';
	end		
	if b(5) == '1', 
		disp('	Convex Beam'); 
      mexcdf('ATTPUT',cdf,'GLOBAL','beam_pattern','CHAR',6,'convex');
      %cdf.beam_pattern = 'convex';
	end		
	junk=bin2dec(b(3:4));
	disp(sprintf('Sensor Configuration #%d',junk+1)); 
   mexcdf('ATTPUT',cdf,'GLOBAL','sensor_configuration','LONG',1,junk+1);
   %cdf.sensor_configuration = nclong(junk+1);
	if b(2) == '0', disp('	Transducer head not attached'); end		
	if b(2) == '1', disp('	Transducer head attached'); end		
   mexcdf('ATTPUT',cdf,'GLOBAL','transducer_attached','LONG',1,b(2));
   %cdf.transducer_attached = nclong(b(2));
	if b(1) ~= '1', 
		disp('	Downward facing beam orientation'); 
      mexcdf('ATTPUT',cdf,'GLOBAL','orientation','CHAR',4,'DOWN');
      %cdf.orientation = 'DOWN';
      orientation='down'; %added alr 01/10/03
	end		
	if b(1) == '1', 
		disp('	Upward facing beam orientation'); 
		mexcdf('ATTPUT',cdf,'GLOBAL','orientation','CHAR',2,'UP');
      %cdf.orientation = 'UP';
       orientation='up'; %added alr 01/10/03
	end		
end;
fld=fld+1;
data(fld)=fread(fid,1,'uchar');
if verbose, 
	disp(sprintf('Hardware Configuration MSB %d',data(fld))); 
	%b=dec2bin(data(fld));
   b=zeros(1,8); %dec2bin(data(fld));
   b(9-length(dec2bin(data(fld))):8)=dec2bin(data(fld));
   b=char(b);
	angles = [15 20 30 0];
	junk=bin2dec(b(7:8));
	disp(sprintf('	Beam angle = %d degrees',angles(junk+1)));
	% note that is beam angle = 0, it is some other than the above
   mexcdf('ATTPUT',cdf,'GLOBAL','beam_angle','LONG',1,angles(junk+1));
   %cdf.beam_angle = nclong(angles(junk+1));
	junk=bin2dec(b(1:4));
	if junk == 4, 
		disp('	4-beam janus configuration'); 
      mexcdf('ATTPUT',cdf,'GLOBAL','janus','CHAR',6,'4 Beam');
      %cdf.janus = '4 Beam';
	end
	if junk == 5, 
		disp('	5-beam janus configuration, 3 demodulators'); 
      mexcdf('ATTPUT',cdf,'GLOBAL','janus','CHAR',15,'5 Beam, 3 demod');
      %cdf.janus = '5 Beam, 3 demod';
	end
	if junk == 15, 
		disp('	4-beam janus configuration, 2 demodulators'); 
      mexcdf('ATTPUT',cdf,'GLOBAL','janus','CHAR',15,'4 Beam, 2 demod');
      %cdf.janus = '4 Beam, 2 demod';
	end
end;
fld=fld+1;
% real (0) or simulated (1) data flag
data(fld)=fread(fid,1,'uchar');
if data(fld),
	if verbose, disp('The data is simulated'); end
else
	if verbose, disp('The data is real'); end
end
mexcdf('ATTPUT',cdf,'GLOBAL','simulated_data','LONG',1,data(fld));
%cdf.simulated_data = nclong(data(fld));
fld=fld+1;
% undefined
data(fld)=fread(fid,1,'uchar');	fld=fld+1;
% number of beams
data(fld)=fread(fid,1,'uchar');	
if verbose, disp(sprintf('Number of beams used to calculate velocity data:  %d',data(fld))); end
mexcdf('ATTPUT',cdf,'GLOBAL','beams_in_velocity_calculation','LONG',1,data(fld));
%cdf.beams_in_velocity_calculation = nclong(data(fld));
fld=fld+1;
nbeams = 4;	%data(fld)
% number of depth cells
data(fld)=fread(fid,1,'uchar');
if verbose, disp(sprintf('Number of depth cells %d',data(fld))); end;
nbins = data(fld);
fld=fld+1;
% pings per ensemble
data(fld)=fread(fid,1,'ushort');
if verbose, disp(sprintf('Pings per ensemble %d',data(fld))); end;
mexcdf('ATTPUT',cdf,'GLOBAL','pings_per_ensemble','LONG',1,data(fld));
%cdf.pings_per_ensemble = nclong(data(fld));
npings = data(fld);
fld=fld+1;
% depth cell length in cm