Changeset 27606

Timestamp:

02/22/23 07:30:10 (2 years ago)

Author:

Mathieu Morlighem

Message:

CHG: enable MOG or other geotif images

Location:

issm/trunk-jpl/src/m/plot

Files:

• plot_overlay.m (modified) (4 diffs)
• radarpower.m (modified) (1 diff)

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

@@ -42 +42 @@
         xlim=getfieldvalue(options,'xlim',[min(x) max(x)])/getfieldvalue(options,'unit',1);
         ylim=getfieldvalue(options,'ylim',[min(y) max(y)])/getfieldvalue(options,'unit',1);
+        if exist(options, 'axis');
+                myaxis = getfieldvalue(options,'axis');
+                xlim = [myaxis(1), myaxis(2)];
+                ylim = [myaxis(3), myaxis(4)];
+        end
         options=addfielddefault(options,'xlim',xlim);
         options=addfielddefault(options,'ylim',ylim);
@@ -48 +53 @@
 contrast = getfieldvalue(options,'contrast',1);  
 radar    = md.radaroverlay.pwr;
-if size(radar,3)>1,
-        disp('WARNING: color image converted to greyscale intensity image');
-        if strcmp(class(radar),'uint8'),
-                radar=double(sum(radar,3))/(255*3);
-        else
-                radar=sum(radar,3)/3;
+
+if ~radaronly
+        if size(radar,3)>1,
+                disp('WARNING: color image converted to greyscale intensity image');
+                if strcmp(class(radar),'uint8'),
+                        radar=double(sum(radar,3))/(255*3);
+                else
+                        radar=sum(radar,3)/3;
+                end
         end
+        if getfieldvalue(options,'backgroundbtw',0)
+                radar(find(radar==0))=1; %Change background from black to white
+        end
+        radar = radar.^(contrast);
+        radar = radar./max(radar(:));
 end
-if getfieldvalue(options,'backgroundbtw',0)
-        radar(find(radar==0))=1; %Change background from black to white
-end
-radar = radar.^(contrast);
-radar = radar./max(radar(:));
 
 if getfieldvalue(options,'backgroundbtw',0)
@@ -101 +109 @@
 
 %Special colormaps that require hsv treatment
-colorm=getfieldvalue(options,'colormap','Rignot');
+colorm=getfieldvalue(options,'colormap','parula');
 if strcmpi(colorm,'Rignot') | strcmpi(colorm,'Seroussi') | strcmpi(colorm,'redblue')
         if strcmpi(colorm,'Rignot'),
@@ -133 +141 @@
         image_rgb=hsv2rgb(image_hsv);
 else
-        colorm = getcolormap(options);
-        len    = size(colorm,1);
+        if radaronly
+                image_rgb = radar;
+        else
+                colorm = getcolormap(options);
+                len    = size(colorm,1);
 
-        ind = ceil((len-1)*(data_grid-data_min)/(data_max - data_min + eps) +1);
-        ind(find(ind>len))=len;
-        image_rgb=zeros(size(data_grid,1),size(data_grid,2),3);
-        r=colorm(:,1); image_rgb(:,:,1)=r(ind); clear r;
-        g=colorm(:,2); image_rgb(:,:,2)=g(ind); clear g;
-        b=colorm(:,3); image_rgb(:,:,3)=b(ind); clear b;
+                ind = ceil((len-1)*(data_grid-data_min)/(data_max - data_min + eps) +1);
+                ind(find(ind>len))=len;
+                image_rgb=zeros(size(data_grid,1),size(data_grid,2),3);
+                r=colorm(:,1); image_rgb(:,:,1)=r(ind); clear r;
+                g=colorm(:,2); image_rgb(:,:,2)=g(ind); clear g;
+                b=colorm(:,3); image_rgb(:,:,3)=b(ind); clear b;
 
-        %Now add radarmap
-        r = image_rgb(:,:,1).*radar;  r(data_nan) = radar(data_nan);  image_rgb(:,:,1) = r;  clear r;
-        g = image_rgb(:,:,2).*radar;  g(data_nan) = radar(data_nan);  image_rgb(:,:,2) = g;  clear g;
-        b = image_rgb(:,:,3).*radar;  b(data_nan) = radar(data_nan);  image_rgb(:,:,3) = b;  clear b;
+                %Now add radarmap
+                r = image_rgb(:,:,1).*radar;  r(data_nan) = radar(data_nan);  image_rgb(:,:,1) = r;  clear r;
+                g = image_rgb(:,:,2).*radar;  g(data_nan) = radar(data_nan);  image_rgb(:,:,2) = g;  clear g;
+                b = image_rgb(:,:,3).*radar;  b(data_nan) = radar(data_nan);  image_rgb(:,:,3) = b;  clear b;
+        end
 end
 

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

@@ -100 +100 @@
         %}}}
 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');
+        filename = getfieldvalue(options,'overlay_image');
+        [filepath,name,ext] = fileparts(filename);
+        if ~exist(filename)
+                error([filename ' not found']);
         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
+        %Is it a geotiff?
+        if strcmp(ext,'.tiff') || strcmp(ext,'.tif')
 
-        %Read and crop file
-        disp('Warning: expecting coordinates in polar stereographic (Std Latitude: 70ºN Meridian: 45º)');
-        im=imread(overlay_image);
-        %adjust contrast and brightness
-        %im=imadjust(im,[a b],[c d]);
-        im=im(toplefty:toplefty+sizey,topleftx:topleftx+sizex);
-        md.radaroverlay.pwr=double(flipud(im));
-        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);
+                %Crop image from xylim
+                tempfilename='./temp.tif';
+                eval(['!gdal_translate -quiet -projwin ' num2str(x0) ' ' num2str(y1) ' ' num2str(x1) ' ' num2str(y0) ' ' filename ' ' tempfilename]);
+
+                %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.pwr=double(im(1:pixelskip:end,1:pixelskip:end,:))/255;
+                md.radaroverlay.x=x0:(x1-x0)/(size(md.radaroverlay.pwr,2)-1):x1;
+                md.radaroverlay.y=y1:-(y1-y0)/(size(md.radaroverlay.pwr,1)-1):y0;
+
+                %Erase image or keep it?
+                if ~getfieldvalue(options,'keep_image',0),
+                        delete(tempfilename);
+                end
+        else
+                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_xlim=getfieldvalue(options,'overlay_xlim');
+                overlay_ylim=getfieldvalue(options,'overlay_ylim');
+                overlay_xposting=getfieldvalue(options,'overlay_xposting');
+                overlay_yposting=getfieldvalue(options,'overlay_yposting');
+                overlay_image=getfieldvalue(options,'overlay_image');
+
+                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(overlay_image);
+                %adjust contrast and brightness
+                %im=imadjust(im,[a b],[c d]);
+                im=im(toplefty:toplefty+sizey,topleftx:topleftx+sizex);
+                md.radaroverlay.pwr=double(flipud(im));
+                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
 end %}}}