Index: /issm/trunk-jpl/src/m/plot/radarpower.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/radarpower.m	(revision 20486)
+++ /issm/trunk-jpl/src/m/plot/radarpower.m	(revision 20487)
@@ -34,6 +34,5 @@
 y0=min(ylim); y1=max(ylim);
 
-%figure out if we should go look for Greenland or Antarctica geotiff, or if user provided one.
-if ~exist(options,'overlay_image'),
+if ~exist(options,'overlay_image'), % no image provided, go look into ModelData for one!{{{
 	if md.mesh.epsg==3413,   %Greenland 
 		%old code {{{ 
@@ -143,21 +142,44 @@
 
 	else
-		error('EPSG code not supported yet (ex: 3413 for UPS Greenland, 3031 for UPS Antarctica)');
-	end
-else
+		if ~exist(options,'geotiff_name'), 
+			error('Need a geotiff for areas outside of Greenland and Antarctica');
+		end
+		
+		%Name of image
+		inputname='./temp.tif';
+		eval(['!gdal_translate -quiet -projwin ' num2str(x0) ' ' num2str(y1) ' ' num2str(x1) ' ' num2str(y0) ' ' geotiff_name ' ' inputname ]);
+
+		%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.x=(x0:(x1-x0)/(size(md.radaroverlay.pwr,2)-1):x1);
+		md.radaroverlay.y=(y0:(y1-y0)/(size(md.radaroverlay.pwr,1)-1):y1);
+
+		%Erase image or keep it?
+		if ~getfieldvalue(options,'keep_image',0),
+			system('rm -rf ./temp.tif');
+		end
+	end 
+	%}}}
+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');
 	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
+	overlay_image=getfieldvalue(options,'overlay_image')
+	overlay_xlim=getfieldvalue(options,'overlay_xlim')
+	x0
+	x1
+	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
 
 	%Read and crop file
@@ -170,7 +192,7 @@
 	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
-
-%Was a triangulation requested for the area of the image that is not covered by the mesh?
+end %}}}
+
+%Was a triangulation requested for the area of the image that is not covered by the mesh? %{{{
 if strcmpi(getfieldvalue(options,'outertriangulation','no'),'yes'),
 
@@ -184,7 +206,7 @@
 
 	%inner hole from mesh segments: 
-	box(2).x=[md.mesh.x(md.mesh.segments(:,1)); md.mesh.x(md.mesh.segments(end,2))];
+	box(2).x=md.mesh.x(md.mesh.segments(:,1));
 	box(2).x=[box(2).x; box(2).x(1)];
-	box(2).y=[md.mesh.y(md.mesh.segments(:,1)); md.mesh.y(md.mesh.segments(end,2))];
+	box(2).y=md.mesh.y(md.mesh.segments(:,1)); 
 	box(2).y=[box(2).y; box(2).y(1)];
 
@@ -200,4 +222,5 @@
 	maxarea=max(GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y));
 	outermd=triangle(model(),'./outertriangulation.exp',sqrt(maxarea));
+	%outermd=bamg(model(),'domain','./outertriangulation.exp','hmin',sqrt(maxarea));
 	
 	%save the triangulation: 
@@ -206,3 +229,3 @@
 	md.radaroverlay.outery=outermd.mesh.y;
 
-end
+end %}}}