Index: /issm/trunk-jpl/src/m/plot/googlemaps.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/googlemaps.m	(revision 15159)
+++ /issm/trunk-jpl/src/m/plot/googlemaps.m	(revision 15160)
@@ -48,4 +48,14 @@
 end
 
+%Get region specific projection parameters
+EPSGgoogle = 'EPSG:3785';   % Mercator       http://www.spatialreference.org/ref/epsg/3785/
+if strcmpi(md.mesh.hemisphere,'n'),
+	EPSGlocal = 'EPSG:3413'; % UPS Greenland  http://www.spatialreference.org/ref/epsg/3413/
+elseif strcmpi(md.mesh.hemisphere,'s'),
+	EPSGlocal = 'EPSG:3031'; % UPS Antarctica http://www.spatialreference.org/ref/epsg/3031/
+else
+	error('field hemisphere should either be ''n'' or ''s''');
+end
+
 %Find optimal zoom
 if exist(options,'zoom'),
@@ -85,6 +95,5 @@
 		position = [num2str(latn) ',' num2str(lonn)];
 		disp(['Google Earth tile: ' num2str(x) '/' num2str(cols-1) ' ' num2str(y) '/' num2str(rows-1) ' (center: ' position ')']);
-		%Google maps API
-		%http://developers.google.com/maps/documentation/staticmaps/
+		%Google maps API: http://developers.google.com/maps/documentation/staticmaps/
 		params = [...
 			'center=' position ...
@@ -105,18 +114,34 @@
 end
 
-%Create coordinates grids
-[gX gY]=meshgrid(ulx:ulx+size(final,2)-1,uly:-1:uly-size(final,1)+1);
-[LAT LON]=pixelstolatlon(gX,gY, zoom);
-if strcmpi(md.mesh.hemisphere,'n'),
-	[X Y]=ll2xy(LAT,LON,+1,45,70);
-elseif strcmpi(md.mesh.hemisphere,'s'),
-	[X Y]=ll2xy(LAT,LON,-1,0,71);
-else
-	error('field hemisphere should either be ''n'' or ''s''');
-end
+%Write image, create geotiff and reproject from mercator to local projection
+imwrite(final,'temp.png','png')
+[ulmx ulmy]=ll2mercator(ullat,ullon);
+[lrmx lrmy]=ll2mercator(lrlat,lrlon);
+
+%Create Geotiff for Mercator projection 
+system(['gdal_translate -of Gtiff -co "tfw=yes"  -a_ullr '...
+	num2str(ulmx,'%15.8f') ' ' num2str(ulmy,'%15.8f') ' ' num2str(lrmx,'%15.8f') ' ' num2str(lrmy,'%15.8f')...
+	' -a_srs "' EPSGgoogle '" "temp.png" "temp.tiff"']);
+delete('temp.png');
+
+%reproject from mercator (EPSG:3785) to UPS Ant (EPSG:3031)
+system(['gdalwarp  -s_srs ' EPSGgoogle ' -t_srs ' EPSGlocal ' temp.tiff temp2.tiff']);
+delete('temp.tiff','temp.tfw');
+
+%Put everything in model
+[status output]=system('gdalinfo temp2.tiff | grep "Upper Left"');
+ul = sscanf(output,'Upper Left  (%f, %f)');
+[status output]=system('gdalinfo temp2.tiff | grep "Lower Right"');
+lr = sscanf(output,'Lower Right (%f, %f)');
+[status output]=system('gdalinfo temp2.tiff | grep "Size is"');
+si = sscanf(output,'Size is %i, %i');
+x_m=linspace(ul(1),lr(1),si(1));
+y_m=linspace(ul(2),lr(2),si(2)); %We need to reverse y_m because the image is read upside down by matlab
+final=imread('temp2.tiff');
+delete('temp2.tiff');
 
 md.radaroverlay.pwr=final;
-md.radaroverlay.x=X;
-md.radaroverlay.y=Y;
+md.radaroverlay.x=x_m;
+md.radaroverlay.y=y_m;
 
 end
Index: /issm/trunk-jpl/src/m/plot/plot_googlemaps.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_googlemaps.m	(revision 15159)
+++ /issm/trunk-jpl/src/m/plot/plot_googlemaps.m	(revision 15160)
@@ -47,7 +47,6 @@
 
 %Retrieve image from md
-X = md.radaroverlay.x;
-Y = md.radaroverlay.y;
-final = md.radaroverlay.pwr;
+[X Y] = meshgrid(md.radaroverlay.x,md.radaroverlay.y);
+final = double(md.radaroverlay.pwr)/double(max(md.radaroverlay.pwr(:))); %rescale between 0 and 1
 
 %Get some options