function [z,refvec] = gtopo30(varargin)
%GTOPO30 Read 30-Arc-Sec global digital elevation data.
%
%   [Z, REFVEC] = GTOPO30(TILENAME) reads the GTOPO30 tile specified by
%   TILENAME and returns the result as a regular data grid. TILENAME is a
%   string which does not include an extension and indicates a GTOPO30 tile
%   in the current directory or on the MATLAB path. If TILENAME is empty or
%   omitted, a file browser will open for interactive selection of the
%   GTOPO30 header file. The data is returned at full resolution with the
%   latitude and longitude limits determined from the GTOPO30 tile. The
%   data grid, Z, is returned as an array of elevations.  Elevations are
%   given in meters above mean sea level using WGS84 as a horizontal datum.
%   REFVEC is the associated referencing vector.
%
%   [Z, REFVEC] = GTOPO30(TILENAME, SAMPLEFACTOR) reads a subset of the
%   elevation data from TILENAME. SAMPLEFACTOR is a scalar integer, which
%   when equal to 1 reads the data at its full resolution. When
%   SAMPLEFACTOR is an integer n greater than one, every nth point is read.
%   If SAMPLEFACTOR is omitted or empty, it defaults to 1. 
%   
%   [Z, REFVEC] = GTOPO30(TILENAME, SAMPLEFACTOR, LATLIM, LONLIM) reads a
%   subset of the elevation data from TILENAME. 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. Longitude limits range from [-180 180]. If
%   LATLIM or LONLIM is empty, the coordinate limits are determined from
%   the file.
%  
%   [Z, REFVEC] = GTOPO30(DIRNAME, SAMPLEFACTOR, LATLIM, LONLIM) reads and
%   concatenates data from multiple tiles within a GTOPO30 CD-ROM or
%   directory structure. The DIRNAME input is a string with the name of the
%   directory which contains the GTOPO30 tile directories or GTOPO30 tiles.
%   Within the tile directories are the uncompressed data files. The
%   DIRNAME for CD-ROMs distributed  by the USGS is the device name of the
%   CD-ROM drive. SAMPLEFACTOR if omitted or empty defaults to 1. LATLIM if
%   omitted or empty defaults to [-90 90]. LONLIM if omitted or empty
%   defaults to [-180 180]. 
%  
%   For details on locating GTOPO30 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#gtopo30') 
%   ">http://www.mathworks.com/support/tech-notes/2100/2101.html</a>
% 
%
%   Example 1
%   ---------
%   % Extract and display a subset of full resolution data for the state of
%   % Massachusetts.
%   % Read the stateline polygon boundary and calculate boundary limits.
%   Massachusetts = shaperead('usastatehi', 'UseGeoCoords', true, ...
%       'Selector',{@(name) strcmpi(name,'Massachusetts'), 'Name'});
%   latlim = [min(Massachusetts.Lat(:)) max(Massachusetts.Lat(:))];
%   lonlim = [min(Massachusetts.Lon(:)) max(Massachusetts.Lon(:))];
%   
%   % Read the GTOPO30 data at full resolution.
%   [Z,refvec] = gtopo30('W100N90',1,latlim,lonlim);
%
%   % Display the data grid and overlay the stateline boundary.
%   figure
%   usamap('Massachusetts');
%   geoshow(Z, refvec, 'DisplayType', 'surface')
%   geoshow([Massachusetts.Lat],[Massachusetts.Lon],'Color','magenta')
%   colormap(demcmap(Z))
%
%   Example 2
%   ---------
%   % Extract every 20th point from a tile. 
%   % Provide an empty filename and select the file interactively. 
%   [z,refvec] = gtopo30([],20);
%
%   Example 3
%   ---------
%   % Extract data for Thailand, an area which straddles two tiles. 
%   % The data is on CD number 3 distributed by the USGS.
%   % The CD-device is 'F:\'
%   latlim = [5.22 20.90]; 
%   lonlim = [96.72 106.38];
%   gtopo30s(latlim,lonlim)
%   % Extract every fifth data point for Thailand.
%   [Z,refvec] = gtopo30('F:\',5,latlim,lonlim);
%
%   Example 4
%   ---------
%   % Extract every 10th point from a column of data 5 degrees around the 
%   % prime meridian. The current directory contains GTOPO30 data.
%   [Z, refvec] = gtopo30(pwd, 10, [], [-5 5]);
%
%   See also GTOPO30S, GLOBEDEM, DTED, SATBATH, TBASE, USGSDEM.

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

%   Ascii header file, binary
%   Data arranged in W-E columns by N-S rows
%   Elevation in meters
%   GTOPO30 files are binary. No line ending conversion should be performed
%   during transfer or decompression.

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

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

% Read the file or directory
if exist([name '.DEM'],'file') == 2
    [z,refvec] = gtopo30FileRead(name,samplefactor,latlim,lonlim);
elseif exist(name,'dir') == 7
    [z,refvec] = gtopo30DirRead(name, samplefactor, latlim, lonlim);
else
    [z,refvec] = gtopo30FileRead(name,samplefactor,latlim,lonlim);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [name, samplefactor, latlim, lonlim] = parseInputs(varargin)
% Parse the command line inputs.

% Set default values
name = '';
samplefactor = [];

% Allow empty args
latlim = [];
lonlim = [];

% Obtain input arguments
switch nargin
    case 1
        name = varargin{1};
    case 2
        name = varargin{1};
        samplefactor = varargin{2};
    case 3
        name = varargin{1};
        samplefactor = varargin{2};
        latlim = varargin{3};
    case 4
        name = varargin{1};
        samplefactor = varargin{2};
        latlim = varargin{3};
        lonlim = varargin{4};
end
if isempty(samplefactor)
    samplefactor = 1;
else
    checkinput(samplefactor,{'double'},{'scalar','integer','finite'}, ...
               mfilename,'SAMPLEFACTOR',2);
end
if ~isempty(name) && ~ischar(name)
    eid = sprintf('%s:%s:invalidType', getcomp, mfilename);
    msg = sprintf('%s\n%s\n\n%s%s', ...
        'Function GTOPO30 expected its first input, TILENAME (or DIRNAME),', ...
        'to be type char','Instead its type was ',class(name));
    error(eid, '%s', msg);
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [z,refvec] = gtopo30DirRead(rd,samplefactor,latlim,lonlim)
% Read and concatenate GTOPO30 (30-arc-sec resolution) digital
% elevation files from a directory.

% Verify the latitude and longiutde limits
if isempty(latlim)
    latlim = [-90 90];
end
if isempty(lonlim)
    lonlim = [-180 180];
end
checklatlonlim(latlim,lonlim,'LATLIM','LONLIM',3,4);

% Get the bounding box filenames, dateline flag, and indices.
[names, XMIN, XMAX, YMIN, YMAX] = getLatLonLimits;
[fname, dateline, do] = getBoundingBoxFiles(latlim, lonlim, names, ...
                                            XMIN, XMAX, YMIN, YMAX);

% Add filesep to directory input if needed
if ~isequal(rd(end),filesep)
    rd(end+1) = filesep;
end

% Check the current directory for files
d = dir([rd '*.DEM']);
haveDirFiles = true;
if ~isempty(d)
    [ffname{1:length(d)}] = deal(d(:).name);
    for i=1:length(ffname)
        ffname{i} = ffname{i}(1:length(ffname{i})-4);
    end
    for i = length(do):-1:1
        if ~any(strcmp(fname{i},ffname))
            fname(i) = [];
            do(i) = [];
        end
    end
else
    % DEM files do not exist in current directory,
    % Append root directory and check to see if required files exist
    haveDirFiles = false;
    for i = length(do):-1:1
        ffname = [rd,fname{i},filesep,fname{i}];
        if ~(exist([ffname,'.DEM'],'file') == 2)
            ffname = [rd,fname{i}];
            if ~(exist([ffname,'.DEM'],'file') == 2)
                fname(i) = [];
                do(i) = [];
            else
                haveDirFiles = true;
            end
        end
    end
end

% Caculate the required data grid
[nrows,ncols] = calcDataGridSize(latlim, lonlim, samplefactor);

% Sort order over which files will be read
lon = []; lat = [];
for i = 1:length(do)
    lond = fname{i}(1);
    if lond == 'W'
        lonv = -1*str2double(fname{i}(2:4));
    else
        lonv = str2double(fname{i}(2:4));
    end
    lon = [lon lonv];
    latd = fname{i}(5);
    if latd == 'S'
        latv = -1*str2double(fname{i}(6:7));
    else
        latv = str2double(fname{i}(6:7));
    end
    lat = [lat latv];
end

% Calculate unique latitudes and longitudes
uniquelons = sort(unique(lon));
uniquelats = fliplr(sort(unique(lat)));
if isempty(uniquelons) || isempty(uniquelats)
    z = fillNaNDataGrid(nrows, ncols);
    refvec = calcRefVec(latlim, lonlim, samplefactor);
    wid = sprintf('%s:%s:noDataInDir', getcomp, mfilename);
    msg = sprintf('%s%s%s\n%s\n','The GTOPO30 directory "', rd, '" ', ...
        'does not contain data for the specified latitude longitude limits.');
    warning(wid,'%s',msg);
    return
end

nfile = cell(length(uniquelats), length(uniquelons));
for i = 1:length(uniquelats)
    for j = 1:length(uniquelons)
        switch sign(uniquelons(j))
            case 1
                lonh = 'E';
            case 0
                lonh = 'W';
            case -1
                lonh = 'W';
        end
        lonv = num2str(uniquelons(j));
        lonv = strrep(lonv,'-','');
        while length(lonv)< 3
            lonv = ['0',lonv];
        end

        % Create part of the file name referring to the longitude
        lonh = [lonh,lonv];
        switch sign(uniquelats(i))
            case 1
                lath = 'N';
            case -1
                lath = 'S';
        end

        % Create part of the file name referring to the latitude
        latv = num2str(uniquelats(i));
        latv = strrep(latv,'-','');
        lath = [lath,latv];

        % Create full file name
        nfile{i,j} = [lonh,lath];
        if haveDirFiles
            nfile{i,j} = [rd,nfile{i,j}];
        else
            nfile{i,j} = [rd,nfile{i,j},filesep,nfile{i,j}];
        end
    end
end

if length(uniquelats) == 1 && length(uniquelons) == 1 
    % single tile
    [z, refvec] = readDataGrid(nfile{1,1}, samplefactor,latlim,lonlim, ...
                                    nrows, ncols, dateline);

elseif length(uniquelats) == 1 && length(uniquelons) > 1
    % single row
    [z, refvec] = readDataGrid(nfile{1,1},samplefactor,latlim,lonlim, ...
                                    nrows, ncols, dateline);
    for j = 2:length(uniquelons)
        z = fillDataGrid(z, refvec, nfile{1,j},  samplefactor, ...
                         latlim, lonlim, dateline);
    end

elseif length(uniquelats) > 1 && length(uniquelons) == 1
    % single column
    [z, refvec] = readDataGrid(nfile{1,1},samplefactor,latlim,lonlim, ...
                                    nrows, ncols, dateline);
    for j = 2:length(uniquelats)
        z = fillDataGrid(z, refvec, nfile{j,1},  samplefactor, ...
                           latlim, lonlim, dateline);
    end

elseif length(uniquelats) > 1 && length(uniquelons) > 1
    % matrix
    z = fillNaNDataGrid(nrows, ncols);
    refvec = calcRefVec(latlim, lonlim, samplefactor);

    for i = 1:length(uniquelats)
        for j = 1:length(uniquelons)
            try
                z = fillDataGrid(z, refvec, nfile{i,j},  samplefactor, ...
                                   latlim, lonlim, dateline);
            catch
                [lastmsg, lastid] = lasterr;
                if ~isequal(lastid, 'map:gtopo30:fileNotFound')
                    rethrow(lasterror)
                end
            end
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [z,refvec] = readAcrossDateline(fname,samplefactor,latlim,lonlim)
% Read data across the dataline

lmin = lonlim(1);
lmax = lonlim(2);
% -180 to lonmax
lonlim(1) = -180;
lonlim(2) = lmax;

z = [];
refvec_left = [];
refvec_right = [];

try
    [z,refvec_left] = readDataFiles(fname,samplefactor,latlim,lonlim);
catch
end

% lon-min to 180
lonlim(1) = lmin;
lonlim(2) = 180;
try
    [tz,refvec_right] = readDataFiles(fname,samplefactor,latlim,lonlim);
    z = [z tz];
catch
end

if ~isempty(refvec_left)
    refvec = refvec_left;
elseif ~isempty(refvec_right)
    refvec = refvec_right;
else
    refvec = [];
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [z,refvec] = gtopo30FileRead(fname,samplefactor,latlim,lonlim)
% Read data from a GTOPO30 file.

% If the requested filename is empty, use the filebrowser to obtain the name
if isempty(fname)
    [filename, path] = uigetfile('*.HDR', 'select the GTOPO30 header file (*.HDR)');
    if filename == 0
        z = [];
        refvec = [];
        return
    end
    fname = ([path filename]);
    fname = fname(1:length(fname)-4);
end

% Does the DEM file exist?
if ~exist([fname,'.DEM'],'file')
    eid = sprintf('%s:%s:fileNotFound', getcomp, mfilename);
    fileMsg = sprintf('%s%s%s%s', 'GTOPO30 file "',  fname,'.DEM"', ...
        ' does not exist.');
    error(eid, '%s', fileMsg);
end

% Verify latitude and longitude limits
[names, XMIN, XMAX, YMIN, YMAX] = getLatLonLimits;
[path, name] = fileparts(fname);
index = strcmp(name,names);

if isempty(latlim)
    latlim = [YMIN(index) YMAX(index)];
    if isempty(latlim)
        latlim = [-90 90];
    end
end
if isempty(lonlim)
    lonlim = [XMIN(index) XMAX(index)];
    if isempty(lonlim)
        lonlim = [-180 180];
    end
end
checklatlonlim(latlim,lonlim,'LATLIM','LONLIM',3,4);

% Get the bounding box list of filenames  and output image size
[filenames, dateline] = getBoundingBoxFiles(latlim, lonlim, names, ...
                                            XMIN, XMAX, YMIN, YMAX);
[nrows,ncols] = calcDataGridSize(latlim, lonlim, samplefactor);

% Check if the requested filename is within the bounding-box
if ~any(strcmp(name,filenames))
    z = fillNaNDataGrid(nrows, ncols);
    refvec = calcRefVec(latlim, lonlim, samplefactor);

    wid = sprintf('%s:%s:noData', getcomp, mfilename);
    msg = sprintf('%s%s%s\n%s\n','The GTOPO30 file, ', fname,  ...
        ', does not contain data for the specified latitude longitude limits.', ...
        'Data for the latitude and longitude limits can be found in files: ');
    files = sprintf('%s\n',filenames{:});
    warning(wid,'%s',[msg files]);
    return
end

% Read the specified file
[z, refvec] = readDataGrid(fname, samplefactor,latlim,lonlim, ...
                           nrows, ncols, dateline);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [z,refvec] = readDataFiles(fname,samplefactor,latlim,lonlim)
%Read 30-Arc-Sec global digital elevation data from file.


%  ---------- Header Information ----------

%  Open ascii header file and read information
filename = [fname '.HDR'];
fid = fopen(filename,'r');
if fid==-1
    eid = sprintf('%s:%s:fileNotFound', getcomp, mfilename);
    msg = sprintf('GTOPO30 header file not found: %s', filename);
    error(eid, '%s', msg)
end

nrows = [];
ncols = [];
nodata = [];
ulxmap = [];
ulymap = [];
xdim = [];
ydim = [];

eof = 0;
while ~eof
    str = fscanf(fid,'%s',1);
    switch lower(str)
        case 'nrows', nrows = fscanf(fid,'%d',1);
        case 'ncols', ncols = fscanf(fid,'%d',1);
        case 'nodata', nodata = fscanf(fid,'%d',1);
        case 'ulxmap', ulxmap = fscanf(fid,'%f',1);
        case 'ulymap', ulymap = fscanf(fid,'%f',1);
        case 'xdim', xdim = fscanf(fid,'%f',1);
        case 'ydim', ydim = fscanf(fid,'%f',1);
        case '', eof = 1;
        otherwise, fscanf(fid,'%s',1);
    end
end
fclose(fid);

% Some of the data we wanted wasn't in the hdr file.
% Read the world file  to get it

if length([nrows ncols nodata ulxmap ulymap xdim ydim]) < 7
    filename = [fname '.BLW'];
    fid = fopen(filename,'r');
    if fid==-1
        filename = [fname '.DMW'];
        fid = fopen(filename,'r');
        if fid==-1
            eid = sprintf('%s:%s:ancillaryFileNotFound', getcomp, mfilename);
            msg = sprintf( ...
                'The GTOPO30 ancillary file (.BLW or .DMW) is not found for prefix: %s', ...
                fname);
            error(eid, '%s', msg)
        end
    end

    xdim = fscanf(fid,'%f',1);
    fscanf(fid,'%f',1);
    fscanf(fid,'%f',1);
    ydim = fscanf(fid,'%f',1); ydim = -ydim;
    ulxmap = fscanf(fid,'%f',1);
    ulymap = fscanf(fid,'%f',1);
    fclose(fid);
end

% Any information still missing?
if length([nrows ncols nodata ulxmap ulymap xdim ydim]) < 7
    eid = sprintf('%s:%s:incompleteHeader', getcomp, mfilename);
    msg = sprintf( ...
        'Incomplete header file or change in header file format for %s', ...
        fname);
    error(eid, '%s', msg);
end

% Other information about the file
precision = 'int16';
machineformat = 'ieee-be';
lato = ulymap;
lono = ulxmap;
dlat = -ydim;
dlon = xdim;

% Convert lat and lonlim to column and row indices
[clim,rlim] = yx2rc(lonlim(:),latlim(:),lono,lato,dlon,dlat);

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

readrows = rlim(1):samplefactor:rlim(2);
readcols = clim(1):samplefactor:clim(2);
readcols = mod(readcols,ncols); readcols(readcols == 0) = ncols;

% Extract the z matrix
filename = [fname '.DEM'];
z = readmtx(filename,nrows,ncols,precision,readrows,readcols,machineformat);
z = flipud(z);
z(z==nodata) = NaN;

% Construct the referencing vector.
[la1,lo1] = rc2yx(rlim,clim,lato,lono,dlat,dlon);
refvec = [abs(1/(dlat*samplefactor)) la1(1)-dlat/2 lo1(1)-dlon/2 ];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function refvec = calcRefVec(latlim, lonlim, samplefactor)
% Calculate the referencing vector

% 30 arc-seconds expressed in degrees.
resolution = 30/3600;

nrows = 6000;
ncols = 4800;

dlat = -resolution;
dlon = resolution;

% Use center of pixel
ulLat = latlim(2);
ulLon = lonlim(1);
lono = ulLon + resolution/2;
lato = ulLat - resolution/2;

% Convert lat and lonlim to column and row indices
[clim,rlim] = yx2rc(lonlim(:),latlim(:),lono,lato,dlon,dlat);

% 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])];

% Construct the referencing vector.
[la1,lo1] = rc2yx(rlim,clim,lato,lono,dlat,dlon);
refvec = [abs(1/(dlat*samplefactor)) la1(1)-dlat/2 lo1(1)-dlon/2 ];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [z refvec] = readDataGrid(filename, samplefactor, ...
                                   latlim, lonlim, nrows, ncols, dateline)
% Read the data grid from filename

% Calculate the refvector
refvec = calcRefVec(latlim, lonlim, samplefactor);

% Fill the data grid with values
z = fillDataGrid(fillNaNDataGrid(nrows, ncols), refvec,  filename, ...
                 samplefactor, latlim, lonlim, dateline);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function z  = fillNaNDataGrid(nrows, ncols)
% Fill a data grid with NaN values

try
    z = nan(nrows, ncols);
catch
    eid = sprintf('%s:%s:nanError', getcomp, mfilename);
    msg = sprintf('%s%d%s%d', 'Unable to create a data grid of size ', ...
                  nrows, ' by ', ncols);
    error(eid, '%s', msg);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function z  = fillDataGrid(z, refvec, nfile, samplefactor, ...
                           latlim, lonlim, dateline)
% Fill the data grid with values from nfile

if dateline
    [tz,trefvec] = readAcrossDateline(nfile, samplefactor, ...
                                      latlim, lonlim);
else
    [tz,trefvec] = readDataFiles(nfile, samplefactor,latlim, lonlim);
end

[row, col, tz] = getRowCol(tz, trefvec, z, refvec);
if ~isempty(row) && ~isempty(col)
    z(row: row+size(tz,1)-1, col:col+size(tz,2)-1) = tz;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [row, col, tz] = getRowCol(tz, trefvec, z, refvec)
% Get the row and column numbers

cr = [1:2;1:2];
cc = cr';
row = [];
col = [];
try
    [lat, lon] = pix2latlon(refvec2mat( trefvec, size(tz)), cr,  cc);
    [row, col] = latlon2pix(refvec2mat(  refvec, size( z)), lat, lon);
    row = round(row);
    col = round(col);
    rindex = min(find(row>0));
    cindex = min(find(col>0));
    index = max(min(cindex),min(rindex));
    row = row(index);
    col = col(index);
    if index ~= 1
        tz(cr(index):end, cc(index):end) = tz;
    end
catch
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [fnames, XMIN, XMAX, YMIN, YMAX] = getLatLonLimits
% Get the latitude and longitude limits of GTOPO30 tiles

% Open the GTOPO30 bounding box file
fid = fopen('gtopo30s.dat','r');
if fid==-1
    eid = sprintf('%s:%s:internalError', getcomp, mfilename);
    error(eid, '%s', 'Unable to open gtopo30s.dat')
end

% Preallocate bounding rectangle data for speed
YMIN = zeros(1,33); 
YMAX = YMIN;
XMIN = YMIN; 
XMAX = YMIN;

% Read names and bounding rectangle limits
for n=1:33
    fnames{n,1} = upper(fscanf(fid,'%s',1));
    YMIN(n) = fscanf(fid,'%d',1);
    YMAX(n) = fscanf(fid,'%d',1);
    XMIN(n) = fscanf(fid,'%d',1);
    XMAX(n) = fscanf(fid,'%d',1);
end
fclose(fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [filenames, dateline, do] = getBoundingBoxFiles(latlim, lonlim, fnames, ...
                                                         XMIN, XMAX, YMIN, YMAX)
% Get the bounding box filenames

% case where dateline is not crossed
if lonlim(1) <= lonlim(2)
    dateline = 0; % tile do not cross dateline
    do = ...
        find( ...
        (...
        (latlim(1) <= YMIN & latlim(2) >= YMAX) | ... % tile is completely within region
        (latlim(1) >= YMIN & latlim(2) <= YMAX) | ... % region is completely within tile
        (latlim(1) >  YMIN & latlim(1) <  YMAX) | ... % min of region is on tile
        (latlim(2) >  YMIN & latlim(2) <  YMAX)   ... % max of region is on tile
        ) ...
        &...
        (...
        (lonlim(1) <= XMIN & lonlim(2) >= XMAX) | ... % tile is completely within region
        (lonlim(1) >= XMIN & lonlim(2) <= XMAX) | ... % region is completely within tile
        (lonlim(1) >  XMIN & lonlim(1) <  XMAX) | ... % min of region is on tile
        (lonlim(2) >  XMIN & lonlim(2) <  XMAX)   ... % max of region is on tile
        )...
        );
else

    % case where the dateline is crossed
    % lonlim(1) > lonlim(2)
    dateline = 1; % tiles  cross dateline
    lmin = lonlim(1);
    lmax = lonlim(2);
    lonlim(2) = 180;

    % do eastern side of the dateline first
    doEAST = ...
        find( ...
        (...
        (latlim(1) <= YMIN & latlim(2) >= YMAX) | ... % tile is completely within region
        (latlim(1) >= YMIN & latlim(2) <= YMAX) | ... % region is completely within tile
        (latlim(1) >  YMIN & latlim(1) <  YMAX) | ... % min of region is on tile
        (latlim(2) >  YMIN & latlim(2) <  YMAX)   ... % max of region is on tile
        ) ...
        &...
        (...
        (lonlim(1) <= XMIN & lonlim(2) >= XMAX) | ... % tile is completely within region
        (lonlim(1) >= XMIN & lonlim(2) <= XMAX) | ... % region is completely within tile
        (lonlim(1) >  XMIN & lonlim(1) <  XMAX) | ... % min of region is on tile
        (lonlim(2) >  XMIN & lonlim(2) <  XMAX)   ... % max of region is on tile
        )...
        );
    % do western side of the dateline second
    lonlim(1) = -180; lonlim(2) = lmax;
    doWEST = ...
        find( ...
        (...
        (latlim(1) <= YMIN & latlim(2) >= YMAX) | ... % tile is completely within region
        (latlim(1) >= YMIN & latlim(2) <= YMAX) | ... % region is completely within tile
        (latlim(1) >  YMIN & latlim(1) <  YMAX) | ... % min of region is on tile
        (latlim(2) >  YMIN & latlim(2) <  YMAX)   ... % max of region is on tile
        ) ...
        &...
        (...
        (lonlim(1) <= XMIN & lonlim(2) >= XMAX) | ... % tile is completely within region
        (lonlim(1) >= XMIN & lonlim(2) <= XMAX) | ... % region is completely within tile
        (lonlim(1) >  XMIN & lonlim(1) <  XMAX) | ... % min of region is on tile
        (lonlim(2) >  XMIN & lonlim(2) <  XMAX)   ... % max of region is on tile
        )...
        );
    % concatenate indices
    do = [doEAST doWEST];
    % restore original values for lonlim
    lonlim(1) = lmin; 
    lonlim(2) = lmax;
end

% Set filenames to only valid files
if ~isempty(do)
    filenames = fnames(do);
else
    eid = sprintf('%s:%s:invalidLimits', getcomp, mfilename);
    error(eid, '%s', 'Invalid latitude or longitude limits.');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [nrows,ncols] = calcDataGridSize(latlim, lonlim, samplefactor)
% Calculate the size of the data grid

% 30 arc-seconds expressed in degrees
resolution = 30/3600;
extractedsize = sizem(latlim,lonlim,1/(samplefactor*resolution));
nrows = extractedsize(1);
ncols = extractedsize(2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function checklatlonlim( lat, lon, lat_var_name, lon_var_name, ...
                         lat_pos, lon_pos)

checkinput(lat,{'double'},{'real','vector','finite'},mfilename, ...
           lat_var_name,lat_pos);
checkinput(lon,{'double'},{'real','vector','finite'},mfilename, ...
           lon_var_name,lon_pos);

if numel(lon) ~= 2 || numel(lat) ~= 2
    eid = sprintf('map:%s:latlonSizeMismatch',mfilename);
    msg1 = sprintf('Function %s expected its %s and %s input arguments,', ...
                    upper(mfilename), num2ordinal(lat_pos), ...
                    num2ordinal(lon_pos));
    msg2 = sprintf('%s and %s, to be vectors of length 2.', ...
                   lat_var_name, lon_var_name);
    error(eid, '%s\n%s', msg1, msg2);
end

if lat(1) > lat(2)
    eid = sprintf('map:%s:latLimitError',mfilename);
    error(eid, '%s', 'Latitude limits must be ascending.');
end