function [xm,ym,modis]=modis(modisgeotif,xlim,ylim)
%MODIS - from modis geotiff, return image
%
%   Usage:
%      [xm,ym,modis]=modis(modisgeotif,xlim,ylim)
%

%find gdal coordinates
x0=min(xlim);
x1=max(xlim);

y0=min(ylim);
y1=max(ylim);

%Get path  to gdal binaries
path_gdal=[issmdir() '/externalpackages/gdal/install/bin/'];

%Was gdal compiled? 
if ~exist([path_gdal 'gdal_translate']),
	error(['modis error message: GDAL library needs to be compiled to use this routine. Compile GDAL in ' issmdir() '/externalpackages/gdal to use this routine.']);
end

inputname='./temp.tif';
system([path_gdal 'gdal_translate -quiet -projwin ' num2str(x0) ' ' num2str(y1) ' ' num2str(x1) ' ' num2str(y0) ' ' modisgeotif ' ' inputname ]);

%Read in temp.tif:
modis=double(flipud(imread('temp.tif','TIFF')));
xm=(x0:(x1-x0)/(size(md.radaroverlay.pwr,2)-1):x1);
ym=(y0:(y1-y0)/(size(md.radaroverlay.pwr,1)-1):y1);

%Erase image
system('rm -rf ./temp.tif');


end
