function [Z,t1,t2] = xyztogr(X, lab)
% [Z,t1,t2]  = xyztogr(X,lab)
% converts a matrix of x (E-W) coordinates, y (N-S) 
% coordinates, and z values to a grid Z
% x and y need not be integers, or range from 1 up
% Z is a matrix of the gridded z values, t1 and t2
% are vectors which can be used as scales for the
% x and y axes.	Also plots contour and mesh diagrams
% To label contours set lab to 1 (auto) or 2 (manual,
% using mouse, otherwise 0 (default)
% uses contvec to determine contour levels
% Written by Gerry Middleton, Feb 1997
if nargin < 2, lab = 0; end
eps = 0.000001;
x1min = min(X(:,1));  % determine maxima and minima for x (x1)
x1max = max(X(:,1));
x2min = min(X(:,2));  % and for y (x2)
x2max = max(X(:,2));
zmin = min(X(:,3));
zmax = max(X(:,3));
xx = sort(X(:,1));	  % determine increments
inc1 = max(diff(xx));
xx = sort(X(:,2));
inc2 = max(diff(xx));
t1 = x1min:inc1:x1max+eps; % make x1 and x2 vectors
t2 = x2min:inc2:x2max+eps; % eps used to avoid rounding errors
ng1 = length(t1);	   % find their size; ng1 is no of cols
ng2 = length(t2);	   % ng2 is no of rows
ci = round((X(:,1) - x1min)/inc1 + 1);	% convert to integer indices
ri = round((X(:,2) - x2min)/inc2 + 1);
Z = full(sparse(ri,ci,X(:,3))); % make array
v = contvec(zmin,zmax);
vc = length(v);
fprintf('Max value in data: %4.1f  ',zmax);
fprintf('Min value  %4.1f  \n',zmin);
fprintf('There are %2.0f equally spaced contour lines \n',vc);
fprintf('Press <Enter> to continue, at each pause \n');
fprintf('\n');
pause; 
c = contour(t1,t2,Z,v);
if lab == 1 
   clabel(c);
end
else if lab == 2
   clabel(c,'manual');
end
xlabel('x'), ylabel('y');
axis('equal');
pause;
figure
mesh(t1,t2,Z);
xlabel('X1'),ylabel('X2');
i = [ng2:-1:1];
Z = Z(i,:);  % reverse order of rows

