source:
issm/oecreview/Archive/14064-14311/ISSM-14194-14195.diff
Last change on this file was 14312, checked in by , 12 years ago | |
---|---|
File size: 6.4 KB |
-
../trunk-jpl/src/m/plot/plot_gridded.m
23 23 %Interpolating data on grid 24 24 [x_m y_m data_grid]=InterpFromMeshToGrid(elements,x,y,data,xlim(1),ylim(2),post,post,round(diff(ylim)/post),round(diff(xlim)/post),NaN); 25 25 if size(data_grid,1)<3 | size(data_grid,2)<3, 26 error('data_grid size too small in plot_gridded, check posting and uni ');26 error('data_grid size too small in plot_gridded, check posting and units'); 27 27 end 28 28 29 29 %Get and change colormap -
../trunk-jpl/src/m/plot/kmlgridded.m
1 function kmlgridded(md,data,varargin) 2 3 %process options 4 options=pairoptions(varargin{:}); 5 6 %process options 7 options=changefieldvalue(options,'coord','latlon'); 8 9 %process mesh and data 10 [x y z elements is2d isplanet]=processmesh(md,[],options); 11 [data datatype]=processdata(md,data,options); 12 13 %check is2d 14 if ~is2d, 15 error('buildgridded error message: gridded not supported for 3d meshes, project on a layer'); 16 end 17 18 %Get xlim and ylim (used to extract radar image) 19 xlim=[min(x) max(x)]; 20 ylim=[min(y) max(y)]; 21 post=getfieldvalue(options,'posting',diff(xlim)/1000); 22 if(diff(xlim)/post>10000), 23 error(['posting too large']); 24 end 25 26 %Interpolating data on grid 27 [x_m y_m data_grid]=InterpFromMeshToGrid(elements,x,y,data,xlim(1),ylim(2),post,post,round(diff(ylim)/post),round(diff(xlim)/post),NaN); 28 if size(data_grid,1)<3 | size(data_grid,2)<3, 29 error('data_grid size too small, check posting and units'); 30 end 31 pos=find(isinf(data_grid)); 32 if ~isempty(pos), 33 disp('Warning: removing Infs from vector (probably log(0)?)'); 34 data_grid(pos)=NaN; 35 end 36 37 %Process data_grid: add white in NaN and correct caxis accordingly 38 data_nan=find(isnan(data_grid)); 39 data_min=min(data_grid(:)); 40 data_max=max(data_grid(:)); 41 if exist(options,'caxis'), 42 caxis_opt=getfieldvalue(options,'caxis'); 43 data_grid(find(data_grid<caxis_opt(1)))=caxis_opt(1); 44 data_grid(find(data_grid>caxis_opt(2)))=caxis_opt(2); 45 data_min=caxis_opt(1); 46 data_max=caxis_opt(2); 47 end 48 49 %Get colormap 50 colorm = getcolormap(options); 51 len = size(colorm,1); 52 ind = ceil((len-1)*(data_grid-data_min)/(data_max - data_min + eps) +1); 53 ind(find(ind>len))=len; 54 ind(find(ind<1) )=1; 55 ind(find(isnan(ind)))=1; 56 image_rgb=zeros(size(data_grid,1),size(data_grid,2),3); 57 r=colorm(:,1); image_rgb(:,:,1)=r(ind); clear r; 58 g=colorm(:,2); image_rgb(:,:,2)=g(ind); clear g; 59 b=colorm(:,3); image_rgb(:,:,3)=b(ind); clear b; 60 61 %Deal with alpha 62 alpha=getfieldvalue(options,'alpha',.8); 63 alphaMatrix = alpha*ones(size(data_grid)); 64 alphaMatrix(data_nan) = 0; 65 66 %write kml 67 kmlfilename=getfieldvalue(options,'kmlfilename','tempfile.kml'); 68 kmlroot=getfieldvalue(options,'kmlroot','./'); 69 kmlimagename=getfieldvalue(options,'kmlimagename','tempimage'); 70 kmlresolution=getfieldvalue(options,'kmlresolution',1); 71 kmlfolder=getfieldvalue(options,'kmlfolder','Ground Overlay'); 72 kmlfolderdescription=getfieldvalue(options,'kmlfolderdescription',''); 73 kmlgroundoverlayname=getfieldvalue(options,'kmlgroundoverlayname','ground overlay'); 74 kmlgroundoverlaydescription=getfieldvalue(options,'kmlgroundoverlaydescription','description'); 75 76 %write png 77 imwrite(image_rgb,[kmlimagename '.png'],'png','alpha',alphaMatrix); 78 clear image_rgb alphaMatrix 79 80 %prepare colorbar 81 iscolorbar=0; 82 if strcmpi(getfieldvalue(options,'colorbar','on'),'on'), 83 X = linspace(0,1,len)'; 84 Xlab = round(linspace(data_min,data_max,len+1)); 85 html = ['<TABLE border=' num2str(1) ' bgcolor=#FFFFFF>',10]; 86 87 for k=len:-1:1 88 f = (Xlab(k)-data_min)/(data_max-data_min); 89 if f<0, f=0; end 90 if f>1, f=1; end 91 polyColor(1,1) = interp1(X,colorm(:,1),f); 92 polyColor(1,2) = interp1(X,colorm(:,2),f); 93 polyColor(1,3) = interp1(X,colorm(:,3),f); 94 polyColorStr(1:2) = dec2hex(round(polyColor(1)*255),2); 95 polyColorStr(3:4) = dec2hex(round(polyColor(2)*255),2); 96 polyColorStr(5:6) = dec2hex(round(polyColor(3)*255),2); 97 html = [html,'<TR><TD width="15px" bgcolor=#',polyColorStr, '> </TD>','<TD bgcolor=#FFFFFF>']; 98 if k==1 99 html=[html,'<= ',num2str(Xlab(k),'%g')]; 100 elseif k==len 101 html=[html,'>= ',num2str(Xlab(k),'%g')]; 102 else 103 html=[html,num2str(Xlab(k),'%g'),' to ',num2str(Xlab(k+1),'%g'),'</TD>']; 104 end 105 html = [html,'</TR>',10]; 106 end 107 html = [html,'</TABLE>']; 108 iscolorbar = 1; 109 end 110 111 %now write kml file 112 fid=fopen([kmlroot '/' kmlfilename],'w'); 113 fprintf(fid,'%s\n','<?xml version="1.0" encoding="UTF-8"?>'); 114 fprintf(fid,'%s\n','<kml xmlns="http://earth.google.com/kml/2.1">'); 115 fprintf(fid,'%s\n','<Document>'); 116 fprintf(fid,'%s%s%s\n','<name>',kmlfilename,'</name>'); 117 if iscolorbar, 118 fprintf(fid,'<Placemark id="colorbar">\n'); 119 fprintf(fid,'%s%s%s\n','<name>','click the icon to see the colorbar','</name>'); 120 fprintf(fid,'%s%s%s\n','<description>','Ground overlay colorbar','</description>'); 121 fprintf(fid,'<visibility>1</visibility>\n'); 122 fprintf(fid,['<description>',10,'<![CDATA[' html ']]>',10,'</description>',10,'\n']); 123 fprintf(fid,['<Style><IconStyle><scale>1</scale><Icon><href>http://maps.google.com/mapfiles/kml/shapes/donut.png</href></Icon></IconStyle><ListStyle></ListStyle></Style><Point id="poly_colorbar">\n']); 124 fprintf(fid,'<altitudeMode>clampToGround</altitudeMode>\n'); 125 fprintf(fid,'<extrude>1</extrude>\n'); 126 fprintf(fid,'<tessellate>1</tessellate>\n'); 127 fprintf(fid,'%s%g,%g%s\n','<coordinates>',max(x),mean(y),'</coordinates>'); 128 fprintf(fid,'</Point>\n'); 129 fprintf(fid,'</Placemark>\n'); 130 end 131 fprintf(fid,'%s\n','<GroundOverlay id="groundoverlay">'); 132 fprintf(fid,'%s%s%s\n','<name>',kmlgroundoverlayname,'</name>'); 133 fprintf(fid,'%s\n','<description>',kmlgroundoverlaydescription,'</description>'); 134 fprintf(fid,'%s%s.%s%s\n','<Icon>',kmlimagename,'png','</Icon>'); 135 fprintf(fid,'%s\n','<LatLonBox>'); 136 fprintf(fid,'%s%f%s\n','<north>',max(y_m),'</north>'); 137 fprintf(fid,'%s%f%s\n','<south>',min(y_m),'</south>'); 138 fprintf(fid,'%s%f%s\n','<east>',max(x_m),'</east>'); 139 fprintf(fid,'%s%f%s\n','<west>',min(x_m),'</west>'); 140 fprintf(fid,'%s\n','<rotation>0</rotation>'); 141 fprintf(fid,'%s\n','</LatLonBox>'); 142 fprintf(fid,'%s\n','</GroundOverlay>'); 143 fprintf(fid,'%s\n','</Document>'); 144 fprintf(fid,'%s\n','</kml>'); 145 fclose(fid);
Note:
See TracBrowser
for help on using the repository browser.