function [Z,refvec] = etopo(varargin)
%ETOPO Read global 5-min or 2-min digital terrain data.
%
%  [Z, REFVEC] = ETOPO reads the ETOPO data for the entire world from the
%  ETOPO data in the current directory. The current directory is searched
%  first for ETOPO2 binary data, followed by ETOPO5 binary data, followed
%  by ETOPO5 ASCII data from the file names etopo5.northern.bat and
%  etopo5.southern.bat. Once a match is found the data is read. The data
%  grid, Z, is returned as an array of elevations. Data values are in whole
%  meters, representing the elevation of the center of each cell.  REFVEC
%  is the associated referencing vector.
%
%  [Z, REFVEC] = ETOPO(SAMPLEFACTOR) reads the data for the entire world,
%  downsampling the data by SAMPLEFACTOR.  SAMPLEFACTOR is a scalar
%  integer, which when equal to 1 gives the data at its full resolution
%  (1080 by 4320 values for ETOPO5 data and 5400 by 10800 values for ETOPO2
%  data).  When SAMPLEFACTOR is an integer n greater than one, every n-th
%  point is returned.  SAMPLEFACTOR must divide evenly into the number of
%  rows and columns of the data file.  If SAMPLEFACTOR is omitted or empty,
%  it defaults to 1. 
%
%  [Z, REFVEC] = ETOPO(SAMPLEFACTOR, LATLIM, LONLIM) reads the data for the
%  part of the world within the specified latitude and longitude limits.
%  The limits of the desired data are specified as two element vectors of
%  latitude, LATLIM, and longitude, LONLIM, in degrees. The elements of
%  LATLIM and LONLIM must be in ascending order.  LONLIM must be specified
%  in the range [0 360] for ETOPO5 data and [-180 180] for ETOPO2 data. If
%  LATLIM is empty the latitude limits are [-90 90]. If LONLIM is empty,
%  the longitude limits are determined by the file type.
%
%  [Z, REFVEC] = ETOPO(DIRECTORY, ...) allows the path for the ETOPO data
%  file to be specified by DIRECTORY rather than the current directory. 
%
%  [Z, REFVEC] = ETOPO(FILE, ...) reads the ETOPO data from FILE, where
%  FILE is a string or a cell array of strings containing the name or names
%  of the ETOPO data files. 
%
%  For details on locating ETOPO data for download over the Internet, see
%  the following documentation at the MathWorks web site: 
%
%  <a href="matlab: 
%  web('http://www.mathworks.com/support/tech-notes/2100/2101.html#etopo5') 
%  ">http://www.mathworks.com/support/tech-notes/2100/2101.html</a>
%
%
%  Example
%  -------
%  % Read and display the ETOPO5 data from the directory 'etopo5' 
%  % downsampled by a factor of 10.
%  [Z, refvec] = etopo('etopo5',10);
%  axesm merc
%  geoshow(Z, refvec, 'DisplayType', 'surface');
%  colormap(demcmap(Z));
%
%  % Read the ETOPO2 binary data downsampled by a factor of 10.
%  [Z, refvec] = etopo('ETOPO2.dos.bin', 10);
%
%  See also: GTOPO30, TBASE, USGSDEM.

%  Copyright 1996-2004 The MathWorks, Inc.
%  $Revision: 1.1.6.1 $  $Date: 2004/12/18 07:46:26 $

% Check number of inputs
checknargin(0,4,nargin,mfilename);

% Parse input arguments
[file, samplefactor, latlim, lonlim, readFcn] = parseInputs(varargin{:});

% Read the data
[Z,refvec] = readFcn(file, samplefactor, latlim, lonlim);

%----------------------------------------------------------------------------
function [file, samplefactor, latlim, lonlim, readFcn] = parseInputs(varargin)

% Parse varargin{1}
if nargin == 0 || (~iscell(varargin{1}) && ~ischar(varargin{1}))
   % No inputs, set to pwd
   dataDir = pwd;
elseif iscell(varargin{1})
   % sort data files so that 'northern' file is always first
   file = sort(varargin{1});
   varargin(1) = [];
   dataDir = '';
elseif  exist(varargin{1},'dir')
   dataDir = varargin{1};
   varargin(1) = [];
else
   dataDir = '';
   file = {varargin{1}};
   varargin(1) = [];
end

% Check for file or directory input
if ~isempty(dataDir)
   file = getFileFromDir(dataDir);
end

% Verify files exist
for i=1:numel(file)
   if ~exist(file{i},'file')
      eid = sprintf('%s:%s:fileNotFound', getcomp, mfilename);
      msg = sprintf('ETOPO data file "%s" does not exist.',file{1} );
      error(eid, '%s', msg);
   end
end

% Assign default read handle
%  and parse for ETOPO5 or ETOPO2 files
switch numel(file)
   case 0
      eid = sprintf('%s:%s:tooFewFiles', getcomp, mfilename);
      error(eid, 'No input file is specified.')
   case 1
      [null, null, ext] = fileparts(file{1});
      if isequal(lower(ext),'.bil')
         % new_etopo5.bil
         readFcn = @etopo5BinaryRead;
      elseif ~isequal(lower(ext),'.bat')
         % Check for ETOPO2 byte size
         list = dir(file{1});
         if ~isempty(list) && list(1).bytes == 116640000
            readFcn = @etopo2BinaryRead;
         else
            readFcn = @etopo5BinaryRead;
         end
      else
         readFcn = @etopo5AsciiRead;
      end
   case 2
      % Cell array input, expecting etopo5*bat files
      [null, null, ext1] = fileparts(file{1});
      [null, null, ext2] = fileparts(file{2});
      if ~isequal(ext1, ext2) || ~isequal(lower(ext1),'.bat')
         eid = sprintf('%s:%s:invalidExtensions', getcomp, mfilename);
         msg = sprintf('Filenames have an invalid extension. Expecting ".bat"');
         error(eid, '%s', msg);
      end
      readFcn = @etopo5AsciiRead;
   otherwise
      eid = sprintf('%s:%s:tooManyFiles', getcomp, mfilename);
      error(eid, 'Too many input files.')
end

% Parse samplefactor, latlim, lonlim
switch numel(varargin)
   case 0
      samplefactor = 1;
      latlim = [-90 90];
      lonlim = [0 360];
   case 1
      samplefactor = varargin{1};
      latlim = [-90 90];
      lonlim = [0 360];
   case 2
      [samplefactor,latlim] = deal(varargin{1:2});
      if isempty(latlim); latlim = [-90 90]; end
      lonlim = [0 360];
   case 3
      [samplefactor,latlim,lonlim] = deal(varargin{1:3});
      if isempty(latlim); latlim = [-90 90]; end
      if isempty(lonlim); lonlim = [0 360]; end
end

% Verify samplefactor, latlim, lonlim
if isempty(samplefactor)
   samplefactor = 1;
elseif ~isscalar(samplefactor)
   eid = sprintf('%s:%s:notScalar', getcomp, mfilename);
   error(eid,'SAMPLEFACTOR must be a scalar');
end

checkinput(latlim,{'double'},{'real','vector','finite'},mfilename, ...
   'LATLIM',3);
checkinput(lonlim,{'double'},{'real','vector','finite'},mfilename, ...
   'LONLIM',4);

if numel(latlim) ~= 2
   eid = sprintf('%s:%s:invalidLatLim', getcomp, mfilename);
   error(eid,'latlim must be a two element vector in units of degrees');
end

if numel(lonlim) ~= 2
   eid = sprintf('%s:%s:invalidLonLim', getcomp, mfilename);
   error(eid,'lonlim must be a two element vector in units of degrees');
end

if latlim(1)>latlim(2)
   eid = sprintf('%s:%s:invalidLatOrder', getcomp, mfilename);
   error(eid,'1st element of latlim must be greater than 2nd')
end

if lonlim(1)>=lonlim(2) 
   eid = sprintf('%s:%s:invalidLonOrder', getcomp, mfilename);
   error(eid,'1st element of lonlim must be greater than 2nd')
end

if latlim(1)<-90 || latlim(2)>90
   eid = sprintf('%s:%s:invalidLatLimits', getcomp, mfilename);
   error(eid,'latlim must be between -90 and 90')
end

if any(lonlim>360)
   eid = sprintf('%s:%s:invalidZeroTo2PiLonLimits', getcomp, mfilename);
   error(eid,'lonlim must be between 0 and 360')
end

if any(lonlim<0) && (any(lonlim<-180) || any(lonlim>180))
   eid = sprintf('%s:%s:invalidLonLimits', getcomp, mfilename);
   error(eid,'lonlim must be between -180 and 180')
end

%----------------------------------------------------------------------------
function  file = getFileFromDir(inputDir)
% Search directory for files in this order:
%   ETOPO2*
%   ETOPO5.DOS
%   ETOPO5.DAT
%   etopo5.northern.bat etopo5.southern.bat

list = dir(inputDir); 
index = find(~cellfun('isempty',regexpi({list(:).name},'etopo2.*.bin')));
if ~isempty(index)
   file{1} = [inputDir filesep list(index(1)).name];
   return
end

index = find(~cellfun('isempty',regexpi({list(:).name},'etopo5.dos')));
if ~isempty(index)
   file{1} = [inputDir filesep list(index(1)).name];
   return
end

index = find(~cellfun('isempty',regexpi({list(:).name},'etopo5.dat')));
if ~isempty(index)
   file{1} = [inputDir filesep list(index(1)).name];
   return
end

index = find(~cellfun('isempty',regexpi({list(:).name},'etopo5.*.bat')));
if ~isempty(index)
   % sort etopo5.northern.bat before etopo5.southern.bat
   file=sort({list(index).name});
   for i=1:numel(file)
       file{i} = [inputDir filesep file{i}];
   end
   return
end

% Error: no data in directory
eid = sprintf('%s:%s:noDataFound', getcomp, mfilename);
msg = sprintf('Directory "%s" does not contain ETOPO data.', inputDir);
error(eid, '%s', msg);
   
%----------------------------------------------------------------------------
function [map,refvec] = etopo5AsciiRead(file,samplefactor,latlim,lonlim)
%ETOPO5ASCIIREAD Read the ASCII file version of ETOPO5

%  Ascii data files (N & S hemispheres)
%  Data arranged in W-E columns (0 to 360) by N-S rows (90 to -90).
%  Elevation in meters

if all(latlim == [-90 90]) && (all(lonlim == [-180 180]) || all(lonlim == [0 360]))
   subset = 0;
else
   subset = 1;
end

sf = samplefactor;
dcell = 5/60;			% 5 minute grid
shift = 0;

if numel(file) == 1
   [path, name] = fileparts(file{1});
   if isempty(path)
      path = pwd;
   end
   if isequal(name,'etopo5.northern')
      file{2} = [path filesep 'etopo5.southern.bat'];
   else
      file{2} = file{1};
      file{1} = [path filesep 'etopo5.northern.bat'];
   end
end

if ~subset
   %  Check to see if samplefactor fits matrix dimensions
   if mod(1080,sf)~=0 || mod(4320,sf)~=0
      eid = sprintf('%s:%s:invalidSamplefactor', getcomp, mfilename);
      error(eid,'SAMPLEFACTOR must divide evenly into 1080 rows and 4320 columns')
   end
   
else

   %  Check to see if data needs to be shifted (-pi to pi)
   if lonlim(1)<0
      shift = 1;
   end

   %  Convert lat and lon limits to row and col limits
   if latlim(2)==90
      rowlim(1) = 1;
   else
      rowlim(1) = floor(-12*(latlim(2)-90)) + 1;
   end
   if latlim(1)==-90
      rowlim(2) = 2160;
   else
      rowlim(2) = ceil(-12*(latlim(1)-90));
   end
   if ~shift
      lon0 = 0;
   else
      lon0 = -180;
   end
   if (~shift && lonlim(1)==0) || (shift && lonlim(1)==-180)
      collim(1) = 1;
   else
      collim(1) = floor(12*(lonlim(1)-lon0)) + 1;
   end
   if (~shift && lonlim(2)==360) || (shift && lonlim(2)==180)
      collim(2) = 4320;
   else
      collim(2) = ceil(12*(lonlim(2)-lon0));
   end

end

%  Calculate row and col indices
srow = [0; (25923:25924:27971995)'];	% start row position indicators
rowindx1 = [];
rowindx2 = [];
if ~subset
   rowindx1 = 1:sf:1080;
   rowindx2 = 1:sf:1080;
   colindx = 1:sf:4320;
   maptop = 90;
   mapleft = 0;
else
   if rowlim(1)<=1080 && rowlim(2)<=1080		% submap in N hemisphere
      rowindx1 = rowlim(1):sf:rowlim(2);
      rowindx2 = [];
   elseif rowlim(1)>=1080 && rowlim(2)>=1080	% submap in S hemisphere
      rowindx1 = [];
      rowindx2 = rowlim(1)-1080:sf:rowlim(2)-1080;
   elseif rowlim(1)<=1080 && rowlim(2)>=1080	% submap in both hemispheres
      rowindx1 = rowlim(1):sf:1080;
      row1 = sf - (1080-rowindx1(length(rowindx1)));
      rowindx2 = row1:sf:rowlim(2)-1080;
   end
   colindx = collim(1):sf:collim(2);
   maptop = 90 - dcell*(rowlim(1)-1);
   mapleft = dcell*(collim(1)-1);
   if shift
      mapleft = dcell*(collim(1)-1) - 180;
   end
end

%  Read ETOPO5 ascii image files
%  Read from bottom to top of map (first row of matrix is bottom of map)

%  Northern hemisphere file
if ~isempty(rowindx1)
   fid = fopen(file{1},'r');
   if fid==-1
      eid = sprintf('%s:%s:northernFileNotFound', getcomp, mfilename);
      msg = sprintf('Filename "%s" not found',file{1});
      error(eid, '%s', msg);
   end
   y = srow(rowindx1);
   new_m = length(y);
   for m=new_m:-1:1
      %if mod(m,10)==0
      %   disp(['N' num2str(m)])
      %end
      fseek(fid,y(m),'bof');
      temp = fscanf(fid,'%d',[1 4320]);
      % pad data if necessary
      if size(temp,2) ~= 4320
         numberMissing = 4320-size(temp,2);
         temp = [temp temp(1:numberMissing)];
      end
      if shift
         temp = [temp(2161:4320) temp(1:2160)];
      end
      mapN(new_m+1-m,:) = temp(colindx);
   end
   fclose(fid);
else
   mapN = [];
end

%  Southern hemisphere file
if ~isempty(rowindx2)
   fid = fopen(file{2},'r');
   if fid==-1
      eid = sprintf('%s:%s:southernFileNotFound', getcomp, mfilename);
      msg = sprintf('Filename "%s" not found',file{2});
      error(eid, '%s', msg);
   end
   y = srow(rowindx2);
   new_m = length(y);
   for m=new_m:-1:1
      %if mod(m,10)==0
      %   disp(['S' num2str(m)])
      %end
      fseek(fid,y(m),'bof');
      temp = fscanf(fid,'%d',[1 4320]);
      if size(temp,2) ~= 4320
         numberMissing = 4320-size(temp,2);
         temp = [temp temp(1:numberMissing)];
      end
      if shift
         temp = [temp(2161:4320) temp(1:2160)];
      end
      mapS(new_m+1-m,:) = temp(colindx);
   end
   fclose(fid);
else
   mapS = [];
end

map = [mapS; mapN];
cellsize = (sf*dcell);
refvec = [1/cellsize maptop+cellsize/2  mapleft-cellsize/2];

%----------------------------------------------------------------------------
function [map,refvec] = etopo5BinaryRead(filename, samplefactor,latlim,lonlim)
%ETOPO5BINARYREAD read the binary version of ETOPO5 from ETOPO5.DOS

if any(lonlim<0) || any (lonlim>360)
   eid = sprintf('%s:%s:invalidEtopo5LonLimits', getcomp, mfilename);
   error(eid,'lonlim must be between 0 and 360')
end
fileList = dir(filename{1});
if isempty(fileList) || fileList(1).bytes ~= 18662400
   if isempty(fileList)
      fileSize = 0;
   else
      fileSize = fileList(1).bytes;
   end
   msg = sprintf('%s\n%s%s%s', ...
                'Input file size does not match expected file size. ', ...
                'Expected 18662400 bytes but file size is ', ...
                num2str(fileSize), ' bytes.');
   eid = sprintf('%s:%s:unexpectedFileSize', getcomp, mfilename);
   error(eid,msg);
end
    
% Hardwired information about the file format
ncols = 4320;
nrows = 2160;
cellsize = 5/60;

[path, name, ext] = fileparts(filename{1});
if isequal(lower([name ext]),'new_etopo5.bil')
   %    machineformat = 'ieee-be';
   %    lono = -180;
   %    lonlim = lonlim - 180;
   msg = sprintf('%s\n%s', ...
      'Reading new_etopo5.bil is no longer supported.', ...
      'Please refer to the tech note for obtaining ETOPO5 data.' );
   error('map:etopo:unsupportedFileName', msg);
elseif isequal(lower(ext),'.dat')
   machineformat = 'ieee-be';
   lono = 0;
else
   machineformat = 'ieee-le';
   lono = 0;
end
[map,refvec] = etopoBinaryRead(filename, samplefactor,latlim,lonlim, ...
                                  nrows, ncols, cellsize,lono,machineformat);

%----------------------------------------------------------------------------
function [map,refvec] = etopo2BinaryRead(filename, samplefactor,latlim,lonlim)
%ETOPO2BINARYREAD read the binary version of ETOPO2 

% Change default
if lonlim == [0 360]
   lonlim = [-180 180];
end
if any(lonlim<-180) || any(lonlim>180)
   eid = sprintf('%s:%s:invalidEtopo2LonLimits', getcomp, mfilename);
   error(eid,'lonlim must be between -180 and 180')
end

% Hardwired information about the file format
ncols = 10800;
nrows = 5400;
cellsize = 2/60;
if ~isempty(findstr('raw',lower(filename{1})))
   machineformat = 'ieee-be';
else
   machineformat = 'ieee-le';
end
lono = -180;
[map,refvec] = etopoBinaryRead(filename, samplefactor,latlim,lonlim, ...
                                  nrows, ncols, cellsize,lono,machineformat);


%--------------------------------------------------------------------------
function [map,refvec] = etopoBinaryRead(filename, samplefactor, ...
   latlim,lonlim, nrows, ncols, cellsize, lon11, machineformat)
%ETOPOBINARYREAD Read binary ETOPO data

% Null data value
NODATA_value = -32768;

% Other information about the file
precision = 'int16';

% offset from cell to cell
offset = samplefactor*cellsize;

% Create a referencing matrix of the entire data
% in order to obtain the row and column indices
% in which to read.
%
% Cell-center latitude  corresponding to data at (1,1)
% (South pole data is not included in ETOPO)
% lon11 is the cell-center longitude corresponding to data at (1,1)
lat11 = 90.0;
 
% Distance between latitude and longitude pixels of the data file.
dLat = -cellsize;  % From row to row moving south by cellsize degree
dLon =  cellsize;  % From column to column moving east by cellsize degree

% Create the referencing matrix
R = makerefmat(lon11, lat11, dLon, dLat);
                             
% Convert lat and lonlim to column and row indices
% (prevent longitude wrapping, subtract cellsize from lonlim max)
if lonlim(2) > 360 - cellsize || (lon11 == -180 && lonlim(2) > 180 -cellsize)
   lonlim(2) = lonlim(2) - cellsize;
end
[rlim, clim] = latlon2pix(R,latlim,lonlim);
rlim = round(rlim);
clim = round(clim);

% Ensure matrix coordinates are within limits
rlim = [max([1,min(rlim)]) min([max(rlim),nrows])];
rlim = sort(flipud(rlim(:))');
clim = [max([1,min(clim)]) min([max(clim),ncols])];
clim = sort(flipud(clim(:))');

readrows = rlim(1):samplefactor:rlim(2);
readcols = clim(1):samplefactor:clim(2);

% Ensure that we don't try to read beyond data limits
% if samplefactor = 1, we try to read an additional column
readcols = mod(readcols,ncols); 
readcols(readcols == 0) = ncols;
if size(readcols,2) > ncols
   readcols = readcols(1:ncols)
end

% Extract the data matrix
map = readmtx(filename{1},nrows,ncols,precision,readrows,readcols,machineformat);

% Re-read the data if the endianness is incorrect
if max(map(:)) > 20000 || min(map(:)) < -20000
  % Invalid values, flip endianness
  fromFormat = machineformat;
  if isequal(machineformat,'ieee-be')
     machineformat = 'ieee-le';
  else
     machineformat = 'ieee-be';
  end
  wid = sprintf('%s:%s:switchMachineFormat', getcomp, mfilename);
  msg = sprintf('%s\n%s%s%s%s\n%s%s%s', ...
                'Data file does not conform to standard naming conventions.', ...
                'Switching machine format from ', fromFormat, ' to ', machineformat, ...
                'and re-reading the file "',filename{1},'".');
  warning(wid, '%s', msg);
  map = readmtx(filename{1},nrows,ncols,precision,readrows,readcols,machineformat);
end

% Flip the data and set no data values
map = flipud(map);
map(map==NODATA_value) = NaN;

% Construct the referencing vector
% Find the min latitude and longitude that is read from the file
[readlat,readlon] = pix2latlon(R,[min(readrows) max(readrows)], ...
                                 [min(readcols) max(readcols)]);

% Create a new referencing matrix
% data is now from south to north and east to west
R = makerefmat(min(readlon), min(readlat), offset, offset);

% Convert to a referencing vector
refvec = refmat2vec(R,size(map));