% First stab at creating an indexed-face-set from a netcdf file
% Matlab must have the netcdf tools installed,
% see http://crusty.er.usgs.gov/mexcdf.html if needed
% Created 11/25/98 by C. Moore
% (use the next 3 lines to initialize netcdf interaction)
%cd ~/bin
%ncstartup
%cd ~/data/temp
clear
for n=1:11
n
%
% open netcdf file, and create the vrml header
%
indir='../iso/lat_lon/';
outdir='data/';
if n<10
 filei=[indir,'980',int2str(n),'_iso20_xy_lf.cdf'];
 fileo=[outdir,'980',int2str(n),'_iso20gs.wrl'];
else
 filei=[indir,'98',int2str(n),'_iso20_xy_lf.cdf'];
 fileo=[outdir,'98',int2str(n),'_iso20gs.wrl'];
end
f=netcdf(filei,'nowrite');
ncclose(filei);
data=f{6}(:);
lat=f{4}(:);
lon=f{5}(:);
%
% Print header to output vrml (.wrl) file
%
fid=fopen(fileo,'w');
fprintf(fid,'#VRML V2.0 utf8\n');
fprintf(fid,'# created by mkiso.m - Chris Moore\n');
fprintf(fid,'# 11/30/98\n');
fprintf(fid,'Transform {\n');
fprintf(fid,'   children [\n');
fprintf(fid,'	Shape {\n');
% uncomment the next line for transparency
fprintf(fid,'   appearance Appearance {material Material{transparency 0.6}}\n');
fprintf(fid,'	geometry IndexedFaceSet {\n');
fprintf(fid,'	   solid FALSE\n');
fprintf(fid,'	   coord Coordinate {\n');
fprintf(fid,'		point [\n');
%
% do not use bad data: in this case, data greater than 1000 are bad
%
BADDATA=1.0e4;
%
% lat/lon to vrml coordinates (adjust scale to taste)
scale=4;
x=0:scale:(length(lon)-1)*scale;
y=0:-scale:-(length(lat)-1)*scale;
%
% Save the points to vrml file
%  note: we 'invert' the data to give depth in the -y direction
%
INVERT = -1.0;
pointnum=ones(length(lon),length(lat));
count=0;
for l=1:length(lat),
 for k=1:length(lon),
  if (data(l,k)<BADDATA)
   fprintf(fid,' %1.3f %1.3f %1.3f,',x(k),INVERT*data(l,k),y(l));
   pointnum(k,l)=count;
   count=count+1;
  end
 end
 fprintf(fid,'\n');
end
%
% Print text for Index section
%
fprintf(fid,'		]\n');
fprintf(fid,'	   }\n');
fprintf(fid,'	   coordIndex [\n');
%
% Generate and save the indices
%
numfacesx=length(lon)-1;
numfacesy=length(lat)-1;
for i=1:numfacesy
 for j=1:numfacesx
  if (data(i,j)<BADDATA)&(data(i,j+1)<BADDATA)&(data(i+1,j+1)<BADDATA)&(data(i+1,j)<BADDATA)
   fprintf(fid,' %d %d %d %d -1,',pointnum(j,i),pointnum(j+1,i),pointnum(j+1,i+1),pointnum(j,i+1));
  end
 end
 fprintf(fid,'\n');
end
%
% Print text for color table section
%
fprintf(fid,'	   ]\n');
fprintf(fid,'	   colorPerVertex TRUE\n');
fprintf(fid,'	   color Color {\n');
fprintf(fid,'	     color [\n');
%
% Print the color table
%
load colortable
for i=1:length(colortable)
 fprintf(fid,' %1.4f %1.4f %1.4f,\n',colortable(i,1),colortable(i,2),colortable(i,3));
end
%
% Print text for color index section
%
fprintf(fid,'	     ]\n');
fprintf(fid,'	   }\n');
fprintf(fid,'	   colorIndex [\n');
%
% Generate and save the scaled color, and indices
%
mincolor=1;
maxcolor=length(colortable)-1;
%q=find(data<BADDATA);
%maxdata=max(max(INVERT*data(q)));
%mindata=min(min(INVERT*data(q)));
% note: use standard max/min to compare different data sets on same color 
%       scale.
maxdata=0;
mindata=-200;
m=(maxcolor-mincolor)/(maxdata-mindata);
b=1-mindata*m;
color=ones(4);
numfacesx=length(lon)-1;
numfacesy=length(lat)-1;
for i=1:numfacesy
 for j=1:numfacesx
  if (data(i,j)<BADDATA)&(data(i,j+1)<BADDATA)&(data(i+1,j+1)<BADDATA)&(data(i+1,j)<BADDATA)
   color(1)=round(m*INVERT*data(i,j)+b);
   color(2)=round(m*INVERT*data(i,j+1)+b);
   color(3)=round(m*INVERT*data(i+1,j+1)+b);
   color(4)=round(m*INVERT*data(i+1,j)+b);
   fprintf(fid,' %d %d %d %d -1,',color(1),color(2),color(3),color(4));
  end
 end
 fprintf(fid,'\n');
end
%
% Print the closing text and close output file
%
fprintf(fid,'	   ]\n');
fprintf(fid,'	 }\n');
fprintf(fid,'     }\n');
fprintf(fid,'  ]\n');
fprintf(fid,'}\n');
fclose(fid);
end