Changeset 2725

Timestamp:

12/11/09 14:07:23 (15 years ago)

Author:

Mathieu Morlighem

Message:

Added overlay for Greenland! (youhou)

Location:

issm/trunk/src/m/classes/public

Files:

• plot/plot_overlay.m (modified) (3 diffs)
• plot/plotdoc.m (modified) (1 diff)
• radarpower.m (modified) (3 diffs)

issm/trunk/src/m/classes/public/plot/plot_overlay.m

@@ -31 +31 @@
 if redo,
         t1=clock; fprintf('%s','Extracting radar image...');
-        md=radarpower(md,xlim,ylim,getfieldvalue(options,'highres',0));
+        md=radarpower(md,getfieldvalue(options,'hem','s'),xlim,ylim,getfieldvalue(options,'highres',0));
         t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
 end
@@ -75 +75 @@
 
 %intensity
-v=radar/max(max(radar));
+v=radar/max(radar(:));
 
 %hue
@@ -87 +87 @@
 
 %Include border
-v((1:border),:)=0;  v((end-border+1:end),:)=0; v(:,1:border)=0;v(:,(end-border+1:end))=0;
+%v((1:border),:)=0;  v((end-border+1:end),:)=0; v(:,1:border)=0;v(:,(end-border+1:end))=0;
 
 %Transform hsv to rgb

issm/trunk/src/m/classes/public/plot/plotdoc.m

@@ -85 +85 @@
 disp('       ''overlay'': yes or no. This will overlay a radar amplitude image behind');
 disp('       ''highres'': resolution of overlayed radar amplitude image (default is 0, high resolution is 1).');
+disp('       ''hem'': specify hemisphere ''n'' or ''s'' (default is ''s'').');
 disp('       ''alpha'': transparency coefficient (the higher, the more transparent). Default is 1.5');
 disp('       ''scaling'': scaling factor used by quiver plots. Default is 0.4');

issm/trunk/src/m/classes/public/radarpower.m

@@ -1 @@
-function md=radarpower(md,xlim,ylim,highres)
+function md=radarpower(md,hemisphere,xlim,ylim,highres)
 %RADARPOWER - overlay a power radar image on an existing mesh
 %
@@ -7 +7 @@
 %
 %   Usage:
-%      md=radarpower(md,xlim,ylim,highres)
+%      md=radarpower(md,hemisphere,xlim,ylim,highres)
 
 global ISSM_DIR
@@ -31 +31 @@
 
 %the geotiff image is either 200m or 1km accuracy. 
-if highres,
-        if ~exist([MODELDATA '/MosaicTiffRsat/amm125m_v2_200m.tif']),
-                error(['radarpower error message: file ' MODELDATA '/MosaicTiffRsat/amm125m_v2_200m.tif not found. Check MODELDATA variable..']);
+if strcmpi(hemisphere,'n'),
+        if ~exist([MODELDATA '/MOG/mog150_greenland_map.jpg']),
+                error(['radarpower error message: file ' MODELDATA '/MOG/mog150_greenland_map.jpg not found. Check MODELDATA variable..']);
         end
-        geotiff_name=[MODELDATA '/MosaicTiffRsat/amm125m_v2_200m.tif'];
+        jpgim=[MODELDATA '/MOG/mog150_greenland_map.jpg'];
+        geom=load([MODELDATA '/MOG/mog150_greenland_map.jpgw'],'ascii');
+        sizex=floor((x1-x0)/geom(1)); % x posting
+        sizey=floor((y1-y0)/geom(4)); % y posting
+        topleftx=floor((x0-geom(5))/geom(1)); % x min
+        toplefty=floor((geom(6)-y1)/geom(4)); % y max
+
+        %Read and crop file
+        tic
+        im=imread(jpgim);
+        toc
+        tic
+        im=im(toplefty:toplefty+sizey,topleftx:topleftx+sizex);
+        toc
+        md.sarpwr=double(flipud(im));
+        md.sarxm=(x0:(x1-x0)/(size(md.sarpwr,2)-1):x1);
+        md.sarym=(y0:(y1-y0)/(size(md.sarpwr,1)-1):y1);
+
 else
-        if ~exist([MODELDATA '/MosaicTiffRsat/amm125m_v2_1km.tif']),
-                error(['radarpower error message: file ' MODELDATA '/MosaicTiffRsat/amm125m_v2_1km.tif not found. Check MODELDATA variable..']);
+        if highres,
+                if ~exist([MODELDATA '/MosaicTiffRsat/amm125m_v2_200m.tif']),
+                        error(['radarpower error message: file ' MODELDATA '/MosaicTiffRsat/amm125m_v2_200m.tif not found. Check MODELDATA variable..']);
+                end
+                geotiff_name=[MODELDATA '/MosaicTiffRsat/amm125m_v2_200m.tif'];
+        else
+                if ~exist([MODELDATA '/MosaicTiffRsat/amm125m_v2_1km.tif']),
+                        error(['radarpower error message: file ' MODELDATA '/MosaicTiffRsat/amm125m_v2_1km.tif not found. Check MODELDATA variable..']);
+                end
+                geotiff_name=[MODELDATA '/MosaicTiffRsat/amm125m_v2_1km.tif'];
         end
-        geotiff_name=[MODELDATA '/MosaicTiffRsat/amm125m_v2_1km.tif'];
+
+        %Name of image
+        size(imread(geotiff_name,'TIFF'))
+        inputname='./temp.tif';
+        system([path_gdal 'gdal_translate -quiet -projwin ' num2str(x0) ' ' num2str(y1) ' ' num2str(x1) ' ' num2str(y0) ' ' geotiff_name ' ' inputname ]);
+
+        %Read in temp.tif:
+        md.sarpwr=double(flipud(imread('temp.tif','TIFF')));
+        md.sarxm=(x0:(x1-x0)/(size(md.sarpwr,2)-1):x1);
+        md.sarym=(y0:(y1-y0)/(size(md.sarpwr,1)-1):y1);
+
+        %Erase image
+        system('rm -rf ./temp.tif');
+
 end
-
-%Name of image
-inputname='./temp.tif';
-
-command=['system(''' path_gdal 'gdal_translate -quiet -projwin ' num2str(x0) ' ' num2str(y1) ' ' num2str(x1) ' ' num2str(y0) ' ' geotiff_name ' ' inputname ''');'];
-eval(command);
-
-%Read in temp.tif:
-md.sarpwr=double(flipud(imread('temp.tif','TIFF')));
-md.sarxm=(x0:(x1-x0)/(size(md.sarpwr,2)-1):x1);
-md.sarym=(y0:(y1-y0)/(size(md.sarpwr,1)-1):y1);
-
-%Erase image
-system('rm -rf ./temp.tif');