Changeset 13127
- Timestamp:
- 08/22/12 14:50:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/m/plot/radarpower.m
r13009 r13127 10 10 % md=radarpower(md) 11 11 12 %If gdal does not work, uncomment the following line13 %setenv('LD_LIBRARY_PATH','/proj/ice/larour/issm/trunk/externalpackages/gdal/install/lib/');14 12 %Parse inputs 15 13 if nargin==1, … … 22 20 end 23 21 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 default22 highres = getfieldvalue(options,'highres',0); 23 xlim = getfieldvalue(options,'xlim',[min(md.mesh.x) max(md.mesh.x)]); 24 ylim = getfieldvalue(options,'ylim',[min(md.mesh.y) max(md.mesh.y)]); 25 posting = getfieldvalue(options,'posting',0); % 0 -> image posting default 28 26 29 27 %find gdal coordinates … … 34 32 if ~exist(options,'overlay_image'), 35 33 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']; 38 76 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');44 77 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 ]); 50 81 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); 56 88 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 64 92 65 93 elseif strcmpi(md.mesh.hemisphere,'s'), … … 94 122 end 95 123 else 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. 97 125 if (~exist(options,'overlay_xlim') | ~exist(options,'overlay_xlim')| ~exist(options,'overlay_xposting')| ~exist(options,'overlay_yposting')), 98 126 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.