function values = nc_varget(ncfile, varname, varargin )
% NC_VARGET:  Retrieve data from a NetCDF variable.
%
% NC_VARGET tries to be intelligent about retrieving the data.
% Since most general matlab operations are done in double precision,
% retrieved numeric data will be cast to double precision, while 
% character data remains just character data.  
%
% Singleton dimensions are removed from the output data.  
%
% A '_FillValue' attribute is honored by flagging those datums as NaN.
% A 'missing_value' attribute is honored by flagging those datums as 
% NaN.  The exception to this is for NC_CHAR variables, as mixing 
% character data and NaN doesn't really seem to work in matlab.
%
% If the named NetCDF variable has valid scale_factor and add_offset 
% attributes, then the data is scaled accordingly.  
%
% USAGE:  values = nc_varget(ncfile,varname);
%         
%     Retrieve an entire dataset.
%
%
% USAGE:  values = nc_varget(ncfile,varname,start,count);
%         
%     Retrieve a contiguous section of a dataset, but not necessarily
%     all of it.
%
%
% USAGE:  values = nc_varget(ncfile,varname,start,count,stride);
%         
%     Retrieve a non-contiguous section of a dataset.
%
%
% PARAMETERS:
% Input:
%    ncfile:  
%       name of netcdf file.  This can also be a DODS url.
%    varname:  
%       name of variable whose data is requested
%    start:  
%       Optional, starting index for the retrieval.  Remember that 
%       NetCDF has zero-based indexing, not one-based indexing.
%    count:  
%       Optional, length of the retrieval along each dimension.  If the
%       count for any dimension is -1, then everything is retrieved up
%       to the end of that dimension.  Requires that the "start" 
%       argument be given as well.
%    stride:  
%       Optional.  If not provided, then contiguous elements are used.
%       For example, if STRIDE were [2 2], then every other element
%       is chosen.  Requires that both the "start" and "count" arguments
%       be given as well.
%
% Output:
%    values:  the requested data
%
%
% EXAMPLE:
% #1.  In this case, the variable in question has rank 2, and has size 
%      500x700.  We want to retrieve starting at row 300, column 250.
%      We want 100 contiguous rows, 200 contiguous columns.
% 
%      vardata = nc_varget ( file, variable_name, [300 250], [100 200] );
%
% In case of an error, an exception is thrown.
%
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% $Id: nc_varget.m 1976 2006-11-01 19:55:24Z johnevans007 $
% $LastChangedDate: 2006-11-01 14:55:24 -0500 (Wed, 01 Nov 2006) $
% $LastChangedRevision: 1976 $
% $LastChangedBy: johnevans007 $
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





if (nargin < 2)
	error ( '%s:  too few input arguments.\n', mfilename );
end
if (nargin > 5)
	error ( '%s:  too many input arguments.\n', mfilename );
end



%
% Error checking on the inputs.
if ~ischar(ncfile)
	error ( '%s:  the filename must be character.\n', mfilename );
end
if ~ischar(varname)
	error ( '%s:  the variable name must be character.\n', mfilename );
end

if nargin == 3
	error ( '%s:  there cannot be three input arguments.\n', mfilename );
elseif nargin == 4
	if ~isnumeric ( varargin{1} )
		error ( '%s:  the ''start'' argument must be numeric.\n', mfilename );
	end
	if ~isnumeric ( varargin{2} )
		error ( '%s:  the ''count'' argument must be numeric.\n', mfilename );
	end
elseif nargin == 5
	if ~isnumeric ( varargin{1} )
		error ( '%s:  the ''start'' argument must be numeric.\n', mfilename );
	end
	if ~isnumeric ( varargin{2} )
		error ( '%s:  the ''count'' argument must be numeric.\n', mfilename );
	end
	if ~isnumeric ( varargin{3} )
		error ( '%s:  the ''stride'' argument must be numeric.\n', mfilename );
	end
end





%
% One optional output argument.
if ( nargout > 1 )
    error ( '%s:  at most 1 output argument is allowed.\n', mfilename );
end



%
% Do we use java instead of mexnc?
if strcmp ( nc_method, 'java' )
	values = nc_varget_java(ncfile, varname, varargin{:} );
	return
end





% One of these is set when we figure out what arguments we
% were given.
USE_GET_VAR = 0;   % get the entire variable?
USE_GET_VARA = 0;  % get a contiguous part of the variable?
USE_GET_VARS = 0;  % get a variable with some stride%
USE_GET_VAR1 = 0;  % get just a single value?

% Don't assume we are given coordinates unless we know otherwise.
start = [];
count = [];
stride = [];

%






%
% Set the 
switch nargin
case 2
	USE_GET_VAR = 1;


case 4
	%
	% 3rd and 4th must be numeric.  start and count.
	USE_GET_VARA = 1;
	start = varargin{1};
	count = varargin{2};

case 5
	%
	% 3rd and 4th must be numeric.  start and count.
	USE_GET_VARS = 1;
	start = varargin{1};
	count = varargin{2};
	stride = varargin{3};

end



%
% Open netCDF file
%
[ncid,status]=mexnc('open',ncfile,'NOWRITE');
if status ~= 0
	ncerr = mexnc('strerror', status);
	mexnc('close',ncid);
	strformat = '%s:  mexnc:open failed on file ''%s'', error msg ''%s''.\n';
	error ( strformat, mfilename, ncfile, ncerr ); 
end


%
% Get variable id
%
[varid, status]=mexnc('inq_varid',ncid,varname);
if status ~= 0
	ncerr = mexnc('strerror', status);
	mexnc('close',ncid);
	strformat = '%s:  mexnc:inq_varid failed on variable ''%s'', file ''%s'', error msg ''%s''.\n';
	error ( strformat, mfilename, varname, ncfile, ncerr ); 
end

[dud,var_type,nvdims,var_dim,dud,status]=mexnc('inq_var',ncid,varid);
if status ~= 0
	ncerr_msg = mexnc ( 'strerror', status );
	mexnc('close',ncid);
	error ( '%s:  INQ_VAR failed on %s, file %s (%s)\n', mfilename, varname, ncerr_msg );
end



%
% Check that the start, count, stride parameters have appropriate lengths.
% Otherwise we get confusing error messages later on.
if ( USE_GET_VAR1 | USE_GET_VARA | USE_GET_VARS )
	if (length(start) ~= nvdims)
		mexnc('close',ncid);
		error ( '%s:  length of the start argument (%d) does not equal the number of dimensions (%d) for %s\n', mfilename, length(start), nvdims, varname );
	end
end
if ( USE_GET_VARA | USE_GET_VARS )
	if (length(count) ~= nvdims)
		mexnc('close',ncid);
		error ( '%s:  length of the count argument (%d) does not equal the number of dimensions (%d) for %s\n', mfilename, length(count), nvdims, varname );
	end
end
if ( USE_GET_VARS )
	if (length(stride) ~= nvdims)
		mexnc('close',ncid);
		error ( '%s:  length of the stride argument (%d) does not equal the number of dimensions (%d) for %s\n', mfilename, length(stride), nvdims, varname );
	end
end


%
% If a singleton, then use GET_VAR1.  This is only because some 
% opendap-enabled mexnc clients have trouble using GET_VAR on 
% singletons.  It is annoying to have to do this, but it works just as
% well.
if nvdims == 0

	%
	% Must be a scalar variable.  
	the_var_size = 1;
	USE_GET_VAR1 = 1;
	USE_GET_VAR = 0;
	start = 0;
	count = 1;

%
% otherwise we need to figure out how big the variable is.
else
	for n=1:nvdims,
	    dimid=var_dim(n);
	    [dud,dim_size,status]=mexnc('inq_dim', ncid, dimid);
	    if ( status ~= 0 )
		ncerr_msg = mexnc ( 'strerror', status );
		mexnc('close',ncid);
		error ( '%s:  INQ_DIM failed on dimid %d, file %s (%s)\n', mfilename, dimid, ncerr_msg );
	    end
	    the_var_size(n)=dim_size;
	end
end


%
% If the user had set non-positive numbers in "count", then
% we replace them with what we need to get the
% rest of the variable.
if USE_GET_VARA | USE_GET_VARS

	
	for j=1:nvdims
		if ( count(j) <= 0 )
			count(j) = the_var_size(j) - start(j);
		end
	end

end



%
% Is the requested size negative?
ind = find(count<0);
if ~isempty(ind)
	mexnc('close',ncid);
	hyperslab_start = sprintf ( '[%d]', start );
	hyperslab_request = sprintf ( '[%d]', count );
	hyperslab_extent = sprintf ( '[%d]', the_var_size );
	strformat = '%s:  requested ''%s'' data hyperslab extent %s + %s is negative, restate to fit existing data bounds %s.\n';
	error ( strformat, mfilename, varname, hyperslab_start, hyperslab_request, hyperslab_extent );
end


if isempty(stride)
	stride = ones(size(count));
end


%
% Make sure that all three of start, count, and stride have
% the same length
if ( length(start) ~= length(count) )
	error ( '%s:  length of start argument must be the same as the count argument.\n', mfilename );
else
	if ( length(start) ~= length(stride) )
		error ( '%s:  length of start argument must be the same as the stride argument.\n', mfilename );
	end
end

%
% Does the hyperslab exceed the dataset size?
if USE_GET_VARA | USE_GET_VARS
	last_datum_location_retrieved = start + (count-1).*stride;
	ind = find (last_datum_location_retrieved >= the_var_size );
	if ~isempty(ind)
		mexnc('close',ncid);
		hyperslab_start = sprintf ( '[%d]', start );
		hyperslab_request = sprintf ( '[%d]', count );
		hyperslab_stride = sprintf ( '[%d]', stride );
		hyperslab_extent = sprintf ( '[%d]', the_var_size );
		format = '%s:  requested ''%s'' data hyperslab start %s, count %s, stride %s exceeds existing data bounds %s.\n';
		error ( format, mfilename, varname, hyperslab_start, hyperslab_request, hyperslab_stride, hyperslab_extent );
	end
end




%
% Get slab, unless the count is zero.
%
% There is a bug in Mexcdf.  If the count of the
% first dimension is zero, it core dumps.  Something to figure
% out later.
if prod(count) == 0
	mexnc('close',ncid);
	return
end


nc_datatype = nc_datatype_string(var_type);


%
% want slowest varying dimension to be # of rows
%if nvdims == 0
if USE_GET_VAR1


	switch ( nc_datatype )
		case 'NC_CHAR'
			[values, status] = mexnc ( 'get_var1_text', ncid, varid, 0 );
		case { 'NC_DOUBLE', 'NC_FLOAT', 'NC_INT', 'NC_SHORT', 'NC_BYTE' }
			[values, status] = mexnc ( 'get_var1_double', ncid, varid, 0 );
		otherwise
			mexnc('close',ncid);
			ncerr_msg = mexnc ( 'strerror', status );
			error ( '%s:  unhandled datatype in GET_VAR1_XXX block, varname %s, file %s (%s)\n', mfilename, varname, ncfile, ncerr_msg );
	end



%
% Get everything.
%elseif nargin == 2
elseif USE_GET_VAR

	switch ( nc_datatype )
		case 'NC_CHAR'
			[values, status] = mexnc ( 'get_var_text', ncid, varid);
		case { 'NC_DOUBLE', 'NC_FLOAT', 'NC_INT', 'NC_SHORT', 'NC_BYTE' }
			[values, status] = mexnc ( 'get_var_double', ncid, varid);
		otherwise
			mexnc('close',ncid);
			ncerr_msg = mexnc ( 'strerror', status );
			error ( '%s:  unhandled datatype in GET_VAR_XXX block, varname %s, file %s (%s)\n', mfilename, varname, ncfile, ncerr_msg );

	end

%elseif nargin == 5
elseif USE_GET_VARS

	switch ( nc_datatype )
		case 'NC_CHAR'
			[values, status] = mexnc ( 'get_vars_text', ncid, varid, start, count, stride );
		case { 'NC_DOUBLE', 'NC_FLOAT', 'NC_INT', 'NC_SHORT', 'NC_BYTE' }
			[values, status] = mexnc ( 'get_vars_double', ncid, varid, start, count, stride );
		otherwise
			mexnc('close',ncid);
			ncerr_msg = mexnc ( 'strerror', status );
			error ( '%s:  unhandled datatype in GET_VARS_XXX block, varname %s, file %s (%s)\n', mfilename, varname, ncfile, ncerr_msg );

	end


elseif USE_GET_VARA

	switch ( nc_datatype )
		case 'NC_CHAR'
			[values, status] = mexnc ( 'get_vara_text', ncid, varid, start, count);
		case { 'NC_DOUBLE', 'NC_FLOAT', 'NC_INT', 'NC_SHORT', 'NC_BYTE' }
			[values, status] = mexnc ( 'get_vara_double', ncid, varid, start, count);
		otherwise
			mexnc('close',ncid);
			ncerr_msg = mexnc ( 'strerror', status );
			error ( '%s:  unhandled datatype in GET_VARA_XXX block, varname %s, file %s (%s)\n', mfilename, varname, ncfile, ncerr_msg );

	end

else

	error ( '%s:  We could not figure out which one of the ''get_var?_xxx'' routines to use.\n', mfilename );

end



if ( status < 0 )
	mexnc('close',ncid);
	error ( '%s:  varget operation failed on %s, %s.\n', mfilename, varname, ncfile );
end

%
% If it's a 1D vector, make it a column vector.
if length(the_var_size) == 1
	values = values(:);
	%values = values';
elseif ( length(the_var_size) > 1 )
	pv = fliplr ( 1:length(the_var_size) );
	values = permute(values,pv);
end                                                                                   


%
% Handle the fill value, if any.  Change those values into NaN.
[dud, dud, status] = mexnc('INQ_ATT', ncid, varid, '_FillValue' );
if ( status == 0 )

	switch ( nc_datatype )
	case 'NC_CHAR'
		%
		% For now, do nothing.  Does a fill value even make sense with char data?
		% If it does, please tell me so.
		%

	case { 'NC_DOUBLE', 'NC_FLOAT', 'NC_INT', 'NC_SHORT', 'NC_BYTE' }
		[fill_value, status] = mexnc ( 'get_att_double', ncid, varid, '_FillValue' );
		fill_value_indices = find(values==fill_value);
		values(fill_value_indices) = NaN;

	otherwise
		mexnc('close',ncid);
		error ( '%s:  unhandled datatype %s\n', mfilename, nc_datatype_string(var_type) );
	end

	if ( status ~= 0 )
		mexnc('close',ncid);
		ncerrmsg = mexnc ( 'strerror', status );
		error ( '%s:  GET_ATT_XXX failed on _FillValue, var %s, file %s (%s).\n', mfilename, varname, ncfile, ncerrmsg );
	end



end









%
% Handle the missing value, if any.  Change those values into NaN.
[dud, dud, status] = mexnc('INQ_ATT', ncid, varid, 'missing_value' );
if ( status == 0 )

	switch ( nc_datatype )
	case 'NC_CHAR'
		%
		% For now, do nothing.  Does a fill value even make sense with char data?
		% If it does, please tell me so.
		%

	case { 'NC_DOUBLE', 'NC_FLOAT', 'NC_INT', 'NC_SHORT', 'NC_BYTE' }
		[missing_value, status] = mexnc ( 'get_att_double', ncid, varid, 'missing_value' );
		missing_value_indices = find(values==missing_value);
		values(missing_value_indices) = NaN;

	otherwise
		mexnc('close',ncid);
		error ( '%s:  unhandled datatype %s\n', mfilename, nc_datatype_string(var_type) );
	end

	if ( status ~= 0 )
		mexnc('close',ncid);
		ncerrmsg = mexnc ( 'strerror', status );
		error ( '%s:  GET_ATT_XXX failed on missing_value, var %s, file %s (%s).\n', mfilename, varname, ncfile, ncerrmsg );
	end



end









%
% Handle the scale factor and add_offsets.   It's certainly possible to 
% have a scale_factor without an add_offset.  Some envisat data does 
% this.
have_scale = 0;
have_addoffset = 0;
[dud, dud, status] = mexnc('INQ_ATT', ncid, varid, 'scale_factor' );
if ( status == 0 )
	have_scale = 1;
end
[dud, dud, status] = mexnc('INQ_ATT', ncid, varid, 'add_offset' );
if ( status == 0 )
	have_addoffset = 1;
end



if have_scale 
	switch ( nc_datatype )
	case 'NC_CHAR'
		%
		% do nothing.

	case { 'NC_DOUBLE', 'NC_FLOAT', 'NC_INT', 'NC_SHORT', 'NC_BYTE' }
		[scale_factor, status] = mexnc ( 'get_att_double', ncid, varid, 'scale_factor' );
		if ( status ~= 0 )
			mexnc('close',ncid);
			error ( '%s:  ATT_GET_DOUBLE failed on scale_factor.\n', mfilename );
		end

		if have_addoffset
			[add_offset, status] = mexnc ( 'get_att_double', ncid, varid, 'add_offset' );
			if ( status ~= 0 )
				mexnc('close',ncid);
				error ( '%s:  ATT_GET_DOUBLE failed on add_offset.\n', mfilename );
			end
		else
			add_offset = 0.0;
		end
		
		values = values * scale_factor + add_offset;
	
	otherwise
		mexnc('close',ncid);
		error ( '%s:  unhandled datatype %s\n', mfilename, nc_datatype_string(var_type) );
	end
end





%
% remove any singleton dimensions.
values = squeeze ( values );



mexnc('close',ncid);


return

















