Index: /issm/trunk-jpl/src/m/plot/plot_googleearth.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_googleearth.m	(revision 15109)
+++ /issm/trunk-jpl/src/m/plot/plot_googleearth.m	(revision 15109)
@@ -0,0 +1,185 @@
+function plot_googleearth(md,data,options,plotlines,plotcols,i)
+%PLOT_OVERLAY - superimpose radar image to a given field
+%
+%   Usage:
+%      plot_googleearth(md,data,options,plotlines,plotcols,i)
+%
+%   See also: PLOTMODEL
+
+%process mesh and data
+[x y z elements 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 xlim and ylim (used to extract radar image)
+xlim=getfieldvalue(options,'xlim',[min(x) max(x)]);
+ylim=getfieldvalue(options,'ylim',[min(y) max(y)]);
+post=getfieldvalue(options,'posting',diff(xlim)/1000);
+
+%Interpolating data on grid
+[x_m y_m data_grid]=InterpFromMeshToGrid(elements,x,y,data,xlim(1),ylim(2),post,post,round(diff(ylim)/post),round(diff(xlim)/post),NaN);
+if size(data_grid,1)<3 | size(data_grid,2)<3,
+	error('data_grid size too small in plot_gridded, check posting and units');
+end
+
+%limits in lat/long
+[ullat,ullon] = xy2ll(xlim(1),ylim(2),+1,45,70);
+[lrlat,lrlon] = xy2ll(xlim(2),ylim(1),+1,45,70);
+
+%Find optimal zoom
+zoom = optimalzoom(ullat,ullon,lrlat,lrlon);
+scale   = 1;
+maxsize = 640;
+
+%convert all these coordinates to pixels
+[ulx, uly]= latlontopixels(ullat, ullon, zoom);
+[lrx, lry]= latlontopixels(lrlat, lrlon, zoom);
+
+%calculate total pixel dimensions of final image
+dx = lrx - ulx;
+dy = uly - lry;
+
+%calculate rows and columns
+cols = ceil(dx/maxsize);
+rows = ceil(dy/maxsize);
+
+%calculate pixel dimensions of each small image
+bottom = 120;
+largura = ceil(dx/cols);
+altura  = ceil(dy/rows);
+alturaplus = altura + bottom;
+
+%Initialize final image
+final = zeros(floor(dy),floor(dx),3);%RGB image
+for x=0:cols-1,
+	for y=0:rows-1,
+		dxn = largura * (0.5 + x);
+		dyn = altura * (0.5 + y);
+		[latn, lonn] = pixelstolatlon(ulx + dxn, uly - dyn - bottom/2, zoom);
+		position = [num2str(latn) ',' num2str(lonn)];
+		disp([ num2str(x) ' ' num2str(y) ' (' position ')']);
+		params = [...
+			'center=' position ...
+			'&zoom=' num2str(zoom)...
+			'&size=' num2str(largura) 'x' num2str(alturaplus)...
+			'&maptype=satellite'...
+			'&sensor=false'...
+			'&scale=' num2str(scale)];
+		url = ['http://maps.google.com/maps/api/staticmap?' params];
+		[X, map]=imread(url,'png');
+		X=ind2rgb(X,map);
+		indx1 = floor(x*largura)+1;
+		indx2 = min(floor(dx),floor(x*largura)+size(X,2));
+		indy1 = floor(y*altura)+1;
+		indy2 = min(floor(dy),floor(y*altura)+size(X,1));
+		final(indy1:indy2,indx1:indx2,:)=X(1:indy2-indy1+1,1:indx2-indx1+1,:);
+	end
+end
+
+%Create model image
+[gX gY]=meshgrid(ulx:ulx+size(final,2)-1,uly:-1:uly-size(final,1)+1);
+[LAT LON]=pixelstolatlon(gX,gY, zoom);
+[X Y]=ll2xy(LAT,LON,+1,45,70);
+data_grid=InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,data,X(:),Y(:),'default',NaN); data_grid=reshape(data_grid,size(X));
+
+%Process data_grid: add white in NaN and correct caxis accordingly
+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(X));
+alpha(find(~data_nan))=0.5;
+alpha=repmat(alpha,[1 1 3]);
+
+final=alpha.*final+(1-alpha).*image_rgb;
+
+%Select plot area 
+subplotmodel(plotlines,plotcols,i,options);
+
+h=imagesc(X(1,:),Y(:,1),final);
+options=addfielddefault(options,'axis','xy equal off'); % default axis
+
+%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
+applyoptions(md,data,options);
+end
+
+function [mx my]=latlontomercator(lat, lon),
+
+	EARTH_RADIUS = 6378137;
+	EQUATOR_CIRCUMFERENCE = 2 * pi * EARTH_RADIUS;
+	INITIAL_RESOLUTION = EQUATOR_CIRCUMFERENCE / 256.0;
+	ORIGIN_SHIFT = EQUATOR_CIRCUMFERENCE / 2.0;
+
+	mx = (lon * ORIGIN_SHIFT) / 180.0;
+	my = log(tan((90 + lat) * pi/360.0))/(pi/180.0);
+	my = (my * ORIGIN_SHIFT) /180.0;
+end
+function [px py]=latlontopixels(lat, lon, zoom),
+
+	EARTH_RADIUS = 6378137;
+	EQUATOR_CIRCUMFERENCE = 2 * pi * EARTH_RADIUS;
+	INITIAL_RESOLUTION = EQUATOR_CIRCUMFERENCE / 256.0;
+	ORIGIN_SHIFT = EQUATOR_CIRCUMFERENCE / 2.0;
+
+    mx = (lon * ORIGIN_SHIFT) / 180.0;
+    my = log(tan((90 + lat) * pi/360.0))/(pi/180.0);
+    my = (my * ORIGIN_SHIFT) /180.0;
+    res = INITIAL_RESOLUTION / (2^zoom);
+    px = (mx + ORIGIN_SHIFT) / res;
+    py = (my + ORIGIN_SHIFT) / res;
+ end
+
+function [lat lon]=pixelstolatlon(px, py, zoom),
+
+	EARTH_RADIUS = 6378137;
+	EQUATOR_CIRCUMFERENCE = 2 * pi * EARTH_RADIUS;
+	INITIAL_RESOLUTION = EQUATOR_CIRCUMFERENCE / 256.0;
+	ORIGIN_SHIFT = EQUATOR_CIRCUMFERENCE / 2.0;
+
+    res = INITIAL_RESOLUTION / (2^zoom);
+    mx = px * res - ORIGIN_SHIFT;
+    my = py * res - ORIGIN_SHIFT;
+    lat = (my / ORIGIN_SHIFT) * 180.0;
+    lat = 180 / pi * (2*atan(exp(lat*pi/180.0)) - pi/2.0);
+    lon = (mx / ORIGIN_SHIFT) * 180.0;
+ end
+function  zoom = optimalzoom(ullat,ullon,lrlat,lrlon)
+
+	EARTH_RADIUS = 6378137;
+	EQUATOR_CIRCUMFERENCE = 2 * pi * EARTH_RADIUS;
+	INITIAL_RESOLUTION = EQUATOR_CIRCUMFERENCE / 256.0;
+
+	optimalsize = 1000; %Number of pixels in final image
+
+	[ulmx ulmy]=latlontomercator(ullat,ullon);
+	[lrmx lrmy]=latlontomercator(lrlat,lrlon);
+	distance = sqrt((lrmx-ulmx)^2 + (lrmy-ulmy)^2);
+
+	zoom1 = floor(log(INITIAL_RESOLUTION*optimalsize/(lrmx-ulmx))/log(2));
+	zoom2 = floor(log(INITIAL_RESOLUTION*optimalsize/(ulmy-lrmy))/log(2));
+
+	zoom=max(zoom1,zoom2);
+
+	zoom = min(max(1,zoom),21);
+end
Index: /issm/trunk-jpl/src/m/plot/plot_manager.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_manager.m	(revision 15108)
+++ /issm/trunk-jpl/src/m/plot/plot_manager.m	(revision 15109)
@@ -157,4 +157,10 @@
 
 %Figure out if this is a semi-transparent plot.
+if exist(options,'googleearth'),
+	plot_googleearth(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);