Index: /issm/trunk/src/m/model/radarpower.m
===================================================================
--- /issm/trunk/src/m/model/radarpower.m	(revision 6669)
+++ /issm/trunk/src/m/model/radarpower.m	(revision 6670)
@@ -1,3 +1,3 @@
-function md=radarpower(md,xlim,ylim,highres)
+function md=radarpower(md,options)
 %RADARPOWER - overlay a power radar image on an existing mesh
 %
@@ -7,6 +7,5 @@
 %
 %   Usage:
-%      md=radarpower(md,xlim,ylim,highres)
-%      md=radarpower(md,xlim,ylim)
+%      md=radarpower(md,options);
 %      md=radarpower(md)
 
@@ -18,11 +17,11 @@
 
 %Parse inputs
-if nargin<4,
-	highres=0;
+if nargin==1,
+	options=pairoptions;
 end
-if nargin<3,
-	xlim=[min(md.x) max(md.x)];
-	ylim=[min(md.y) max(md.y)];
-end
+
+highres=getfieldvalue(options,'highres',0);
+xlim=getfieldvalue(options,'xlim',[min(md.x) max(md.x)]);
+ylim=getfieldvalue(options,'ylim',[min(md.y) max(md.y)]);
 
 %find gdal coordinates
@@ -30,50 +29,76 @@
 y0=min(ylim); y1=max(ylim);
 
-%the geotiff image is either 200m or 1km accuracy. 
-if strcmpi(md.hemisphere,'n'),
-	if ~exist([MODELDATA '/MOG/mog150_greenland_map.jpg']),
-		error(['radarpower error message: file ' MODELDATA '/MOG/mog150_greenland_map.jpg not found. Check MODELDATA variable..']);
+%figure out if we should go look for Greenland or Antarctica geotiff, or if user provided one.
+if ~exist(options,'overlay_image'),
+	if strcmpi(md.hemisphere,'n'),
+		if ~exist([MODELDATA '/MOG/mog150_greenland_map.jpg']),
+			error(['radarpower error message: file ' MODELDATA '/MOG/mog150_greenland_map.jpg not found. Check MODELDATA variable..']);
+		end
+		jpgim=[MODELDATA '/MOG/mog150_greenland_map.jpg'];
+		geom=load([MODELDATA '/MOG/mog150_greenland_map.jpgw'],'ascii');
+		sizex=floor((x1-x0)/geom(1)); % x posting
+		sizey=floor((y1-y0)/geom(4)); % y posting
+		topleftx=floor((x0-geom(5))/geom(1)); % x min
+		toplefty=floor((geom(6)-y1)/geom(4)); % y max
+
+		%Read and crop file
+		disp('Warning: expecting coordinates in polar stereographic (Std Latitude: 70ºN Meridian: 45º)');
+		im=imread(jpgim);
+		im=im(toplefty:toplefty+sizey,topleftx:topleftx+sizex);
+		md.sarpwr=double(flipud(im));
+		md.sarxm=(x0:(x1-x0)/(size(md.sarpwr,2)-1):x1);
+		md.sarym=(y0:(y1-y0)/(size(md.sarpwr,1)-1):y1);
+
+	elseif strcmpi(md.hemisphere,'s'),
+		if highres,
+			if ~exist([MODELDATA '/MosaicTiffRsat/amm125m_v2_200m.tif']),
+				error(['radarpower error message: file ' MODELDATA '/MosaicTiffRsat/amm125m_v2_200m.tif not found. Check MODELDATA variable..']);
+			end
+			geotiff_name=[MODELDATA '/MosaicTiffRsat/amm125m_v2_200m.tif'];
+		else
+			if ~exist([MODELDATA '/MosaicTiffRsat/amm125m_v2_1km.tif']),
+				error(['radarpower error message: file ' MODELDATA '/MosaicTiffRsat/amm125m_v2_1km.tif not found. Check MODELDATA variable..']);
+			end
+			geotiff_name=[MODELDATA '/MosaicTiffRsat/amm125m_v2_1km.tif'];
+		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:
+		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');
+
+	else
+		error('field hemisphere should either be ''n'' or ''s''');
 	end
-	jpgim=[MODELDATA '/MOG/mog150_greenland_map.jpg'];
-	geom=load([MODELDATA '/MOG/mog150_greenland_map.jpgw'],'ascii');
-	sizex=floor((x1-x0)/geom(1)); % x posting
-	sizey=floor((y1-y0)/geom(4)); % y posting
-	topleftx=floor((x0-geom(5))/geom(1)); % x min
-	toplefty=floor((geom(6)-y1)/geom(4)); % y max
+else
+	%ok, 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
+
 
 	%Read and crop file
 	disp('Warning: expecting coordinates in polar stereographic (Std Latitude: 70ºN Meridian: 45º)');
-	im=imread(jpgim);
+	im=imread(overlay_image);
 	im=im(toplefty:toplefty+sizey,topleftx:topleftx+sizex);
 	md.sarpwr=double(flipud(im));
 	md.sarxm=(x0:(x1-x0)/(size(md.sarpwr,2)-1):x1);
 	md.sarym=(y0:(y1-y0)/(size(md.sarpwr,1)-1):y1);
-
-elseif strcmpi(md.hemisphere,'s'),
-	if highres,
-		if ~exist([MODELDATA '/MosaicTiffRsat/amm125m_v2_200m.tif']),
-			error(['radarpower error message: file ' MODELDATA '/MosaicTiffRsat/amm125m_v2_200m.tif not found. Check MODELDATA variable..']);
-		end
-		geotiff_name=[MODELDATA '/MosaicTiffRsat/amm125m_v2_200m.tif'];
-	else
-		if ~exist([MODELDATA '/MosaicTiffRsat/amm125m_v2_1km.tif']),
-			error(['radarpower error message: file ' MODELDATA '/MosaicTiffRsat/amm125m_v2_1km.tif not found. Check MODELDATA variable..']);
-		end
-		geotiff_name=[MODELDATA '/MosaicTiffRsat/amm125m_v2_1km.tif'];
-	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:
-	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');
-
-else
-	error('field hemisphere should either be ''n'' or ''s''');
 end