function [X,Y,Z,NaNLocations,Z1,Z2]=xyz2grid(x,y,z,method,z1,z2)
%
% xyz2grid creates gridded X, Y, Z, data from xyz data
% x,y,z are vectors of x station locations, y station locations, and 
% observed values respectively
% method is 'nearest' 'cubic' 'linear' 'v4'
% z2, z2 are additional vectors at the same locations
% the user is prompted for limits and sampling intervals
%
% [X,Y,Z,NaNLocations, Z1,Z2]=xyz2grid(x,y,z,method,z1,z2)


%***********************************************************************************
%
%                              Jim Merriam
%                       University of Saskatchewan
%                          jim.merriam@usask.ca
%                              (Dec, 1996)  
%                            (last May 2004)
%***********************************************************************************
%

figure(2),clf
figure(1),clf
plot(x,y,'k.','MarkerSize',7),title('STATIONS')
   % 
   % this block defines the grid dimensions
   % 
disp(' DEFINE THE MAP WINDOW')
disp(' SELECTIONS SHOULD BE WITHIN THE BOUNDARIES IN THE xyz DATA')
disp(' ')
disp(' DEFINE A RECTANGLE BY CLICKING ON BOTTOM LEFT AND TOP RIGHT CORNERS')
[xm,ym]=ginput(2);
xmin=xm(1);
xmax=xm(2);
ymin=ym(1);
ymax=ym(2);
% truncate the data set
i=find(x>=xmin & x <=xmax & y>=ymin & y <= ymax);
disp(['THERE ARE ',num2str(length(x)),' SAMPLES IN THE MAP AREA'])

disp(['THERE ARE ',num2str(length(i)),' SAMPLES IN THE DEFINED AREA'])
x=x(i);y=y(i);z=z(i);
if nargin==5;z1=z1(i);elseif nargin==5;z2=z2(i);end
% plot stations in new map area
figure(2),clf,
h1=axes('position',[.06 .06 .85 .85]);plot(x,y,'k.','MarkerSize',4),title('STATIONS IN SELECTED AREA')
set(gca,'fontsize',8)

disp(' DATA CAN BE INTERPOLATED FROM THE AREA SHOWN ')   
ans=input('WOULD YOU LIKE TO CREATE A MAPABLE AREA WITHIN THIS AREA? [y]/n','s');
ans
if isempty(ans); ans='y';end
ans
if ans=='y'|ans=='Y'
    disp(' DEFINE A MAP AREA WITHIN THE AREA SHOWN ')
    disp(' BE AWARE OF EXTARPOLATION ISSUES AT THE MAP EDGES')
  [xm,ym]=ginput(2);
  xmin=xm(1);
xmax=xm(2);
ymin=ym(1);
ymax=ym(2);
end

disp(' ')
   %
   % this block sets the x and y grid spacing 
   %
   
disp('      DEFINE THE GRID SPACINGS')
disp('GRID SPACINGS SHOULD BE SIMILAR TO SPACINGS IN THE xyz DATA')
disp('USE THE MAXIMUM, MINIMUM, AND MOST COMMON SPACINGS IN x AND y and VECTOR LENGTH')
disp('AS A GUIDE TO DETERMINE THE GRID SPACINGS')
disp(' ')
r=sqrt(x.^2+y.^2);
[nr,rr]=hist(abs(r(1:length(r)-1)-r(2:length(r)) ),30);
axes('position',[.69 .66 .2 .2]),bar(rr,nr);
set(gca,'fontsize',8)
xlabel('SEPARATION','fontname','times','fontsize',7)
ylabel('# OF STATIONS','fontname','times','fontsize',7);
title('HISTOGRAM','fontname','times','fontsize',8)
ir=find(nr==max(nr));
rdcommon=rr(ir);
rdmin=min(abs(r(2:length(r))-r(1:length(r)-1)));
rdmax=max(abs(r(2:length(r))-r(1:length(x)-1)));
disp(['THE MOST COMMON VECTOR SPACING IS ', num2str(rdcommon)])

[n,xx]=hist(abs(x(1:length(x)-1)-x(2:length(x)) ),10);
i=find(n==max(n));
xdcommon=xx(i);
xdmin=min(abs(x(2:length(x))-x(1:length(x)-1)));
xdmax=min(abs(x(2:length(x))-x(1:length(x)-1)));
disp(['THE MINIMUM, MAXIMUM AND MOST COMMON x AXIS SPACINGS ARE: ',num2str(xdmin),' ',num2str(xdmax),' ',num2str(xdcommon)]);
xd=input('WHAT x AXIS SPACING TO USE? ');
disp(' ')
[n,yy]=hist(abs(y(1:length(y)-1)-y(2:length(y)) ),10);
i=find(n==max(n));
ydcommon=yy(i);
ydmin=min(abs(y(2:length(y))-y(1:length(y)-1)));
ydmax=max(abs(y(2:length(y))-y(1:length(y)-1)));
disp([' THE MINIMUM, MAXIMUM AND MOST COMMON y AXIS SPACINGS ARE: ',num2str(ydmin),' ',num2str(ydmax),' ',num2str(ydcommon)]);
yd=input('WHAT y AXIS SPACING TO USE? ');
   %
   %

[X,Y]=meshgrid([xmin:xd:xmax],[ymin:yd:ymax]);
[ny,nx]=size(X);
disp(['THERE ARE ',num2str(nx),' GRID POINTS IN x DIRECTION'])
disp(['THERE ARE ',num2str(ny),' GRID POINTS IN y DIRECTION'])
   
Z=griddata(x,y,z,X,Y,method);
if nargin==5;Z1=griddata(x,y,z1,X,Y,method);elseif nargin==5;Z2=griddata(x,y,z2,X,Y,method);end
figure(1)
contour(X,Y,Z,20),hold on,title('CIRCLES ARE ORIGINAL STATIONS, DOTS ARE GRID POINTS ')
set(get(gca,'title'),'fontname','times','fontsize',14);
plot(x,y,'ko','MarkerSize',5)
plot(X,Y,'k.','MarkerSize',5)

if(any(any(isnan(Z)))),disp('WARNING NaN(s) ENCOUNTERED IN Z! '),end

[N_in,M_in]=size(Z); % original size
N=N_in;
M=M_in;
%%%%%%%%%%%%%%%%%%%%%%THIS BLOCK REPLACES NaNs WITH NEAREST NEIGHBOUR
  ZZ=reshape(Z,N*M,1); %
  XX=reshape(X,N*M,1); %
  YY=reshape(Y,N*M,1); %
  
  damnNaN=find(isnan(ZZ));
  NaNLocations=[XX(damnNaN) YY(damnNaN)];
  disp(' ')
S=[' WARNING ',sprintf('%6d',length(damnNaN)),' NaNs WERE FOUND'];
disp(S)
disp('LOCATIONS ARE AVAILABLE IN NaNLocations')
disp(' TO PLOT  plot(NanLocations(:,1),NanLocations(:,2),''r.'' ')