function [map,maplegend] = gtopo30(fname,scalefactor,latlim,lonlim);

%GTOPO30 30-Arc-Sec global digital elevation data extraction
%
% [map,maplegend] = GTOPO30(filename,scalefactor) reads the GTOPO30
% files and returns the result as a regular matrix map.  The filename is
% given as a string which does not include the extension.  If the
% files are not found on the Matlab path, they may be selected interactively.
% Scalefactor is an integer, which when equal to 1 gives the data
% at its full resolution.  When scalefactor is an integer n larger
% than one, every nth point is returned.  The scalefactor must divide
% evenly into the number of rows and columns in the data file. The map data
% is returned as an array of elevations and associated regular matrix map
% legend. Elevations are given in meters above mean sea level using WGS 84
% as a horizontal datum.
%
% [map,maplegend] = GTOPO30(fname,scalefactor,latlim,lonlim)
% allows a subset of the map data to be read. The limits of the desired
% data are specified as vectors of latitude and longitude in degrees.
% The elements of latlim and lonlim must be in ascending order.
%
% The data is available over the Internet via anonymous ftp from
% <ftp://edcftp.cr.usgs.gov/pub/data/gtopo30/global>. The data and
% some documentation is also available over the World-Wide-Web from
% <http://edcwww.cr.usgs.gov/landdaac/gtopo30/gtopo30.html> and
% <http://edcwww.cr.usgs.gov/landdaac/gtopo30/README.html>.
%
% See also: DTED, SATBATH, TBASE, USGSDEM

%  Copyright 1996-1998 by Systems Planning and Analysis, Inc. and The MathWorks, Inc.
% $Revision: 1.5 $
%  Written by:  A. Kim, W. Stumpf
%  Version 1.0 (11/30/96)


%  Ascii header file, binary
%  Data arranged in W-E columns by N-S rows
%  Elevation in meters

if nargin==2
	subset = 0;
elseif nargin==4
	subset = 1;
else
	error('Incorrect number of arguments')
end

sf = scalefactor;
tolerance = 1e-6;		% accounts for floating point error

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

%  Open ascii header file and read information

filename = [fname '.HDR'];
fid = fopen(filename,'r');
if fid==-1
	[filename, path] = uigetfile('*.HDR', 'select the GTOPO30 header file (*.HDR)');
	if filename == 0 ; return; end
	filename = [path filename];
	fid = fopen(filename,'r');
	fname = filename(1:length(filename)-4);
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
			[filename, path] = uigetfile('*', 'select the GTOPO30 world file (*.BLW or *.DMW)');
			if filename == 0 ; return; end
			filename = [path filename];
			fid = fopen(filename,'r');
		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);
end

% Any information still missing?

if length([nrows ncols nodata ulxmap ulymap xdim ydim]) < 7
	error('Incomplete header file or change in header file format')
end

%  Define border of map

maplatlim(1) = ulymap - (nrows-1)*ydim - ydim/2;
maplatlim(2) = ulymap + ydim/2;
maplonlim(1) = ulxmap - xdim/2;
maplonlim(2) = ulxmap + (ncols-1)*xdim + xdim/2;

linebreak = sprintf('\n');

if ~subset

%  Check to see if nrows and ncols fit scalefactor
	incl = mod((nrows-1),sf)==0 & mod((ncols-1),sf)==0;
	excl = mod(nrows,sf)==0 & mod(ncols,sf)==0;
	if ~incl & ~excl
		error(['samplefactor must divide evenly into ' num2str(nrows) ...
		       ' rows and ' num2str(ncols) ' columns'])
	end

else

%  Check to see if latlim and lonlim within map limits
	if  latlim(1)>maplatlim(2)+tolerance | ...
		latlim(2)<maplatlim(1)-tolerance | ...
		lonlim(1)>maplonlim(2)+tolerance | ...
		lonlim(2)<maplonlim(1)-tolerance
			warning([ ...
			'Requested latitude or longitude limits are off the map' linebreak ...
			' latlim for this dataset is ' ...
			 mat2str( [maplatlim(1) maplatlim(2)],3) linebreak ...
		    ' lonlim for this dataset is '...
			 mat2str( [maplonlim(1) maplonlim(2)],3) ...
			 ])
			 map=[];maplegend = [];
			 return
	end

	warn = 0;
	if latlim(1)<maplatlim(1)-tolerance ; latlim(1)=maplatlim(1);warn = 1; end
	if latlim(2)>maplatlim(2)+tolerance ; latlim(2)=maplatlim(2);warn = 1; end
	if lonlim(1)<maplonlim(1)-tolerance ; lonlim(1)=maplonlim(1);warn = 1; end
	if lonlim(2)>maplonlim(2)+tolerance ; lonlim(2)=maplonlim(2);warn = 1; end
	if warn
		warning([ ...
			'Requested latitude or longitude limits exceed map limits' linebreak ...
			' latlim for this dataset is ' ...
			 mat2str( [maplatlim(1) maplatlim(2)],3) linebreak ...
		    ' lonlim for this dataset is '...
			 mat2str( [maplonlim(1) maplonlim(2)],3) ...
			 ])
	end

%  Convert lat and lon limits to row and col limits
	invydim = 1/ydim;
	if abs(latlim(2)-maplatlim(2))<tolerance
		rowlim(1) = 1;
	else
		temp = -invydim*(latlim(2)-maplatlim(2));
		if abs(temp-round(temp))<tolerance
			rowlim(1) = round(temp) + 1;
		else
			rowlim(1) = floor(temp) + 1;
		end
	end
	if abs(latlim(1)-maplatlim(1))<tolerance
		rowlim(2) = nrows;
	else
		temp = -invydim*(latlim(1)-maplatlim(2));
		if abs(temp-round(temp))<tolerance
			rowlim(2) = round(temp);
		else
			rowlim(2) = ceil(temp);
		end
	end
	invxdim = 1/xdim;
	if abs(lonlim(1)-maplonlim(1))<tolerance
		collim(1) = 1;
	else
		temp = invxdim*(lonlim(1)-maplonlim(1));
		if abs(temp-round(temp))<tolerance
			collim(1) = round(temp) + 1;
		else
			collim(1) = floor(temp) + 1;
		end
	end
	if abs(lonlim(2)-maplonlim(2))<tolerance
		collim(2) = ncols;
	else
		temp = invxdim*(lonlim(2)-maplonlim(1));
		if abs(temp-round(temp))<tolerance
			collim(2) = round(temp);
		else
			collim(2) = ceil(temp);
		end
	end

end

%  ---------- Elevation Data ----------

%  Open binary image file and read data
filename = [fname '.BIL'];
fid = fopen(filename,'rb','ieee-be');
if fid==-1
	filename = [fname '.DEM'];
	fid = fopen(filename,'rb','ieee-be');
	if fid==-1
		error('GTOPO30 binary image file not found')
	end
end
if ~subset
	firstrow = 0;
	lastrow = 2*sf*floor((nrows-1)/sf)*ncols;
	colindx = 1:sf:sf*floor((ncols-1)/sf)+1;
	maptop = maplatlim(2);
	mapleft = maplonlim(1);
else
	firstrow = 2*(rowlim(1)-1)*ncols;
	lastrow = firstrow + 2*sf*floor((rowlim(2)-rowlim(1))/sf)*ncols;
	colindx = collim(1):sf:collim(1)+sf*floor((collim(2)-collim(1))/sf);
	maptop =  ulymap - ydim*(rowlim(1)-1) + ydim/2;
	mapleft =  ulxmap + xdim*(collim(1)-1) - xdim/2;
end
srow = firstrow:2*sf*ncols:lastrow;		% start row position indicators
rows = length(srow);
%  Read from bottom to top of map (first row of matrix is bottom of map)
for m=rows:-1:1
%	if mod(m,10)==0
%		disp(m)
%	end
	fseek(fid,srow(m),'bof');
	temp = fread(fid,[1,ncols],'int16');
	temp(find(temp==nodata)) = nan;
	map(rows-m+1,:) = temp(colindx);
end
maplegend = [1/(sf*xdim) maptop mapleft];
fclose(fid);
