Index: /issm/trunk/src/m/kml/kml_mesh_write.m
===================================================================
--- /issm/trunk/src/m/kml/kml_mesh_write.m	(revision 6278)
+++ /issm/trunk/src/m/kml/kml_mesh_write.m	(revision 6278)
@@ -0,0 +1,298 @@
+%
+%  write a kml file of the mesh from the model.
+%
+%  []=kml_mesh_write(filek,md,params)
+%
+%  where the required input is:
+%    filek         (char, name of .kml file)
+%    md            (model, model class object)
+%
+%  the optional input is:
+%    params        (string/numeric, parameter names and values)
+%
+%  and the optional input is:
+%    data          (numeric, results data)
+%    alt           (numeric, altitude for polygons, default 10000)
+%    lwidth        (numeric, line width in pixels, default 1)
+%    popac         (numeric, polygon opacity, default 0.50)
+%    cmin          (numeric, minimum of color scale)
+%    cmax          (numeric, maximum of color scale)
+%    cmap          (char or numeric, colormap definition)
+%
+function []=kml_mesh_write(varargin)
+
+if ~nargin
+    help kml_mesh_write
+    return
+end
+
+%%  process input data
+
+iarg=1;
+if (nargin >= 1)
+    filek=varargin{1};
+end
+if ~exist('filek' ,'var') || isempty(filek)
+    filek=input('kml file to write?  ','s');
+end
+[pathstr,name,ext,versn] = fileparts(filek);
+if isempty(ext)
+    ext='.kml';
+end
+filek2=fullfile(pathstr,[name ext versn]);
+
+display(sprintf('Opening kml file ''%s''.',filek2));
+fid=fopen(sprintf('%s',filek2),'w');
+if (fid < 0)
+    error('File ''%s'' could not be opened.',filek2);
+end
+
+iarg=iarg+1;
+if (nargin >= 2)
+    md=varargin{2};
+end
+if ~exist('md','var') || isempty(md) || ~isa(md,'model')
+    error(['Model ''' inputname(iarg) ''' is unrecognized class ''' class(md) '''.']);
+end
+
+%  parameters
+
+iarg=iarg+1;
+while (iarg <= nargin-1)
+    if ischar(varargin{iarg})
+        eval([varargin{iarg} '=varargin{iarg+1};']);
+        if (numel(varargin{iarg+1}) <= 20)
+            disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']);
+        else
+            disp([varargin{iarg} '=' string_size(varargin{iarg+1}) ' ' class(varargin{iarg+1}) ';']);
+        end
+        if strcmpi(varargin{iarg},'data')
+            cdata=inputname(iarg+1);
+        end
+    else
+        error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']);
+    end
+    iarg=iarg+2;
+end
+
+if exist('data','var') && ~isempty(data)
+    if     (numel(data)==md.numberofelements)
+        edata=data;
+    elseif (numel(data)==md.numberofgrids)
+        display('Averaging nodal data to element data.');
+        edata=zeros(1,md.numberofelements);
+        for i=1:size(md.elements,1)
+            for j=1:size(md.elements,2)
+                edata(i)=edata(i)+data(md.elements(i,j));
+            end
+            edata(i)=edata(i)/size(md.elements,2);
+        end
+    else
+        error(['Data has incorrect number of ' num2str(numel(data)) ' elements.']);
+    end
+end
+
+%%  write kml file
+
+%  write header data
+
+fprintf(fid,'<?xml version="1.0" encoding="UTF-8"?>\n');
+fprintf(fid,'<kml xmlns="http://www.opengis.net/kml/2.2">\n');
+fprintf(fid,'  <Document>\n');
+fprintf(fid,'    <name>ISSM Mesh: %s</name>\n',md.name);
+fprintf(fid,'    <open>1</open>\n');
+fprintf(fid,'    <description>');
+ifirst=true;
+for i=1:numel(md.notes)
+    if ~isempty(md.notes{i})
+        if ~ifirst
+            fprintf(fid,'\n');
+        end
+        ifirst=false;
+        fprintf(fid,'%s',md.notes{i});
+    end
+end
+fprintf(fid,'</description>\n');
+
+%  write style templates for defaults and for each color of the matlab
+%  colormap (note that matlab colormap format is rgb, where each varies
+%  from 0 to 1, whereas the kml color format is aabbggrr, where each
+%  varies from 00 to ff.)
+
+if ~exist('lwidth','var')
+    lwidth=1;
+end
+if ~exist('popac','var')
+    popac=0.50;
+end
+
+fprintf(fid,'    <Style id="BlackLineRandomPoly">\n');
+fprintf(fid,'      <LineStyle>\n');
+fprintf(fid,'        <color>ff000000</color>\n');
+fprintf(fid,'        <colorMode>normal</colorMode>\n');
+fprintf(fid,'        <width>%g</width>\n',lwidth);
+fprintf(fid,'      </LineStyle>\n');
+fprintf(fid,'      <PolyStyle>\n');
+fprintf(fid,'        <color>%02xffffff</color>\n',round(popac*255));
+fprintf(fid,'        <colorMode>random</colorMode>\n');
+fprintf(fid,'      </PolyStyle>\n');
+fprintf(fid,'    </Style>\n');
+fprintf(fid,'    <Style id="BlackLineEmptyPoly">\n');
+fprintf(fid,'      <LineStyle>\n');
+fprintf(fid,'        <color>ff000000</color>\n');
+fprintf(fid,'        <colorMode>normal</colorMode>\n');
+fprintf(fid,'        <width>%g</width>\n',lwidth);
+fprintf(fid,'      </LineStyle>\n');
+fprintf(fid,'      <PolyStyle>\n');
+fprintf(fid,'        <color>00ffffff</color>\n');
+fprintf(fid,'        <colorMode>random</colorMode>\n');
+fprintf(fid,'      </PolyStyle>\n');
+fprintf(fid,'    </Style>\n');
+fprintf(fid,'    <Style id="RedLineRedPoly">\n');
+fprintf(fid,'      <LineStyle>\n');
+fprintf(fid,'        <color>ff0000ff</color>\n');
+fprintf(fid,'        <colorMode>normal</colorMode>\n');
+fprintf(fid,'        <width>%g</width>\n',lwidth);
+fprintf(fid,'      </LineStyle>\n');
+fprintf(fid,'      <PolyStyle>\n');
+fprintf(fid,'        <color>%02x0000ff</color>\n',round(popac*255));
+fprintf(fid,'        <colorMode>random</colorMode>\n');
+fprintf(fid,'      </PolyStyle>\n');
+fprintf(fid,'    </Style>\n');
+if exist('edata','var')
+
+%  colormap command operates on a figure, so create an invisible one
+%  (could also directly call colormaps, e.g. jet(64), but risky)
+
+    hfig=figure('Visible','off');
+    if exist('cmap','var')
+        colormap(cmap)
+    end
+    cmap=colormap;
+    close(hfig)
+    
+    disp(['Writing ' num2str(size(cmap,1)) ' Matlab colors as KML style templates.']);
+    for i=1:size(cmap,1)
+        fprintf(fid,'    <Style id="MatlabColor%d">\n',i);
+        fprintf(fid,'      <LineStyle>\n');
+        fprintf(fid,'        <color>ff000000</color>\n');
+        fprintf(fid,'        <colorMode>normal</colorMode>\n');
+        fprintf(fid,'        <width>%g</width>\n',lwidth);
+        fprintf(fid,'      </LineStyle>\n');
+        fprintf(fid,'      <PolyStyle>\n');
+        fprintf(fid,'        <color>%02x%02x%02x%02x</color>\n',round(popac*255),...
+            round(cmap(i,3)*255),round(cmap(i,2)*255),round(cmap(i,1)*255));
+        fprintf(fid,'        <colorMode>normal</colorMode>\n');
+        fprintf(fid,'      </PolyStyle>\n');
+        fprintf(fid,'    </Style>\n');
+    end
+end
+
+%  write folder for mesh
+
+fprintf(fid,'    <Folder>\n');
+if exist('cdata','var') && ~isempty(cdata)
+    fprintf(fid,'      <name>Data: %s</name>\n',cdata);
+else
+    fprintf(fid,'      <name>Mesh</name>\n');
+end
+fprintf(fid,'      <visibility>1</visibility>\n');
+fprintf(fid,'      <description>Elements=%d, Grids=%d</description>\n',...
+    md.numberofelements,md.numberofgrids);
+
+%  write each element as a polygon
+
+if exist('edata','var')
+    if ~exist('cmin','var')
+        cmin=min(min(edata));
+    end
+    if ~exist('cmax','var')
+        cmax=max(max(edata));
+    end
+end
+
+if ~exist('alt','var')
+    alt=10000;
+end
+
+disp(['Writing ' num2str(size(md.elements,1)) ' tria elements as KML polygons.']);
+for i=1:size(md.elements,1)
+    fprintf(fid,'      <Placemark>\n');
+    fprintf(fid,'        <name>Element %d</name>\n',i);
+    fprintf(fid,'        <visibility>1</visibility>\n');
+    if exist('edata','var')
+        fprintf(fid,'        <description>Element data: %g</description>\n',edata(i));
+        imap = fix((edata(i)-cmin)/(cmax-cmin)*size(cmap,1))+1;
+        if     (imap >= 1) && (imap <= size(cmap,1))
+            fprintf(fid,'        <styleUrl>#MatlabColor%d</styleUrl>\n',imap);
+        elseif (edata(i) == cmax)
+            fprintf(fid,'        <styleUrl>#MatlabColor%d</styleUrl>\n',size(cmap,1));
+        else
+            fprintf(fid,'        <styleUrl>#BlackLineEmptyPoly</styleUrl>\n');
+        end
+    else
+        fprintf(fid,'        <styleUrl>#BlackLineRandomPoly</styleUrl>\n');
+    end
+    fprintf(fid,'        <Polygon>\n');
+    fprintf(fid,'          <extrude>1</extrude>\n');
+    fprintf(fid,'          <altitudeMode>relativeToGround</altitudeMode>\n');
+    fprintf(fid,'          <outerBoundaryIs>\n');
+    fprintf(fid,'            <LinearRing>\n');
+    fprintf(fid,'              <coordinates>\n');
+    for j=1:size(md.elements,2)
+        [lat(j),long(j)]=mapxy(md.x(md.elements(i,j)),md.y(md.elements(i,j)),'s');
+        fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(j),lat(j),alt);
+    end
+    fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(1),lat(1),alt);
+
+    fprintf(fid,'              </coordinates>\n');
+    fprintf(fid,'            </LinearRing>\n');
+    fprintf(fid,'          </outerBoundaryIs>\n');
+    fprintf(fid,'        </Polygon>\n');
+    fprintf(fid,'      </Placemark>\n');
+end
+fprintf(fid,'    </Folder>\n');
+
+%  write folder for partition segments
+
+if md.npart
+    [xseg,yseg]=flagedges(md.elements,md.x,md.y,md.part);
+    fprintf(fid,'    <Folder>\n');
+    fprintf(fid,'      <name>Segments</name>\n');
+    fprintf(fid,'      <visibility>1</visibility>\n');
+    fprintf(fid,'      <description>Partitions=%d, Segments=%d</description>\n',...
+        md.npart,size(xseg,1));
+
+%  write each element as a polygon
+
+    disp(['Writing ' num2str(size(xseg,1)) ' partition segments as KML linestrings.']);
+    for i=1:size(xseg,1)
+        fprintf(fid,'      <Placemark>\n');
+        fprintf(fid,'        <name>Segment %d</name>\n',i);
+        fprintf(fid,'        <visibility>1</visibility>\n');
+        fprintf(fid,'        <styleUrl>#RedLineRedPoly</styleUrl>\n');
+        fprintf(fid,'        <LineString>\n');
+        fprintf(fid,'          <extrude>1</extrude>\n');
+        fprintf(fid,'          <tessellate>1</tessellate>\n');
+        fprintf(fid,'          <altitudeMode>relativeToGround</altitudeMode>\n');
+        fprintf(fid,'          <coordinates>\n');
+        for j=1:2
+            [lat(j),long(j)]=mapxy(xseg(i,j),yseg(i,j),'s');
+            fprintf(fid,'            %0.16g,%0.16g,%0.16g\n',long(j),lat(j),alt);
+        end
+
+        fprintf(fid,'          </coordinates>\n');
+        fprintf(fid,'        </LineString>\n');
+        fprintf(fid,'      </Placemark>\n');
+    end
+    fprintf(fid,'    </Folder>\n');
+end
+
+%  write trailer data
+
+fprintf(fid,'  </Document>\n');
+fprintf(fid,'</kml>\n');
+
+fclose(fid);
+display('End of file successfully written.');
+
