Changeset 13127


Ignore:
Timestamp:
08/22/12 14:50:00 (13 years ago)
Author:
Mathieu Morlighem
Message:

NEW: use MOG geotif (500 or 100 m res) instead of jpeg

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r13009 r13127  
    1010%      md=radarpower(md)
    1111
    12 %If gdal does not work, uncomment the following line
    13 %setenv('LD_LIBRARY_PATH','/proj/ice/larour/issm/trunk/externalpackages/gdal/install/lib/');
    1412%Parse inputs
    1513if nargin==1,
    … …  
    2220end
    2321
    24 highres=getfieldvalue(options,'highres',0);
    25 xlim=getfieldvalue(options,'xlim',[min(md.mesh.x) max(md.mesh.x)]);
    26 ylim=getfieldvalue(options,'ylim',[min(md.mesh.y) max(md.mesh.y)]);
    27 posting=getfieldvalue(options,'posting',0); % 0 -> image posting default
     22highres = getfieldvalue(options,'highres',0);
     23xlim    = getfieldvalue(options,'xlim',[min(md.mesh.x) max(md.mesh.x)]);
     24ylim    = getfieldvalue(options,'ylim',[min(md.mesh.y) max(md.mesh.y)]);
     25posting = getfieldvalue(options,'posting',0); % 0 -> image posting default
    2826
    2927%find gdal coordinates
    … …  
    3432if ~exist(options,'overlay_image'),
    3533        if strcmpi(md.mesh.hemisphere,'n'),
    36                 if ~exist([jplsvn() '/projects/ModelData/MOG/mog150_greenland_map.jpg']),
    37                         error(['radarpower error message: file ' jplsvn() '/projects/ModelData/MOG/mog150_greenland_map.jpg not found.']);
     34                %if ~exist([jplsvn() '/projects/ModelData/MOG/mog150_greenland_map.jpg']),
     35                %       error(['radarpower error message: file ' jplsvn() '/projects/ModelData/MOG/mog150_greenland_map.jpg not found.']);
     36                %end
     37                %name = 'mog150_greenland_map';
     38                %name = 'mog100_hp1_v10';
     39                %%name = 'mog500_hp1_v10';
     40                %jpgim=[jplsvn() '/projects/ModelData/MOG/' name '.jpg'];
     41                %geom=load([jplsvn() '/projects/ModelData/MOG/' name '.jpgw'],'ascii');
     42                %%jpgim='/u/astrid-r1b/morlighe/issmjpl/projects/MorlighemGRL2012/Data/Mosaic_amp_asar2010.jpg';
     43                %%geom=load('/u/astrid-r1b/morlighe/issmjpl/projects/MorlighemGRL2012/Data/Mosaic_amp_asar2010.jpgw');
     44                %jpgim='/u/astrid-r1b/morlighe/issmjpl/projects/MorlighemGRL2012/Data/Russel_asar2010.png';
     45                %geom=load('/u/astrid-r1b/morlighe/issmjpl/projects/MorlighemGRL2012/Data/Russel_asar2010.pngw');
     46
     47                %%geom:   xposting nbcols nbrows yposting xmin ymax
     48                %xmin=max(geom(5),x0);
     49                %xmax=min(geom(5)+geom(1)*geom(2),x1);
     50                %ymin=max(geom(6)-geom(3)*geom(4),y0);
     51                %ymax=min(geom(6),y1);
     52
     53                %firstcol=max(1,floor((xmin-geom(5))/geom(1))); %x min
     54                %firstrow=max(1,floor((geom(6)-ymax)/geom(4))); %y max
     55                %numcols=floor((xmax-xmin)/geom(1)); % x posting
     56                %numrows=floor((ymax-ymin)/geom(4)); % y posting
     57                %pixelskip=max(1,ceil(posting/geom(1)));
     58
     59                %%Read and crop file
     60                %disp('Warning: expecting coordinates in polar stereographic (Std Latitude: 70ºN Meridian: 45º)');
     61                %im=imread(jpgim);
     62                %im=im(firstrow:firstrow+numrows-1,firstcol:firstcol+numcols-1);
     63                %md.radaroverlay.pwr=double(flipud(im(1:pixelskip:end,1:pixelskip:end)));
     64                %md.radaroverlay.x=(xmin:(xmax-xmin)/(size(md.radaroverlay.pwr,2)-1):xmax);
     65                %md.radaroverlay.y=(ymin:(ymax-ymin)/(size(md.radaroverlay.pwr,1)-1):ymax);
     66                if highres,
     67                        if ~exist([jplsvn() '/projects/ModelData/MOG/mog100g_r2_hp1.tif']),
     68                                error(['radarpower error message: file ' jplsvn() '/projects/ModelData/MOG/mog100g_r2_hp1.tif not found.']);
     69                        end
     70                        geotiff_name=[jplsvn() '/projects/ModelData/MOG/mog100_r2_hp1.tif'];
     71                else
     72                        if ~exist([jplsvn() '/projects/ModelData/MOG/mog500g_r2_hp1.tif']),
     73                                error(['radarpower error message: file ' jplsvn() '/projects/ModelData/MOG/mog500g_r2_hp1.tif not found.']);
     74                        end
     75                        geotiff_name=[jplsvn() '/projects/ModelData/MOG/mog500_r2_hp1.tif'];
    3876                end
    39                 name = 'mog150_greenland_map';
    40                 %name = 'mog100_hp1_v10';
    41                 %name = 'mog500_hp1_v10';
    42                 jpgim=[jplsvn() '/projects/ModelData/MOG/' name '.jpg'];
    43                 geom=load([jplsvn() '/projects/ModelData/MOG/' name '.jpgw'],'ascii');
    4477
    45                 %geom:   xposting nbcols nbrows yposting xmin ymax
    46                 xmin=max(geom(5),x0);
    47                 xmax=min(geom(5)+geom(1)*geom(2),x1);
    48                 ymin=max(geom(6)-geom(3)*geom(4),y0);
    49                 ymax=min(geom(6),y1);
     78                %Name of image
     79                inputname='./temp.tif';
     80                eval(['!gdal_translate -quiet -projwin ' num2str(x0) ' ' num2str(y1) ' ' num2str(x1) ' ' num2str(y0) ' ' geotiff_name ' ' inputname ]);
    5081
    51                 firstcol=max(1,floor((xmin-geom(5))/geom(1))); %x min
    52                 firstrow=max(1,floor((geom(6)-ymax)/geom(4))); %y max
    53                 numcols=floor((xmax-xmin)/geom(1)); % x posting
    54                 numrows=floor((ymax-ymin)/geom(4)); % y posting
    55                 pixelskip=max(1,ceil(posting/geom(1)));
     82                %Read in temp.tif:
     83                im=imread('temp.tif','TIFF');
     84                pixelskip=max(1,ceil(posting/((x1-x0)/(size(im,2)))));
     85                md.radaroverlay.pwr=double(flipud(im(1:pixelskip:end,1:pixelskip:end)));
     86                md.radaroverlay.x=(x0:(x1-x0)/(size(md.radaroverlay.pwr,2)-1):x1);
     87                md.radaroverlay.y=(y0:(y1-y0)/(size(md.radaroverlay.pwr,1)-1):y1);
    5688
    57                 %Read and crop file
    58                 disp('Warning: expecting coordinates in polar stereographic (Std Latitude: 70ºN Meridian: 45º)');
    59                 im=imread(jpgim);
    60                 im=im(firstrow:firstrow+numrows-1,firstcol:firstcol+numcols-1);
    61                 md.radaroverlay.pwr=double(flipud(im(1:pixelskip:end,1:pixelskip:end)));
    62                 md.radaroverlay.x=(xmin:(xmax-xmin)/(size(md.radaroverlay.pwr,2)-1):xmax);
    63                 md.radaroverlay.y=(ymin:(ymax-ymin)/(size(md.radaroverlay.pwr,1)-1):ymax);
     89                %Erase image
     90                system('rm -rf ./temp.tif');
     91
    6492
    6593        elseif strcmpi(md.mesh.hemisphere,'s'),
    … …  
    94122        end
    95123else
    96         %ok, user provided an image. check we also have overlay_xlim and overlay_ylim  options, to know what range of coordinates the image covers.
     124        %user provided an image. check we also have overlay_xlim and overlay_ylim  options, to know what range of coordinates the image covers.
    97125        if (~exist(options,'overlay_xlim') | ~exist(options,'overlay_xlim')| ~exist(options,'overlay_xposting')| ~exist(options,'overlay_yposting')),
    98126                error('radarpower error message: please provide overlay_xlim, overlay_ylim, overlay_xposting and overlay_yposting options together with overlay_image option');
Note: See TracChangeset for help on using the changeset viewer.