function [z,lon,lat,inan] = srtmread(xlim,ylim,varargin)
% srtmread
%
% [z,lon,lat,ifill] = srtmread(xlim,ylim)
%
% [z,lon,lat,ifill] = srtmread(xlim,ylim,skip)

skip = 1;

datadir =  './srtm/';

if length(varargin) == 1
	skip = varargin{1};
end

% Array dimennsions
ii   = 1:skip:1201;
indx = 1:length(ii);
n    = length(indx);

% lon-lat
ilon = floor(xlim(1)):ceil(xlim(2)-1);
ilat = floor(ylim(1)):ceil(ylim(2)-1);
z = zeros(length(ilat).*n,length(ilon).*n,'uint16');
lon = linspace(ilon(1),ilon(end)+1,(n-1)*length(ilon) + 1);
lat = linspace(ilat(1),ilat(end)+1,(n-1)*length(ilat) + 1);

for i = 1:length(ilon)
for j = 1:length(ilat)
	% Get filename
	if ilat(j) > 0
		file = sprintf('%s/N%2.2dW%3.3d.hgt',datadir,abs(ilat(j)),abs(ilon(i)));
	else
		file = sprintf('%s/S%2.2dW%3.3d.hgt',datadir,abs(ilat(j)),abs(ilon(i)));
	end
	disp(file);

	% Read the data
	if exist(file,'file')
		fi = fopen(file,'r','ieee-be');

		if skip == 1
			x  = flipud(fread(fi,[1201,1201],'uint16')');
			z((indx) + (j-1)*(n-1),(indx) + (i-1)*(n-1)) = x(ii,ii);
		else
			% seems to work for most (but maybe not all) skips 
			for k = 1:skip:1201
				fseek(fi,(k-1)*1201*2,-1);
				ij = 1+ (j)*(n-1) - (k-1)./skip;
				z(ij, (indx) + (i-1)*(n-1)) =  fread(fi,n,'uint16',(skip-1).*2);
			end
		end
		fclose(fi);
	else 
		z((indx) + (j-1)*(n-1),(indx) + (i-1)*(n-1)) = 0;
	end
end
end

% 
%z(abs(z)==32768)=NaN;

% Fill in missing data regions
%[z,bw] = roifill(z, abs(z)==32768 );
z = single(z);
inan = z>10000;
z1   = z;
z2   = z;
for i = 1:size(z,1)
	if any(inan(i,:))
		igood = find(~inan(i,:));
		ibad =  find(inan(i,:));
		z1(i,ibad) = interp1(igood,z(i,igood),ibad,'linear','extrap');
	end
end
for i = 1:size(z,2)
	if any(inan(:,i))
		igood = find(~inan(:,i));
		ibad =  find(inan(:,i));
		z2(ibad,i) = interp1(igood,z(igood,i),ibad,'linear','extrap');
	end
end
z(inan) = 0.5.*(z1(inan)+z2(inan));

% Get data between limits
ix = find(lon>=xlim(1) & lon<=xlim(2));
iy = find(lat>=ylim(1) & lat<=ylim(2));
z = z(iy,ix);
lon = lon(ix);
lat = lat(iy);
inan = inan(iy,ix);
