Index: /issm/trunk-jpl/src/m/plot/radarpower.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/radarpower.m	(revision 13126)
+++ /issm/trunk-jpl/src/m/plot/radarpower.m	(revision 13127)
@@ -10,6 +10,4 @@
 %      md=radarpower(md)
 
-%If gdal does not work, uncomment the following line
-%setenv('LD_LIBRARY_PATH','/proj/ice/larour/issm/trunk/externalpackages/gdal/install/lib/');
 %Parse inputs
 if nargin==1,
@@ -22,8 +20,8 @@
 end
 
-highres=getfieldvalue(options,'highres',0);
-xlim=getfieldvalue(options,'xlim',[min(md.mesh.x) max(md.mesh.x)]);
-ylim=getfieldvalue(options,'ylim',[min(md.mesh.y) max(md.mesh.y)]);
-posting=getfieldvalue(options,'posting',0); % 0 -> image posting default
+highres = getfieldvalue(options,'highres',0);
+xlim    = getfieldvalue(options,'xlim',[min(md.mesh.x) max(md.mesh.x)]);
+ylim    = getfieldvalue(options,'ylim',[min(md.mesh.y) max(md.mesh.y)]);
+posting = getfieldvalue(options,'posting',0); % 0 -> image posting default
 
 %find gdal coordinates
@@ -34,32 +32,62 @@
 if ~exist(options,'overlay_image'),
 	if strcmpi(md.mesh.hemisphere,'n'),
-		if ~exist([jplsvn() '/projects/ModelData/MOG/mog150_greenland_map.jpg']),
-			error(['radarpower error message: file ' jplsvn() '/projects/ModelData/MOG/mog150_greenland_map.jpg not found.']);
+		%if ~exist([jplsvn() '/projects/ModelData/MOG/mog150_greenland_map.jpg']),
+		%	error(['radarpower error message: file ' jplsvn() '/projects/ModelData/MOG/mog150_greenland_map.jpg not found.']);
+		%end
+		%name = 'mog150_greenland_map';
+		%name = 'mog100_hp1_v10';
+		%%name = 'mog500_hp1_v10';
+		%jpgim=[jplsvn() '/projects/ModelData/MOG/' name '.jpg'];
+		%geom=load([jplsvn() '/projects/ModelData/MOG/' name '.jpgw'],'ascii');
+		%%jpgim='/u/astrid-r1b/morlighe/issmjpl/projects/MorlighemGRL2012/Data/Mosaic_amp_asar2010.jpg';
+		%%geom=load('/u/astrid-r1b/morlighe/issmjpl/projects/MorlighemGRL2012/Data/Mosaic_amp_asar2010.jpgw');
+		%jpgim='/u/astrid-r1b/morlighe/issmjpl/projects/MorlighemGRL2012/Data/Russel_asar2010.png';
+		%geom=load('/u/astrid-r1b/morlighe/issmjpl/projects/MorlighemGRL2012/Data/Russel_asar2010.pngw');
+
+		%%geom:   xposting nbcols nbrows yposting xmin ymax
+		%xmin=max(geom(5),x0);
+		%xmax=min(geom(5)+geom(1)*geom(2),x1);
+		%ymin=max(geom(6)-geom(3)*geom(4),y0);
+		%ymax=min(geom(6),y1);
+
+		%firstcol=max(1,floor((xmin-geom(5))/geom(1))); %x min
+		%firstrow=max(1,floor((geom(6)-ymax)/geom(4))); %y max
+		%numcols=floor((xmax-xmin)/geom(1)); % x posting
+		%numrows=floor((ymax-ymin)/geom(4)); % y posting
+		%pixelskip=max(1,ceil(posting/geom(1)));
+
+		%%Read and crop file
+		%disp('Warning: expecting coordinates in polar stereographic (Std Latitude: 70ºN Meridian: 45º)');
+		%im=imread(jpgim);
+		%im=im(firstrow:firstrow+numrows-1,firstcol:firstcol+numcols-1);
+		%md.radaroverlay.pwr=double(flipud(im(1:pixelskip:end,1:pixelskip:end)));
+		%md.radaroverlay.x=(xmin:(xmax-xmin)/(size(md.radaroverlay.pwr,2)-1):xmax);
+		%md.radaroverlay.y=(ymin:(ymax-ymin)/(size(md.radaroverlay.pwr,1)-1):ymax);
+		if highres,
+			if ~exist([jplsvn() '/projects/ModelData/MOG/mog100g_r2_hp1.tif']),
+				error(['radarpower error message: file ' jplsvn() '/projects/ModelData/MOG/mog100g_r2_hp1.tif not found.']);
+			end
+			geotiff_name=[jplsvn() '/projects/ModelData/MOG/mog100_r2_hp1.tif'];
+		else
+			if ~exist([jplsvn() '/projects/ModelData/MOG/mog500g_r2_hp1.tif']),
+				error(['radarpower error message: file ' jplsvn() '/projects/ModelData/MOG/mog500g_r2_hp1.tif not found.']);
+			end
+			geotiff_name=[jplsvn() '/projects/ModelData/MOG/mog500_r2_hp1.tif'];
 		end
-		name = 'mog150_greenland_map';
-		%name = 'mog100_hp1_v10';
-		%name = 'mog500_hp1_v10';
-		jpgim=[jplsvn() '/projects/ModelData/MOG/' name '.jpg'];
-		geom=load([jplsvn() '/projects/ModelData/MOG/' name '.jpgw'],'ascii');
 
-		%geom:   xposting nbcols nbrows yposting xmin ymax
-		xmin=max(geom(5),x0);
-		xmax=min(geom(5)+geom(1)*geom(2),x1);
-		ymin=max(geom(6)-geom(3)*geom(4),y0);
-		ymax=min(geom(6),y1);
+		%Name of image
+		inputname='./temp.tif';
+		eval(['!gdal_translate -quiet -projwin ' num2str(x0) ' ' num2str(y1) ' ' num2str(x1) ' ' num2str(y0) ' ' geotiff_name ' ' inputname ]);
 
-		firstcol=max(1,floor((xmin-geom(5))/geom(1))); %x min
-		firstrow=max(1,floor((geom(6)-ymax)/geom(4))); %y max
-		numcols=floor((xmax-xmin)/geom(1)); % x posting
-		numrows=floor((ymax-ymin)/geom(4)); % y posting
-		pixelskip=max(1,ceil(posting/geom(1)));
+		%Read in temp.tif:
+		im=imread('temp.tif','TIFF');
+		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);
 
-		%Read and crop file
-		disp('Warning: expecting coordinates in polar stereographic (Std Latitude: 70ºN Meridian: 45º)');
-		im=imread(jpgim);
-		im=im(firstrow:firstrow+numrows-1,firstcol:firstcol+numcols-1);
-		md.radaroverlay.pwr=double(flipud(im(1:pixelskip:end,1:pixelskip:end)));
-		md.radaroverlay.x=(xmin:(xmax-xmin)/(size(md.radaroverlay.pwr,2)-1):xmax);
-		md.radaroverlay.y=(ymin:(ymax-ymin)/(size(md.radaroverlay.pwr,1)-1):ymax);
+		%Erase image
+		system('rm -rf ./temp.tif');
+
 
 	elseif strcmpi(md.mesh.hemisphere,'s'),
@@ -94,5 +122,5 @@
 	end
 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.
+	%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');