source: issm/oecreview/Archive/14064-14311/ISSM-14194-14195.diff@ 14312

Last change on this file since 14312 was 14312, checked in by Mathieu Morlighem, 12 years ago

Added 14064-14311

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