%call this function by
% output = griddata1(lat,lon,data,lat1,lon1,'nearest');
function [xq,yq,vq] = griddata1(varargin)
%GRIDDATA Interpolates scattered data - generally to produce gridded data
%   Vq = griddata(X,Y,V,Xq,Yq) fits a surface of the form V = F(X,Y) to the
%   scattered data in (X, Y, V). The coordinates of the data points are 
%   defined by the vectors (X,Y) and V defines the corresponding values.
%   griddata interpolates the surface F at the query points (Xq,Yq) and 
%   returns the values in Vq. The query points (Xq, Yq) generally represent
%   a grid obtained from NDGRID or MESHGRID, hence the name GRIDDATA.
%
%   Vq = griddata(X,Y,Z,V,Xq,Yq,Zq) fits a hyper-surface of the form
%   V = F(X,Y,Z) to the scattered data in (X, Y, Z, V). The coordinates of 
%   the data points are defined by the vectors (X,Y,Z) and V defines the 
%   corresponding values. griddata interpolates the surface F at the query
%   points (Xq,Yq,Zq) and returns the values in Vq.   
%
%   Vq = griddata(X,Y,V, xq, yq) where xq is a row vector and yq is a
%   column vector, expands (xq, yq) via [Xq, Yq] = meshgrid(xq,yq).
%   [Xq, Yq, Vq] = griddata(X,Y,V, xq, yq) returns the grid coordinates
%   arrays in addition.
%   Note: The syntax for implicit meshgrid expansion of (xq, yq) will be 
%   removed in a future release.
%
%   GRIDDATA(..., METHOD) where METHOD is one of
%       'nearest'   - Nearest neighbor interpolation
%       'linear'    - Linear interpolation (default)
%       'natural'   - Natural neighbor interpolation
%       'cubic'     - Cubic interpolation (2D only)
%       'v4'        - MATLAB 4 griddata method (2D only)
%   defines the interpolation method. The 'nearest' and 'linear' methods 
%   have discontinuities in the zero-th and first derivatives respectively, 
%   while the 'cubic' and 'v4' methods produce smooth surfaces.  All the 
%   methods except 'v4' are based on a Delaunay triangulation of the data.
%
%   Example 1:
%      % Interpolate a 2D scattered data set over a uniform grid
%      xy = -2.5 + 5*gallery('uniformdata',[200 2],0);
%      x = xy(:,1); y = xy(:,2);
%      v = x.*exp(-x.^2-y.^2);
%      [xq,yq] = meshgrid(-2:.2:2, -2:.2:2);
%      vq = griddata(x,y,v,xq,yq);
%      mesh(xq,yq,vq), hold on, plot3(x,y,v,'o'), hold off
%
%   Example 2:
%      % Interpolate a 3D data set over a grid in the x-y (z=0) plane
%      xyz = -1 + 2*gallery('uniformdata',[5000 3],0);
%      x = xyz(:,1); y = xyz(:,2); z = xyz(:,3);
%      v = x.^2 + y.^2 + z.^2;
%      d = -0.8:0.05:0.8;
%      [xq,yq,zq] = meshgrid(d,d,0);
%      vq = griddata(x,y,z,v,xq,yq,zq);
%      surf(xq,yq,vq);
%
%   See also scatteredInterpolant, GRIDDATAN, MESHGRID, NDGRID, DELAUNAY, 
%   INTERPN.

%   Copyright 1984-2015 The MathWorks, Inc. 

narginchk(5,9);

numarg = nargin;
method = 'linear';
if iscell(varargin{numarg})
    error(message('MATLAB:griddata:DeprecatedOptions'));
elseif ischar(varargin{numarg}) || (isstring(varargin{numarg}) && isscalar(varargin{numarg}))
    method = varargin{numarg};
    method = lower(method);
    numarg = numarg-1;
end

if ~any(strcmp(method, {'nearest', 'linear', 'natural', 'cubic', 'v4'}))
    error(message('MATLAB:griddata:UnknownMethod'));
end

if  numarg == 5
    numdims = 2;
elseif numarg == 7
    numdims = 3;
else
    error(message('MATLAB:griddata:InvalidNumInputArgs'));
end

for i=1:(2*numdims)+1
    if (i ~= (numdims+1) && ~isreal(varargin{i}) )
        error(message('MATLAB:griddata:InvalidCoordsComplex'));
    elseif ~isnumeric(varargin{i})
        error(message('MATLAB:griddata:InvalidInputArgs'));
    end
end


for i=1:numarg
    if ndims(varargin{i}) > numdims
        error(message('MATLAB:griddata:HigherDimArray'));
    elseif ( issparse(varargin{i}) )
        error(message('MATLAB:griddata:InvalidDataSparse'));
    end
end

if  numarg == 5
    % potentially 2D validate the data
    % The xyzchk generates a meshgrid - support for this will be removed
    % in a future release.
    x = varargin{1};
    y = varargin{2};
    v = varargin{3};
    xq = varargin{4};
    yq = varargin{5};
    [msg,x,y,~,xq,yq] = xyzchk(x,y,v,xq,yq);
    if ~isempty(msg), error(message(msg.identifier)); end
    inputargs = {x,y,v,xq,yq};
elseif numarg == 7
    % Potentially 3D, check support for the method
    inputargs = varargin;
    if strcmp(method, 'cubic')
        error(message('MATLAB:griddata:CubicMethod3D'));
    elseif strcmp(method, 'v4')
        error(message('MATLAB:griddata:V4Method3D'));
    end
end

switch method
    case 'nearest'
        %change 
        vq = useScatteredInterp(inputargs, numarg, method, 'none'); %change nearest->none
    case {'linear', 'natural'}
        vq = useScatteredInterp(inputargs, numarg, method, 'none');
    case 'cubic'
        vq = cubic(x,y,v,xq,yq);
    case 'v4'
        vq = gdatav4(x,y,v,xq,yq);
end

if nargout<=1, xq = vq; end

%------------------------------------------------------------

function [x, y, v] = mergepoints2D(x,y,v)

% Sort x and y so duplicate points can be averaged

%Need x,y and z to be column vectors
sz = numel(x);
x = reshape(x,sz,1);
y = reshape(y,sz,1);
v = reshape(v,sz,1);
myepsx = eps(0.5 * (max(x) - min(x)))^(1/3);
myepsy = eps(0.5 * (max(y) - min(y)))^(1/3);


% look for x, y points that are indentical (within a tolerance)
% average out the values for these points
if isreal(v)
    xyv = builtin('_mergesimpts', [y, x, v], [myepsy, myepsx, Inf], 'average');
    x = xyv(:,2);
    y = xyv(:,1);
    v = xyv(:,3);
else
    % if z is imaginary split out the real and imaginary parts
    xyv = builtin('_mergesimpts', [y, x, real(v), imag(v)], ...
        [myepsy, myepsx, Inf, Inf], 'average');
    x = xyv(:,2);
    y = xyv(:,1);
    % re-combine the real and imaginary parts
    v = xyv(:,3) + 1i*xyv(:,4);
end
% give a warning if some of the points were duplicates (and averaged out)
if sz>numel(x)
    warning(message('MATLAB:griddata:DuplicateDataPoints'));
end

%------------------------------------------------------------

function vq = useScatteredInterp(inargs, numarg, method, emeth)

% Reference (nearest, linear):
%    David F. Watson, "Contouring: A guide to the analysis and display
%       of spacial data", Pergamon, 1994.
%
% Reference (natural):
%    Sibson, R. (1981). "A brief description of natural neighbor
%       interpolation (Chapter 2)".  In V. Barnett. Interpreting
%       Multivariate Data.  Chichester: John Wiley. pp. 21--36.

if numarg == 5
    F = scatteredInterpolant(inargs{1}(:),inargs{2}(:),inargs{3}(:), ...
        method, emeth);
    vq = F(inargs{4},inargs{5});
elseif numarg == 7
    F = scatteredInterpolant(inargs{1}(:),inargs{2}(:),inargs{3}(:), ...
        inargs{4}(:), method, emeth);
    vq = F(inargs{5},inargs{6},inargs{7});
end

%------------------------------------------------------------

function vq = cubic(x,y,v,xq,yq)
%TRIANGLE Triangle-based cubic interpolation

%   Reference: T. Y. Yang, "Finite Element Structural Analysis",
%   Prentice Hall, 1986.  pp. 446-449.
%
%   Reference: David F. Watson, "Contouring: A guide
%   to the analysis and display of spacial data", Pergamon, 1994.

% Triangulate the data

[x, y, v] = mergepoints2D(x,y,v);

dt = delaunayTriangulation(x,y);
scopedWarnOff = warning('off', 'MATLAB:triangulation:EmptyTri2DWarnId');
restoreWarnOff = onCleanup(@()warning(scopedWarnOff));
dtt = dt.ConnectivityList;
if isempty(dtt)
    warning(message('MATLAB:griddata:EmptyTriangulation'));
    vq = [];
    return
end

tri = dt.ConnectivityList;
% Find the enclosing triangle (t)
siz = size(xq);
t = dt.pointLocation(xq(:),yq(:));
t = reshape(t,siz);

if(isreal(v))
    vq = cubicmx(x,y,v,xq,yq,tri,t);
else
    vre = real(v);
    vim = imag(v);
    vqre = cubicmx(x,y,vre,xq,yq,tri,t);
    vqim = cubicmx(x,y,vim,xq,yq,tri,t);
    vq = complex(vqre,vqim);
end

%------------------------------------------------------------

function vq = gdatav4(x,y,v,xq,yq)
%GDATAV4 MATLAB 4 GRIDDATA interpolation

%   Reference:  David T. Sandwell, Biharmonic spline
%   interpolation of GEOS-3 and SEASAT altimeter
%   data, Geophysical Research Letters, 2, 139-142,
%   1987.  Describes interpolation using value or
%   gradient of value in any dimension.

[x, y, v] = mergepoints2D(x,y,v);

xy = x(:) + 1i*y(:);

% Determine distances between points
d = abs(xy - xy.');

% Determine weights for interpolation
g = (d.^2) .* (log(d)-1);   % Green's function.
% Fixup value of Green's function along diagonal
g(1:size(d,1)+1:end) = 0;
weights = g \ v(:);

[m,n] = size(xq);
vq = zeros(size(xq));
xy = xy.';

% Evaluate at requested points (xq,yq).  Loop to save memory.
for i=1:m
    for j=1:n
        d = abs(xq(i,j) + 1i*yq(i,j) - xy);
        g = (d.^2) .* (log(d)-1);   % Green's function.
        % Value of Green's function at zero
        g(d==0) = 0;
        vq(i,j) = g * weights;        
    end
end