Changeset 21219


Ignore:
Timestamp:
09/22/16 09:30:07 (9 years ago)
Author:
Mathieu Morlighem
Message:

CHG: simplified radarpower and added list of multiple possible paths for geotiff images

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r20776 r21219  
    2020end
    2121
    22 geotiff_name = getfieldvalue(options,'geotiff_name','');
    2322highres = getfieldvalue(options,'highres',0);
    2423xlim    = getfieldvalue(options,'xlim',[min(md.mesh.x) max(md.mesh.x)]);
    … …  
    3534
    3635if ~exist(options,'overlay_image'), % no image provided, go look into ModelData for one!{{{
    37         if md.mesh.epsg==3413,   %Greenland
    38                 %old code {{{
    39                 %if ~exist(['/u/astrid-r1b/ModelData/MOG/mog150_greenland_map.jpg']),
    40                 %       error(['radarpower error message: file ' '/u/astrid-r1b/ModelData/MOG/mog150_greenland_map.jpg not found.']);
    41                 %end
    42                 %name = 'mog150_greenland_map';
    43                 %name = 'mog100_hp1_v10';
    44                 %%name = 'mog500_hp1_v10';
    45                 %jpgim=['/u/astrid-r1b/ModelData/MOG/' name '.jpg'];
    46                 %geom=load(['/u/astrid-r1b/ModelData/MOG/' name '.jpgw'],'ascii');
    47                 %%jpgim='/u/astrid-r1b/morlighe/issmjpl/projects/MorlighemGRL2012/Data/Mosaic_amp_asar2010.jpg';
    48                 %%geom=load('/u/astrid-r1b/morlighe/issmjpl/projects/MorlighemGRL2012/Data/Mosaic_amp_asar2010.jpgw');
    49                 %jpgim='/u/astrid-r1b/morlighe/issmjpl/projects/MorlighemGRL2012/Data/Russel_asar2010.png';
    50                 %geom=load('/u/astrid-r1b/morlighe/issmjpl/projects/MorlighemGRL2012/Data/Russel_asar2010.pngw');
     36        if exist(options,'geotiff_name'),
     37                paths = {getfieldvalue(options,'geotiff_name')};
     38        elseif md.mesh.epsg==3031, %Antarctica
     39                        if highres,
     40                                paths = {'/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif','/home/ModelData/Antarctica/MosaicTiffRsat/amm125m_v2_200m.tif'};
     41                        else
     42                                paths = {'/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif','/home/ModelData/Antarctica/MosaicTiffRsat/amm125m_v2_1km.tif'};
     43                        end
     44        elseif md.mesh.epsg==3413,   %Greenland
     45                if highres,
     46                        paths = {'/u/astrid-r1b/ModelData/MOG/mog100_r2_hp1.tif','/home/ModelData/Greenland/MOG/mog100_r2_hp1.tif'};
     47                else
     48                        paths = {'/u/astrid-r1b/ModelData/MOG/mog500_r2_hp1.tif','/home/ModelData/Greenland/MOG/mog500_r2_hp1.tif'};
     49                end
     50        else
     51                error('Need to provide geotiff for areas outside of Greenland and Antarctica');
     52        end
    5153
    52                 %%geom:   xposting nbcols nbrows yposting xmin ymax
    53                 %xmin=max(geom(5),x0);
    54                 %xmax=min(geom(5)+geom(1)*geom(2),x1);
    55                 %ymin=max(geom(6)-geom(3)*geom(4),y0);
    56                 %ymax=min(geom(6),y1);
    57 
    58                 %firstcol=max(1,floor((xmin-geom(5))/geom(1))); %x min
    59                 %firstrow=max(1,floor((geom(6)-ymax)/geom(4))); %y max
    60                 %numcols=floor((xmax-xmin)/geom(1)); % x posting
    61                 %numrows=floor((ymax-ymin)/geom(4)); % y posting
    62                 %pixelskip=max(1,ceil(posting/geom(1)));
    63 
    64                 %%Read and crop file
    65                 %disp('Warning: expecting coordinates in polar stereographic (Std Latitude: 70ºN Meridian: 45º)');
    66                 %im=imread(jpgim);
    67                 %im=im(firstrow:firstrow+numrows-1,firstcol:firstcol+numcols-1);
    68                 %md.radaroverlay.pwr=double(flipud(im(1:pixelskip:end,1:pixelskip:end)));
    69                 %md.radaroverlay.x=(xmin:(xmax-xmin)/(size(md.radaroverlay.pwr,2)-1):xmax);
    70                 %md.radaroverlay.y=(ymin:(ymax-ymin)/(size(md.radaroverlay.pwr,1)-1):ymax); % }}}
    71                 if ~exist(options,'geotiff_name'),
    72                         if highres,
    73                                 if ~exist(['/u/astrid-r1b/ModelData/MOG/mog100_r2_hp1.tif']),
    74                                         error(['radarpower error message: file ' '/u/astrid-r1b/ModelData/MOG/mog100_r2_hp1.tif not found.']);
    75                                 end
    76                                 geotiff_name=['/u/astrid-r1b/ModelData/MOG/mog100_r2_hp1.tif'];
    77                         else
    78                                 if ~exist(['/u/astrid-r1b/ModelData/MOG/mog500_r2_hp1.tif']),
    79                                         error(['radarpower error message: file ' '/u/astrid-r1b/ModelData/MOG/mog500_r2_hp1.tif not found.']);
    80                                 end
    81                                 geotiff_name=['/u/astrid-r1b/ModelData/MOG/mog500_r2_hp1.tif'];
    82                         end
     54        %Find file from list
     55        found = false;
     56        for i=1:numel(paths),
     57                if exist(paths{i},'file'),
     58                        geotiff_name = paths{i}; found = true;
    8359                end
    84 
    85                 %Name of image
    86                 inputname='./temp.tif';
    87                 eval(['!gdal_translate -quiet -projwin ' num2str(x0) ' ' num2str(y1) ' ' num2str(x1) ' ' num2str(y0) ' ' geotiff_name ' ' inputname ]);
    88 
    89                 %Read in temp.tif:
    90                 im=imread('temp.tif','TIFF');
    91                 %adjust contrast and brightness
    92                 %im=imadjust(im,[a b],[c d]);
    93                 pixelskip=max(1,ceil(posting/((x1-x0)/(size(im,2)))));
    94                 if size(im,3)==1,
    95                         md.radaroverlay.pwr=double(flipud(im(1:pixelskip:end,1:pixelskip:end)));
    96                 else
    97                         md.radaroverlay.pwr=double(im(1:pixelskip:end,1:pixelskip:end,1:pixelskip:end));
    98                         md.radaroverlay.pwr(:,:,1)=flipud(md.radaroverlay.pwr(:,:,1));
    99                         md.radaroverlay.pwr(:,:,2)=flipud(md.radaroverlay.pwr(:,:,2));
    100                         md.radaroverlay.pwr(:,:,3)=flipud(md.radaroverlay.pwr(:,:,3));
    101                 end
    102                 md.radaroverlay.x=(x0:(x1-x0)/(size(md.radaroverlay.pwr,2)-1):x1);
    103                 md.radaroverlay.y=(y0:(y1-y0)/(size(md.radaroverlay.pwr,1)-1):y1);
    104 
    105                 %Erase image or keep it?
    106                 if ~getfieldvalue(options,'keep_image',0),
    107                         system('rm -rf ./temp.tif');
    108                 end
    109         elseif md.mesh.epsg==3031, %Antarctica
    110                 if ~exist(options,'geotiff_name'),
    111                         if highres,
    112                                 if ~exist(['/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif']),
    113                                         error(['radarpower error message: file ' '/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif not found.']);
    114                                 end
    115                                 geotiff_name=['/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif'];
    116                         else
    117                                 if ~exist(['/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif']),
    118                                         error(['radarpower error message: file ' '/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif not found.']);
    119                                 end
    120                                 geotiff_name=['/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif'];
    121                         end
    122                 end
    123 
    124                 %Name of image
    125                 inputname='./temp.tif';
    126                 eval(['!gdal_translate -quiet -projwin ' num2str(x0) ' ' num2str(y1) ' ' num2str(x1) ' ' num2str(y0) ' ' geotiff_name ' ' inputname ]);
    127 
    128                 %Read in temp.tif:
    129                 im=imread('temp.tif','TIFF');
    130                 %adjust contrast and brightness
    131                 %im=imadjust(im,[a b],[c d]);
    132                 pixelskip=max(1,ceil(posting/((x1-x0)/(size(im,2)))));
    133                 md.radaroverlay.pwr=double(flipud(im(1:pixelskip:end,1:pixelskip:end)));
    134                 md.radaroverlay.x=(x0:(x1-x0)/(size(md.radaroverlay.pwr,2)-1):x1);
    135                 md.radaroverlay.y=(y0:(y1-y0)/(size(md.radaroverlay.pwr,1)-1):y1);
    136 
    137                 %Erase image or keep it?
    138                 if ~getfieldvalue(options,'keep_image',0),
    139                         system('rm -rf ./temp.tif');
    140                 end
     60        end
     61        if ~found,
     62                error('could not find radar image');
     63        end
    14164
    14265
    143         else
    144                 if ~exist(options,'geotiff_name'),
    145                         error('Need a geotiff for areas outside of Greenland and Antarctica');
    146                 end
    147                
    148                 %Name of image
    149                 inputname='./temp.tif';
    150                 eval(['!gdal_translate -quiet -projwin ' num2str(x0) ' ' num2str(y1) ' ' num2str(x1) ' ' num2str(y0) ' ' geotiff_name ' ' inputname ]);
     66        %Crop radar image from xylim
     67        filename='./temp.tif';
     68        eval(['!gdal_translate -quiet -projwin ' num2str(x0) ' ' num2str(y1) ' ' num2str(x1) ' ' num2str(y0) ' ' geotiff_name ' ' filename ]);
    15169
    152                 %Read in temp.tif:
    153                 im=imread('temp.tif','TIFF');
    154                 %adjust contrast and brightness
    155                 %im=imadjust(im,[a b],[c d]);
    156                 pixelskip=max(1,ceil(posting/((x1-x0)/(size(im,2)))));
    157                 md.radaroverlay.pwr=double(flipud(im(1:pixelskip:end,1:pixelskip:end)));
    158                 md.radaroverlay.x=(x0:(x1-x0)/(size(md.radaroverlay.pwr,2)-1):x1);
    159                 md.radaroverlay.y=(y0:(y1-y0)/(size(md.radaroverlay.pwr,1)-1):y1);
     70        %Read in temp.tif:
     71        im=imread('temp.tif','TIFF');
     72        %adjust contrast and brightness
     73        %im=imadjust(im,[a b],[c d]);
     74        pixelskip=max(1,ceil(posting/((x1-x0)/(size(im,2)))));
     75        md.radaroverlay.pwr=double(flipud(im(1:pixelskip:end,1:pixelskip:end)));
     76        md.radaroverlay.x=(x0:(x1-x0)/(size(md.radaroverlay.pwr,2)-1):x1);
     77        md.radaroverlay.y=(y0:(y1-y0)/(size(md.radaroverlay.pwr,1)-1):y1);
    16078
    161                 %Erase image or keep it?
    162                 if ~getfieldvalue(options,'keep_image',0),
    163                         system('rm -rf ./temp.tif');
    164                 end
    165         end
     79        %Erase image or keep it?
     80        if ~getfieldvalue(options,'keep_image',0),
     81                system('rm -rf ./temp.tif');
     82        end
    16683        %}}}
    16784else %user provided image {{{
Note: See TracChangeset for help on using the changeset viewer.