Changeset 6670

Timestamp:

11/26/10 12:35:04 (14 years ago)

Author:

Eric.Larour

Message:

More flexible radarpower

File:

• issm/trunk/src/m/model/radarpower.m (modified) (4 diffs)

issm/trunk/src/m/model/radarpower.m

@@ -1 @@
-function md=radarpower(md,xlim,ylim,highres)
+function md=radarpower(md,options)
 %RADARPOWER - overlay a power radar image on an existing mesh
 %
@@ -7 +7 @@
 %
 %   Usage:
-%      md=radarpower(md,xlim,ylim,highres)
-%      md=radarpower(md,xlim,ylim)
+%      md=radarpower(md,options);
 %      md=radarpower(md)
 
@@ -18 +17 @@
 
 %Parse inputs
-if nargin<4,
-        highres=0;
+if nargin==1,
+        options=pairoptions;
 end
-if nargin<3,
-        xlim=[min(md.x) max(md.x)];
-        ylim=[min(md.y) max(md.y)];
-end
+
+highres=getfieldvalue(options,'highres',0);
+xlim=getfieldvalue(options,'xlim',[min(md.x) max(md.x)]);
+ylim=getfieldvalue(options,'ylim',[min(md.y) max(md.y)]);
 
 %find gdal coordinates
@@ -30 +29 @@
 y0=min(ylim); y1=max(ylim);
 
-%the geotiff image is either 200m or 1km accuracy. 
-if strcmpi(md.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..']);
+%figure out if we should go look for Greenland or Antarctica geotiff, or if user provided one.
+if ~exist(options,'overlay_image'),
+        if strcmpi(md.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
+                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
+                disp('Warning: expecting coordinates in polar stereographic (Std Latitude: 70ºN Meridian: 45º)');
+                im=imread(jpgim);
+                im=im(toplefty:toplefty+sizey,topleftx:topleftx+sizex);
+                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);
+
+        elseif strcmpi(md.hemisphere,'s'),
+                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
+
+                %Name of image
+                inputname='./temp.tif';
+                eval(['!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');
+
+        else
+                error('field hemisphere should either be ''n'' or ''s''');
         end
-        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
+else
+        %ok, user provided an image. check we also have overlay_xlim and overlay_ylim  options, to know what range of coordinates the image covers.
+        if (~exist(options,'overlay_xlim') | ~exist(options,'overlay_xlim')| ~exist(options,'overlay_xposting')| ~exist(options,'overlay_yposting')),
+                error('radarpower error message: please provide overlay_xlim, overlay_ylim, overlay_xposting and overlay_yposting options together with overlay_image option');
+        end
+        overlay_image=getfieldvalue(options,'overlay_image');
+        overlay_xlim=getfieldvalue(options,'overlay_xlim');
+        overlay_ylim=getfieldvalue(options,'overlay_ylim');
+        overlay_xposting=getfieldvalue(options,'overlay_xposting');
+        overlay_yposting=getfieldvalue(options,'overlay_yposting');
+
+        sizex=floor((x1-x0)/overlay_xposting);
+        sizey=floor((y1-y0)/overlay_yposting);
+        topleftx=floor((x0-overlay_xlim(1))/overlay_xposting); % x min
+        toplefty=floor((overlay_ylim(2)-y1)/overlay_yposting); % y max
+
 
         %Read and crop file
         disp('Warning: expecting coordinates in polar stereographic (Std Latitude: 70ºN Meridian: 45º)');
-        im=imread(jpgim);
+        im=imread(overlay_image);
         im=im(toplefty:toplefty+sizey,topleftx:topleftx+sizex);
         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);
-
-elseif strcmpi(md.hemisphere,'s'),
-        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
-
-        %Name of image
-        inputname='./temp.tif';
-        eval(['!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');
-
-else
-        error('field hemisphere should either be ''n'' or ''s''');
 end