contourfill

 CONTOURFILL Filled contour plot
        CONTOURFILL(Z) is a filled contour plot of matrix Z treating the values in Z
        as heights above a plane.
        CONTOURFILL(X,Y,Z), where X and Y are vectors, specifies the X- and Y-
        axes used on the plot. X and Y can also be matrices of the same
        size as Z, in which case they specify a surface in an identical
        manner as SURFACE.
        CONTOURFILL(Z,N) and CONTOURFILL(X,Y,Z,N) draw N contour lines, 
        overriding the default automatic value.
        CONTOURFILL(Z,V) and CONTOURFILL(X,Y,Z,V) draw LENGTH(V) contour lines 
        at the values specified in vector V.

0001 function [CS,H] = contourfill(arg1,arg2,arg3,arg4)
0002 % CONTOURFILL Filled contour plot
0003 %        CONTOURFILL(Z) is a filled contour plot of matrix Z treating the values in Z
0004 %        as heights above a plane.
0005 %        CONTOURFILL(X,Y,Z), where X and Y are vectors, specifies the X- and Y-
0006 %        axes used on the plot. X and Y can also be matrices of the same
0007 %        size as Z, in which case they specify a surface in an identical
0008 %        manner as SURFACE.
0009 %        CONTOURFILL(Z,N) and CONTOURFILL(X,Y,Z,N) draw N contour lines,
0010 %        overriding the default automatic value.
0011 %        CONTOURFILL(Z,V) and CONTOURFILL(X,Y,Z,V) draw LENGTH(V) contour lines
0012 %        at the values specified in vector V.
0013 %
0014 
0015 
0016 %  Author: R. Pawlowicz (IOS)  rich@ios.bc.ca
0017 %          14/12/94
0018 
0019 
0020 if (nargin == 4),
0021    x = arg1;
0022    y = arg2;
0023    z = arg3;
0024    nv = arg4;
0025    if (size(y,1)==1), y=y'; end;
0026    if (size(x,2)==1), x=x'; end;
0027    [mz,nz] = size(z);
0028 elseif (nargin == 3),
0029   x = arg1;
0030   y = arg2;
0031   z = arg3;
0032   nv = [];
0033   if (size(y,1)==1), y=y'; end;
0034   if (size(x,2)==1), x=x'; end;
0035   [mz,nz] = size(z);
0036 elseif (nargin == 2),
0037   [mz,nz] = size(arg1);
0038   x = 1:nz;
0039   y = [1:mz]';
0040   z = arg1;
0041   nv = arg2;
0042 elseif (nargin == 1),
0043   [mz,nz] = size(arg1);
0044   x = 1:nz;
0045   y = [1:mz]';
0046   z = arg1;
0047   nv = [];
0048 end
0049 
0050 i = find(finite(z));
0051 minz = min(min(z(i)));
0052 maxz = max(max(z(i)));
0053 
0054 % Generate default contour levels if they aren't specified
0055 if (max(size(nv)) <= 1)
0056   if isempty(nv)
0057     CS=contourc([minz maxz ; minz maxz]);
0058   else
0059     CS=contourc([minz maxz ; minz maxz],nv);
0060   end
0061 
0062   % Find the levels
0063   ii = 1;
0064   nv = [];
0065   while (ii < size(CS,2)),
0066     nv=[nv CS(1,ii)];
0067     ii = ii + CS(2,ii) + 1;
0068   end
0069 end
0070 
0071 % Handle interior holes correctly
0072 draw_min=0;
0073 if any(nv<=minz),
0074  draw_min=1;
0075 end;
0076 
0077 % Get the unique levels
0078 
0079 nv = sort([minz nv maxz]);
0080 zi = [1, find(diff(nv))+1];
0081 nv = nv(zi);
0082 
0083 % Surround the matrix by a very low region to get closed contours, and
0084 % replace any Nan with low numbers as well.
0085 
0086 
0087 zz=[ NaN+ones(1,nz+2) ; NaN+ones(mz,1) z NaN+ones(mz,1) ; NaN+ones(1,nz+2)];
0088 kk=find(isnan(zz(:)));
0089 zz(kk)=minz-1e4*(maxz-minz)+zeros(size(kk));
0090 
0091 xx = [2*x(:,1)-x(:,2), x, 2*x(:,nz)-x(:,nz-1)];
0092 yy = [2*y(1,:)-y(2,:); y; 2*y(mz,:)-y(mz-1,:)];
0093 if (min(size(yy))==1),
0094  CS=contoursurf(xx,yy,zz,nv);
0095 else
0096  CS=contoursurf(xx([ 1 1:mz mz],:),yy(:,[1 1:nz nz]),zz,nv);
0097 end;
0098 
0099 % Find the indices of the curves in the c matrix, and get the
0100 % area of closed curves in order to draw patches correctly.
0101 ii = 1;
0102 ncurves = 0;
0103 I = [];
0104 Area=[];
0105 while (ii < size(CS,2)),
0106   nl=CS(2,ii);
0107   ncurves = ncurves + 1;
0108   I(ncurves) = ii;
0109   x=CS(1,ii+[1:nl]);  % First patch
0110   y=CS(2,ii+[1:nl]);
0111   Area(ncurves)=sum( diff(x).*(y(1:nl-1)+y(2:nl))/2 );
0112   ii = ii + nl + 1;
0113 end
0114 
0115 plot(CS(1,2),CS(2,2),'-');
0116 
0117 % Plot patches in order of decreasing size. This makes sure that
0118 % all the leves get drawn, not matter if we are going up a hill or
0119 % down into a hole. When going down we shift levels though, you can
0120 % tell whether we are going up or down by checking the sign of the
0121 % area (since curves are oriented so that the high side is always
0122 % the same side). Lowest curve is largest and encloses higher data
0123 % always.
0124 
0125 H=[];
0126 [FA,IA]=sort(-abs(Area));
0127 
0128 for jj=IA,
0129  nl=CS(2,I(jj));
0130  lev=CS(1,I(jj));
0131  if (lev ~=minz | draw_min ),
0132    x=CS(1,I(jj)+[1:nl]);  
0133    y=CS(2,I(jj)+[1:nl]); 
0134    if (sign(Area(jj)) ~=sign(Area(IA(1))) ),
0135      kk=find(nv==lev);
0136      if (kk>1+sum(nv<=minz)*(~draw_min)), 
0137       lev=nv(kk-1);
0138      else 
0139       lev=NaN;         % missing data section
0140      end;
0141    end;
0142 
0143    if (finite(lev)),
0144      H=[H;patch(x,y,lev,'facecolor','flat','edgecolor','none')];
0145    else
0146      H=[H;patch(x,y,lev,'facecolor',get(gcf,'color'),'edgecolor','none')];
0147    end;
0148  end;
0149 end;
0150  
0151 
0152  
0153