Index: /issm/trunk-jpl/src/m/plot/landsatmap.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/landsatmap.m	(revision 24478)
+++ /issm/trunk-jpl/src/m/plot/landsatmap.m	(revision 24478)
@@ -0,0 +1,137 @@
+% Explain
+%  upload landsatmap to md.radaroverlay
+%
+% Usage
+%  md = landsatmap(md);
+%  md = landsatmap(md,'highres',1);
+function md = landsatmap(md,varargin),
+
+% check input variables
+if nargin == 1 ,% {{{
+	options = pairoptions;
+else
+	options = pairoptions(varargin{:});
+end% }}}
+
+%process mesh and data
+[x y z elements is2d isplanet]=processmesh(md,[],options);
+
+%check is2d
+if ~is2d,
+   error('buildgridded error message: gridded not supported for 3d meshes, project on a layer');
+end
+
+% get xlim, and ylim
+xlim    = getfieldvalue(options,'xlim',[min(x) max(x)])/getfieldvalue(options,'unit',1);
+ylim    = getfieldvalue(options,'ylim',[min(y) max(y)])/getfieldvalue(options,'unit',1);
+highres = getfieldvalue(options,'highres',0);
+
+if md.mesh.epsg == 3031 % Antarctica region {{{
+	if highres, % high resolution geotiff file
+		if 1, disp('   LIMA with geotiff'), % {{{
+			disp('WARNING : this image shoud be collected with geocoded tif file');
+			% find merged mosaic landsat image {{{
+			limapath = {'/drive/project_inwoo/issm/Data/LIMA/AntarcticaLandsat.tif';
+				'/home/DATA/ICESHEET/LIMA/AntarcticaLandsat.tif'};
+			pos = zeros(length(limapath),1);
+			for ii = 1:length(limapath)
+				if exist(limapath{ii}), pos(ii) = 1; end
+			end
+			limapath = limapath{find(pos)};
+			fprintf('   LIMA path is %s\n', limapath);
+			% }}}
+
+			% read image
+			im = imread(limapath);
+
+			% Region of LIMA data set
+			info = gdalinfo(limapath); % get geotiff info
+			xm = info.xmin + info.dx*[0:info.nx-1];
+			ym = info.ymax - info.dy*[0:info.ny-1];
+
+			%disp('   find region of model at LIMA');
+			offset = 1e+4;
+			posx = find((xm > xlim(1)-offset).* (xm < xlim(2)+offset));
+			posy = find((ym > ylim(1)-offset).* (ym < ylim(2)+offset));
+		end % }}}
+	else
+		if 1, disp('   LIMA with jp2'), % {{{
+			% find merged mosaic landsat image {{{
+			limapath = {'/drive/project_inwoo/issm/Data/LIMA/jp2_100pct/00000-20080314-144756363.jp2';
+				'/data/project_inwoo/issm/Data/LIMA/jp2_100pct/00000-20080314-144756363.jp2';
+				'/home/DATA/ICESHEET/LIMA/jp2_100pct/00000-20080314-144756363.jp2'};
+			pos = zeros(length(limapath),1);
+			for ii = 1:length(limapath)
+				if exist(limapath{ii}), pos(ii) = 1; end
+			end
+			
+			if sum(pos) == 0,
+				fprintf('download website : https://lima.usgs.gov/fullcontinent.php\n');
+				error('Landsat image at Antarctic region should be downloaded at above website');
+			end
+			limapath = limapath{find(pos)};
+			fprintf('   LIMA path is %s\n', limapath);
+			% }}}
+
+			% Resolution and coordinates of upper left corner:
+			xres = 240.010503438;
+			yres = -240.000000516;
+			xul = -2668154.98388;
+			yul = 2362214.96998;
+
+			% Arrays of pixel coordinates:
+			xm = [xul:xres:2.813684914643920e+06];
+			ym = [yul:yres:-2.294505040031947e+06];
+
+			% reduction level 3 corresponds to very 2^3 = 8 points 
+			rlevel = 2;
+			xm = xm(1:2^rlevel:end);
+			ym = ym(1:2^rlevel:end);
+			im = imread(limapath,'reductionlevel',rlevel);
+
+			%disp('   find region of model at LIMA');
+			offset = 1e+4;
+			posx = find((xm > xlim(1)-offset).* (xm < xlim(2)+offset));
+			posy = find((ym > ylim(1)-offset).* (ym < ylim(2)+offset));
+		end % }}}
+		if 0, disp('   LIMA with reduced tiff'), % {{{
+			% find merged mosaic landsat image {{{
+			limapath = {'/drive/project_inwoo/issm/Data/LIMA/tiff_90pct/00000-20080319-092059124.tif'};
+			pos = zeros(length(limapath),1);
+			for ii = 1:length(limapath)
+				if exist(limapath{ii}), pos(ii) = 1; end
+			end
+			
+			if sum(pos) == 0,
+				fprintf('download website : https://lima.usgs.gov/fullcontinent.php\n');
+				error('Landsat image at Antarctic region should be downloaded at above website');
+			end
+			limapath = limapath{find(pos)};
+			fprintf('   LIMA path is %s\n', limapath);
+			% }}}
+
+			% read image
+			im = imread(limapath);
+
+			% Region of LIMA data set
+			info = gdalinfo(limapath); % get geotiff info
+			xm = info.xmin + info.dx*[0:info.nx-1];
+			ym = info.ymax - info.dy*[0:info.ny-1];
+
+			%disp('   find region of model at LIMA');
+			offset = 1e+4;
+			posx = find((xm > xlim(1)-offset).* (xm < xlim(2)+offset));
+			posy = find((ym > ylim(1)-offset).* (ym < ylim(2)+offset));
+		end % }}}
+	end
+
+	% update region of radaroverlay
+	md.radaroverlay.x = xm(posx);
+	md.radaroverlay.y = ym(posy);
+	md.radaroverlay.pwr = im(posy, posx,:);
+
+	% }}}
+else
+	error('Check md.mesh.epsg, available LIMA regeion is at Antarctica (EPSG:3031)');
+end
+
Index: /issm/trunk-jpl/src/m/plot/plot_landsat.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_landsat.m	(revision 24478)
+++ /issm/trunk-jpl/src/m/plot/plot_landsat.m	(revision 24478)
@@ -0,0 +1,152 @@
+% Explain
+%  This funtion loads Landsat Image Mosaic Antarctica (LIMA) for background image.
+%
+% Usage
+%  plot_landsat(md,data,options,plotlines,plotcols,i),
+%
+function plot_landsat(md,data,options,plotlines,plotcols,i),
+
+%process mesh and data
+[x2d y2d z2d elements2d is2d isplanet]=processmesh(md,[],options);
+[data datatype]=processdata(md,data,options);
+
+%check is2d
+if ~is2d,
+   error('buildgridded error message: gridded not supported for 3d meshes, project on a layer');
+end
+
+%Get some options
+transparency = getfieldvalue(options,'transparency',.3);
+highres =  getfieldvalue(options,'highres',0);
+
+% get xlim, and ylim
+xlim=getfieldvalue(options,'xlim',[min(x2d) max(x2d)])/getfieldvalue(options,'unit',1);
+ylim=getfieldvalue(options,'ylim',[min(y2d) max(y2d)])/getfieldvalue(options,'unit',1);
+
+if md.mesh.epsg == 3031 & (isempty(md.radaroverlay.pwr) | isempty(md.radaroverlay.x) | isempty(md.radaroverlay.y) | length(size(md.radaroverlay.pwr)) < 3), % Antarctica region {{{
+	if highres, 
+		disp('   LIMA with geotiff'), % {{{
+		disp('WARNING : this image shoud be collected with geocoded tif file');
+		% find merged mosaic landsat image {{{
+		limapath = {'/drive/project_inwoo/issm/Data/LIMA/AntarcticaLandsat.tif'};
+		pos = zeros(length(limapath),1);
+		for ii = 1:length(limapath)
+			if exist(limapath{ii}), pos(ii) = 1; end
+		end
+		limapath = limapath{find(pos)};
+		fprintf('   LIMA path is %s\n', limapath);
+		% }}}
+
+		% read image
+		im = imread(limapath);
+
+		% Region of LIMA data set
+		info = gdalinfo(limapath); % get geotiff info
+		xm = info.xmin + info.dx*[0:info.nx-1];
+		ym = info.ymax - info.dy*[0:info.ny-1];
+
+		% find region of model at LIMA
+		offset = 1e+4;
+		posx = find((xm > xlim(1)-offset).* (xm < xlim(2)+offset));
+		posy = find((ym > ylim(1)-offset).* (ym < ylim(2)+offset));
+		% }}}
+	else
+		disp('   LIMA with reduced tiff'),
+		% find merged mosaic landsat image {{{
+		limapath = {'/drive/project_inwoo/issm/Data/LIMA/tiff_90pct/00000-20080319-092059124.tif'};
+		pos = zeros(length(limapath),1);
+		for ii = 1:length(limapath)
+			if exist(limapath{ii}), pos(ii) = 1; end
+		end
+		
+		if sum(pos) == 0,
+			fprintf('download website : https://lima.usgs.gov/fullcontinent.php\n');
+			error('Landsat image at Antarctic region should be downloaded at above website');
+		end
+		limapath = limapath{find(pos)};
+		fprintf('   LIMA path is %s\n', limapath);
+		% }}}
+
+		% read image
+		im = imread(limapath);
+
+		% Region of LIMA data set
+		info = gdalinfo(limapath); % get geotiff info
+		xm = info.xmin + info.dx*[0:info.nx-1];
+		ym = info.ymax - info.dy*[0:info.ny-1];
+
+		% find region of model at LIMA
+		offset = 1e+4;
+		posx = find((xm > xlim(1)-offset).* (xm < xlim(2)+offset));
+		posy = find((ym > ylim(1)-offset).* (ym < ylim(2)+offset));
+		% }}}
+	end
+
+	% update region of radaroverlay
+	md.radaroverlay.x = xm(posx);
+	md.radaroverlay.y = ym(posy);
+	md.radaroverlay.pwr = im(posy, posx,:);
+	% }}}
+elseif length(size(md.radaroverlay.pwr)) == 3
+	% it already contains LIMA image.
+elseif md.mesh.epsg == 3431 & (isempty(md.radaroverlay.pwr) | isempty(md.radaroverlay.x) | isempty(md.radaroverlay.y) | length(size(md.radaroverlay.pwr)) < 3), % Greenladn region 
+	error('Greenland region is not yet available.');
+else
+	error('Check md.mesh.epsg, available Landsat regeion is at Antarctica (EPSG:3031)');
+end
+
+%Process image from model
+final = double(md.radaroverlay.pwr)/double(max(md.radaroverlay.pwr(:))); %rescale between 0 and 1
+
+%Prepare grid
+if size(md.radaroverlay.x,1)==1 | size(md.radaroverlay.x,2)==1,
+   x_m = md.radaroverlay.x;
+   y_m = md.radaroverlay.y;
+   data_grid=InterpFromMeshToGrid(elements2d,x2d/getfieldvalue(options,'unit',1),y2d/getfieldvalue(options,'unit',1),data,x_m,y_m,NaN);
+   %data_grid=InterpFromMeshToGrid(md.mesh.elements,md.mesh.x/getfieldvalue(options,'unit',1),md.mesh.y/getfieldvalue(options,'unit',1),data,x_m,y_m,NaN);
+else
+   X = md.radaroverlay.x;
+   Y = md.radaroverlay.y;
+   data_grid=InterpFromMeshToMesh2d(elements2d,x2d,y2d,data,X(:),Y(:),'default',NaN); data_grid=reshape(data_grid,size(X));
+   %data_grid=InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,data,X(:),Y(:),'default',NaN); data_grid=reshape(data_grid,size(X));
+   x_m=X(1,:); y_m=Y(:,1);
+end
+
+data_nan=isnan(data_grid);
+if exist(options,'caxis'),
+   caxis_opt=getfieldvalue(options,'caxis');
+   data_grid(find(data_grid<caxis_opt(1)))=caxis_opt(1);
+   data_grid(find(data_grid>caxis_opt(2)))=caxis_opt(2);
+   data_min=caxis_opt(1);
+   data_max=caxis_opt(2);
+else
+   data_min=min(data_grid(:));
+   data_max=max(data_grid(:));
+end
+colorm = getcolormap(options);
+image_rgb = ind2rgb(uint16((data_grid - data_min)*(length(colorm)/(data_max-data_min))),colorm);
+
+alpha=ones(size(data_grid));
+alpha(find(~data_nan))=transparency;
+alpha=repmat(alpha,[1 1 3]);
+
+final=alpha.*final+(1-alpha).*image_rgb;
+
+%Select plot area 
+subplotmodel(plotlines,plotcols,i,options);
+
+h=imagesc(x_m*getfieldvalue(options,'unit',1),y_m*getfieldvalue(options,'unit',1),final);
+
+%last step: mesh gridded?
+if exist(options,'edgecolor'),
+	A=elements(:,1); B=elements(:,2); C=elements(:,3); 
+	patch('Faces',[A B C],'Vertices', [x y z],'FaceVertexCData',data_grid(1)*ones(size(x)),'FaceColor','none','EdgeColor',getfieldvalue(options,'edgecolor'));
+end
+
+%Apply options
+if ~isnan(data_min),
+	options=changefieldvalue(options,'caxis',[data_min data_max]); % force caxis so that the colorbar is ready
+end
+options=addfielddefault(options,'axis','xy equal off'); % default axis
+applyoptions(md,data,options);
+end
Index: /issm/trunk-jpl/src/m/plot/plot_manager.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_manager.m	(revision 24477)
+++ /issm/trunk-jpl/src/m/plot/plot_manager.m	(revision 24478)
@@ -166,4 +166,10 @@
 
 %Figure out if this is a semi-transparent plot.
+if getfieldvalue(options,'landsat',0),
+	plot_landsat(md,data,options,nlines,ncols,i);
+	return;
+end
+
+%Figure out if this is a semi-transparent plot.
 if exist(options,'gridded'),
 	plot_gridded(md,data,options,nlines,ncols,i);