Index: /issm/trunk-jpl/src/m/plot/plot_overlay.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_overlay.m	(revision 27605)
+++ /issm/trunk-jpl/src/m/plot/plot_overlay.m	(revision 27606)
@@ -42,4 +42,9 @@
 	xlim=getfieldvalue(options,'xlim',[min(x) max(x)])/getfieldvalue(options,'unit',1);
 	ylim=getfieldvalue(options,'ylim',[min(y) max(y)])/getfieldvalue(options,'unit',1);
+	if exist(options, 'axis');
+		myaxis = getfieldvalue(options,'axis');
+		xlim = [myaxis(1), myaxis(2)];
+		ylim = [myaxis(3), myaxis(4)];
+	end
 	options=addfielddefault(options,'xlim',xlim);
 	options=addfielddefault(options,'ylim',ylim);
@@ -48,17 +53,20 @@
 contrast = getfieldvalue(options,'contrast',1);  
 radar    = md.radaroverlay.pwr;
-if size(radar,3)>1,
-	disp('WARNING: color image converted to greyscale intensity image');
-	if strcmp(class(radar),'uint8'),
-		radar=double(sum(radar,3))/(255*3);
-	else
-		radar=sum(radar,3)/3;
+
+if ~radaronly
+	if size(radar,3)>1,
+		disp('WARNING: color image converted to greyscale intensity image');
+		if strcmp(class(radar),'uint8'),
+			radar=double(sum(radar,3))/(255*3);
+		else
+			radar=sum(radar,3)/3;
+		end
 	end
+	if getfieldvalue(options,'backgroundbtw',0)
+		radar(find(radar==0))=1; %Change background from black to white
+	end
+	radar = radar.^(contrast);
+	radar = radar./max(radar(:));
 end
-if getfieldvalue(options,'backgroundbtw',0)
-	radar(find(radar==0))=1; %Change background from black to white
-end
-radar = radar.^(contrast);
-radar = radar./max(radar(:));
 
 if getfieldvalue(options,'backgroundbtw',0)
@@ -101,5 +109,5 @@
 
 %Special colormaps that require hsv treatment
-colorm=getfieldvalue(options,'colormap','Rignot');
+colorm=getfieldvalue(options,'colormap','parula');
 if strcmpi(colorm,'Rignot') | strcmpi(colorm,'Seroussi') | strcmpi(colorm,'redblue')
 	if strcmpi(colorm,'Rignot'),
@@ -133,18 +141,22 @@
 	image_rgb=hsv2rgb(image_hsv);
 else
-	colorm = getcolormap(options);
-	len    = size(colorm,1);
+	if radaronly
+		image_rgb = radar;
+	else
+		colorm = getcolormap(options);
+		len    = size(colorm,1);
 
-	ind = ceil((len-1)*(data_grid-data_min)/(data_max - data_min + eps) +1);
-	ind(find(ind>len))=len;
-	image_rgb=zeros(size(data_grid,1),size(data_grid,2),3);
-	r=colorm(:,1); image_rgb(:,:,1)=r(ind); clear r;
-	g=colorm(:,2); image_rgb(:,:,2)=g(ind); clear g;
-	b=colorm(:,3); image_rgb(:,:,3)=b(ind); clear b;
+		ind = ceil((len-1)*(data_grid-data_min)/(data_max - data_min + eps) +1);
+		ind(find(ind>len))=len;
+		image_rgb=zeros(size(data_grid,1),size(data_grid,2),3);
+		r=colorm(:,1); image_rgb(:,:,1)=r(ind); clear r;
+		g=colorm(:,2); image_rgb(:,:,2)=g(ind); clear g;
+		b=colorm(:,3); image_rgb(:,:,3)=b(ind); clear b;
 
-	%Now add radarmap
-	r = image_rgb(:,:,1).*radar;  r(data_nan) = radar(data_nan);  image_rgb(:,:,1) = r;  clear r;
-	g = image_rgb(:,:,2).*radar;  g(data_nan) = radar(data_nan);  image_rgb(:,:,2) = g;  clear g;
-	b = image_rgb(:,:,3).*radar;  b(data_nan) = radar(data_nan);  image_rgb(:,:,3) = b;  clear b;
+		%Now add radarmap
+		r = image_rgb(:,:,1).*radar;  r(data_nan) = radar(data_nan);  image_rgb(:,:,1) = r;  clear r;
+		g = image_rgb(:,:,2).*radar;  g(data_nan) = radar(data_nan);  image_rgb(:,:,2) = g;  clear g;
+		b = image_rgb(:,:,3).*radar;  b(data_nan) = radar(data_nan);  image_rgb(:,:,3) = b;  clear b;
+	end
 end
 
Index: /issm/trunk-jpl/src/m/plot/radarpower.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/radarpower.m	(revision 27605)
+++ /issm/trunk-jpl/src/m/plot/radarpower.m	(revision 27606)
@@ -100,28 +100,58 @@
 	%}}}
 else %user provided image {{{
+
 	%user provided an image. check we also have overlay_xlim and overlay_ylim  options, to know what range of coordinates the image covers.
-	if (~exist(options,'overlay_xlim') | ~exist(options,'overlay_xlim')| ~exist(options,'overlay_xposting')| ~exist(options,'overlay_yposting')),
-		error('radarpower error message: please provide overlay_xlim, overlay_ylim, overlay_xposting and overlay_yposting options together with overlay_image option');
+	filename = getfieldvalue(options,'overlay_image');
+	[filepath,name,ext] = fileparts(filename);
+	if ~exist(filename)
+		error([filename ' not found']);
 	end
-	overlay_image=getfieldvalue(options,'overlay_image');
-	overlay_xlim=getfieldvalue(options,'overlay_xlim');
-	overlay_ylim=getfieldvalue(options,'overlay_ylim');
-	overlay_xposting=getfieldvalue(options,'overlay_xposting');
-	overlay_yposting=getfieldvalue(options,'overlay_yposting');
 
-	sizex=floor((x1-x0)/overlay_xposting);
-	sizey=floor((y1-y0)/overlay_yposting);
-	topleftx=floor((x0-overlay_xlim(1))/overlay_xposting); % x min
-	toplefty=floor((overlay_ylim(2)-y1)/overlay_yposting); % y max
+	%Is it a geotiff?
+	if strcmp(ext,'.tiff') || strcmp(ext,'.tif')
 
-	%Read and crop file
-	disp('Warning: expecting coordinates in polar stereographic (Std Latitude: 70ºN Meridian: 45º)');
-	im=imread(overlay_image);
-	%adjust contrast and brightness
-	%im=imadjust(im,[a b],[c d]);
-	im=im(toplefty:toplefty+sizey,topleftx:topleftx+sizex);
-	md.radaroverlay.pwr=double(flipud(im));
-	md.radaroverlay.x=(x0:(x1-x0)/(size(md.radaroverlay.pwr,2)-1):x1);
-	md.radaroverlay.y=(y0:(y1-y0)/(size(md.radaroverlay.pwr,1)-1):y1);
+		%Crop image from xylim
+		tempfilename='./temp.tif';
+		eval(['!gdal_translate -quiet -projwin ' num2str(x0) ' ' num2str(y1) ' ' num2str(x1) ' ' num2str(y0) ' ' filename ' ' tempfilename]);
+
+		%Read in temp.tif:
+		im=imread('temp.tif','TIFF');
+		%adjust contrast and brightness
+		%im=imadjust(im,[a b],[c d]);
+		pixelskip=max(1,ceil(posting/((x1-x0)/(size(im,2)))));
+		%md.radaroverlay.pwr=double(flipud(im(1:pixelskip:end,1:pixelskip:end,:)));
+		md.radaroverlay.pwr=double(im(1:pixelskip:end,1:pixelskip:end,:))/255;
+		md.radaroverlay.x=x0:(x1-x0)/(size(md.radaroverlay.pwr,2)-1):x1;
+		md.radaroverlay.y=y1:-(y1-y0)/(size(md.radaroverlay.pwr,1)-1):y0;
+
+		%Erase image or keep it?
+		if ~getfieldvalue(options,'keep_image',0),
+			delete(tempfilename);
+		end
+	else
+		if (~exist(options,'overlay_xlim') | ~exist(options,'overlay_xlim')| ~exist(options,'overlay_xposting')| ~exist(options,'overlay_yposting')),
+			error('radarpower error message: please provide overlay_xlim, overlay_ylim, overlay_xposting and overlay_yposting options together with overlay_image option');
+		end
+		overlay_xlim=getfieldvalue(options,'overlay_xlim');
+		overlay_ylim=getfieldvalue(options,'overlay_ylim');
+		overlay_xposting=getfieldvalue(options,'overlay_xposting');
+		overlay_yposting=getfieldvalue(options,'overlay_yposting');
+		overlay_image=getfieldvalue(options,'overlay_image');
+
+		sizex=floor((x1-x0)/overlay_xposting);
+		sizey=floor((y1-y0)/overlay_yposting);
+		topleftx=floor((x0-overlay_xlim(1))/overlay_xposting); % x min
+		toplefty=floor((overlay_ylim(2)-y1)/overlay_yposting); % y max
+
+		%Read and crop file
+		disp('Warning: expecting coordinates in polar stereographic (Std Latitude: 70ºN Meridian: 45º)');
+		im=imread(overlay_image);
+		%adjust contrast and brightness
+		%im=imadjust(im,[a b],[c d]);
+		im=im(toplefty:toplefty+sizey,topleftx:topleftx+sizex);
+		md.radaroverlay.pwr=double(flipud(im));
+		md.radaroverlay.x=(x0:(x1-x0)/(size(md.radaroverlay.pwr,2)-1):x1);
+		md.radaroverlay.y=(y0:(y1-y0)/(size(md.radaroverlay.pwr,1)-1):y1);
+	end
 end %}}}