function md=radarpower(md,highres)
%RADARPOWER - overlay a power radar image on an existing mesh
%
%   This routine will overlay a power radar image on an existing mesh.
%   The power amplitude will be output to vel for now.
%   In the future, think about a field to hold this value.
%
%   Usage:
%      md=radarpower(md,highres)

global ISSM_DIR
global MODELDATA

%find gdal coordinates
x0=min(md.x);
x1=max(md.x);

y0=min(md.y);
y1=max(md.y);

%Get path  to gdal binaries
path_gdal=[ISSM_DIR '/Libs/gdal/install'];

%Was gdal compiled? 
if exist([path_gdal 'gdal_translate'])==0,
	error(['radarpower error message: GDAL library needs to be compiled to use this routine. Compile GDAL in ' ISSM_DIR '/Libs/gdal to use this routine.']);
end

%the geotiff image is either 200m or 1km accuracy. 
if highres,
	geotiff_name=[MODELDATA 'MosaicTiffRsat/amm125m_v2_200m.tif'];
else
	geotiff_name=[MODELDATA 'MosaicTiffRsat/amm125m_v2_1km.tif'];
end

%Name of image
inputname='./temp.tif';

command=['system(''' path_gdal 'gdal_translate -quiet -projwin ' num2str(x0) ' ' num2str(y1) ' ' num2str(x1) ' ' num2str(y0) ' ' geotiff_name ' ' inputname ''');'];
eval(command);

%Read in temp.tif:
md.sarpwr=double(flipud(imread('temp.tif','TIFF')));
md.sarxm=(x0:(x1-x0)/(size(md.sarpwr,2)-1):x1);
md.sarym=(y0:(y1-y0)/(size(md.sarpwr,1)-1):y1);

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