Changeset 20487

Timestamp:

04/09/16 13:58:08 (9 years ago)

Author:

Eric.Larour

Message:

CHG: added options when geotiffname is not for Antarctica or Greenland.

File:

• issm/trunk-jpl/src/m/plot/radarpower.m (modified) (6 diffs)

issm/trunk-jpl/src/m/plot/radarpower.m ¶

@@ -34 +34 @@
 y0=min(ylim); y1=max(ylim);
 
-%figure out if we should go look for Greenland or Antarctica geotiff, or if user provided one.
-if ~exist(options,'overlay_image'),
+if ~exist(options,'overlay_image'), % no image provided, go look into ModelData for one!{{{
         if md.mesh.epsg==3413,   %Greenland 
                 %old code {{{ 
@@ -143 +142 @@
 
         else
-                error('EPSG code not supported yet (ex: 3413 for UPS Greenland, 3031 for UPS Antarctica)');
-        end
-else
+                if ~exist(options,'geotiff_name'), 
+                        error('Need a geotiff for areas outside of Greenland and Antarctica');
+                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:
+                im=imread('temp.tif','TIFF');
+                %adjust contrast and brightness
+                %im=imadjust(im,[a b],[c d]);
+                pixelskip=max(1,ceil(posting/((x1-x0)/(size(im,2)))));
+                md.radaroverlay.pwr=double(flipud(im(1:pixelskip:end,1:pixelskip:end)));
+                md.radaroverlay.x=(x0:(x1-x0)/(size(md.radaroverlay.pwr,2)-1):x1);
+                md.radaroverlay.y=(y0:(y1-y0)/(size(md.radaroverlay.pwr,1)-1):y1);
+
+                %Erase image or keep it?
+                if ~getfieldvalue(options,'keep_image',0),
+                        system('rm -rf ./temp.tif');
+                end
+        end 
+        %}}}
+else %user provided image {{{
         %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
+        overlay_image=getfieldvalue(options,'overlay_image')
+        overlay_xlim=getfieldvalue(options,'overlay_xlim')
+        x0
+        x1
+        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
@@ -170 +192 @@
         md.radaroverlay.x=(x0:(x1-x0)/(size(md.radaroverlay.pwr,2)-1):x1);
         md.radaroverlay.y=(y0:(y1-y0)/(size(md.radaroverlay.pwr,1)-1):y1);
-end
-
-%Was a triangulation requested for the area of the image that is not covered by the mesh?
+end %}}}
+
+%Was a triangulation requested for the area of the image that is not covered by the mesh? %{{{
 if strcmpi(getfieldvalue(options,'outertriangulation','no'),'yes'),
 
@@ -184 +206 @@
 
         %inner hole from mesh segments: 
-        box(2).x=[md.mesh.x(md.mesh.segments(:,1)); md.mesh.x(md.mesh.segments(end,2))];
+        box(2).x=md.mesh.x(md.mesh.segments(:,1));
         box(2).x=[box(2).x; box(2).x(1)];
-        box(2).y=[md.mesh.y(md.mesh.segments(:,1)); md.mesh.y(md.mesh.segments(end,2))];
+        box(2).y=md.mesh.y(md.mesh.segments(:,1)); 
         box(2).y=[box(2).y; box(2).y(1)];
 
@@ -200 +222 @@
         maxarea=max(GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y));
         outermd=triangle(model(),'./outertriangulation.exp',sqrt(maxarea));
+        %outermd=bamg(model(),'domain','./outertriangulation.exp','hmin',sqrt(maxarea));
         
         %save the triangulation: 
@@ -206 +229 @@
         md.radaroverlay.outery=outermd.mesh.y;
 
-end
+end %}}}