Changeset 15160

Timestamp:

05/30/13 10:00:20 (12 years ago)

Author:

Mathieu Morlighem

Message:

CHG: rotate Google map image so that we usethe same projection (requires new gdal and proj.4...)

Location:

issm/trunk-jpl/src/m/plot

Files:

• googlemaps.m (modified) (3 diffs)
• plot_googlemaps.m (modified) (1 diff)

issm/trunk-jpl/src/m/plot/googlemaps.m

@@ -48 +48 @@
 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 +95 @@
                 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 +114 @@
 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

issm/trunk-jpl/src/m/plot/plot_googlemaps.m

@@ -47 +47 @@
 
 %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