matlab_class/gcmfaces_calc/gcmfaces_interp_2d.m

function [val]=gcmfaces_interp_2d(fld,lon,lat,varargin);
%object:    linearly interpolate field to given positions
%inputs:    fld is a 2D gcmfaces field
%           lon,lat are position vectors
%optional:  doNearFill (1 by default) to use the nearest 
%               neighbor to complement linear interp
%outputs:   val is the vector of interpolated values
%
%pre-requisite: generate the delaunay triangulation using gcmfaces_bindata
%assumption: fld should show NaN for missing values

warning('off','MATLAB:dsearch:DeprecatedFunction');

global mygrid mytri myenv;

%inputs and pre-requisites
if ~isfield(myenv,'useDelaunayTri');
    myenv.useDelaunayTri=~isempty(which('DelaunayTri'));
end;

if isempty(mytri); 
    error('missing triangulation (mytri; from gcmfaces_bindata)');
end;

if nargin>3;
    doNearFill=varargin{1};
else;
    doNearFill=1;
end;

%switch longitude range to -180+180 or 0-360 according to grid
if max(mygrid.XC)<0;
    lon(find(lon>180))=lon(find(lon>180))-360;
end;

if max(mygrid.XC)>180;
    lon(find(lon<0))=lon(find(lon<0))+360;
end;

%do the actual interpolation    
if ~myenv.useDelaunayTri;%(code from old griddata.m)

    % Find the nearest triangle (t)
    x=convert2array(mygrid.XC); x=x(mytri.kk);
    y=convert2array(mygrid.YC); y=y(mytri.kk);
    VEC=convert2array(fld); VEC=VEC(mytri.kk);
    t = tsearch(x,y,mytri.TRI,lon',lat')';%the order of dims matters!!
    
    % Only keep the relevant triangles.
    out = find(isnan(t));
    if ~isempty(out), t(out) = ones(size(out)); end
    tri = mytri.TRI(t(:),:);
    
    % Compute Barycentric coordinates (w).  P. 78 in Watson.
    del = (x(tri(:,2))-x(tri(:,1))) .* (y(tri(:,3))-y(tri(:,1))) - ...
        (x(tri(:,3))-x(tri(:,1))) .* (y(tri(:,2))-y(tri(:,1)));
    w(:,3) = ((x(tri(:,1))-lon(:)).*(y(tri(:,2))-lat(:)) - ...
        (x(tri(:,2))-lon(:)).*(y(tri(:,1))-lat(:))) ./ del;
    w(:,2) = ((x(tri(:,3))-lon(:)).*(y(tri(:,1))-lat(:)) - ...
        (x(tri(:,1))-lon(:)).*(y(tri(:,3))-lat(:))) ./ del;
    w(:,1) = ((x(tri(:,2))-lon(:)).*(y(tri(:,3))-lat(:)) - ...
        (x(tri(:,3))-lon(:)).*(y(tri(:,2))-lat(:))) ./ del;
    
    val = sum(VEC(tri) .* w,2);    
    val=reshape(val,size(lon));
    
else;%(use TriScatteredInterp)
    
    %get interpolant
    VEC=convert2array(fld); VEC=VEC(mytri.kk);
    F = TriScatteredInterp(mytri.TRI, VEC);
    %do interpolation
    val=F(lon,lat);
    
end;

if doNearFill;%use the nearest neighbor to complement linear interp
    kk=gcmfaces_bindata(lon(:),lat(:));
    ARR=convert2array(fld);
    val2=ARR(kk);
    val(isnan(val))=val2(isnan(val));
end;


1	gforget	1.1	function [val]=gcmfaces_interp_2d(fld,lon,lat,varargin);
2			%object: linearly interpolate field to given positions
3			%inputs: fld is a 2D gcmfaces field
4			% lon,lat are position vectors
5			%optional: doNearFill (1 by default) to use the nearest
6			% neighbor to complement linear interp
7			%outputs: val is the vector of interpolated values
8			%
9			%pre-requisite: generate the delaunay triangulation using gcmfaces_bindata
10			%assumption: fld should show NaN for missing values
11
12			warning('off','MATLAB:dsearch:DeprecatedFunction');
13
14			global mygrid mytri myenv;
15
16			%inputs and pre-requisites
17			if ~isfield(myenv,'useDelaunayTri');
18			myenv.useDelaunayTri=~isempty(which('DelaunayTri'));
19			end;
20
21			if isempty(mytri);
22			error('missing triangulation (mytri; from gcmfaces_bindata)');
23			end;
24
25			if nargin>3;
26			doNearFill=varargin{1};
27			else;
28			doNearFill=1;
29			end;
30
31			%switch longitude range to -180+180 or 0-360 according to grid
32			if max(mygrid.XC)<0;
33			lon(find(lon>180))=lon(find(lon>180))-360;
34			end;
35
36			if max(mygrid.XC)>180;
37			lon(find(lon<0))=lon(find(lon<0))+360;
38			end;
39
40			%do the actual interpolation
41			if ~myenv.useDelaunayTri;%(code from old griddata.m)
42
43			% Find the nearest triangle (t)
44			x=convert2array(mygrid.XC); x=x(mytri.kk);
45			y=convert2array(mygrid.YC); y=y(mytri.kk);
46			VEC=convert2array(fld); VEC=VEC(mytri.kk);
47			t = tsearch(x,y,mytri.TRI,lon',lat')';%the order of dims matters!!
48
49			% Only keep the relevant triangles.
50			out = find(isnan(t));
51			if ~isempty(out), t(out) = ones(size(out)); end
52			tri = mytri.TRI(t(:),:);
53
54			% Compute Barycentric coordinates (w). P. 78 in Watson.
55			del = (x(tri(:,2))-x(tri(:,1))) .* (y(tri(:,3))-y(tri(:,1))) - ...
56			(x(tri(:,3))-x(tri(:,1))) .* (y(tri(:,2))-y(tri(:,1)));
57			w(:,3) = ((x(tri(:,1))-lon(:)).*(y(tri(:,2))-lat(:)) - ...
58			(x(tri(:,2))-lon(:)).*(y(tri(:,1))-lat(:))) ./ del;
59			w(:,2) = ((x(tri(:,3))-lon(:)).*(y(tri(:,1))-lat(:)) - ...
60			(x(tri(:,1))-lon(:)).*(y(tri(:,3))-lat(:))) ./ del;
61			w(:,1) = ((x(tri(:,2))-lon(:)).*(y(tri(:,3))-lat(:)) - ...
62			(x(tri(:,3))-lon(:)).*(y(tri(:,2))-lat(:))) ./ del;
63
64			val = sum(VEC(tri) .* w,2);
65			val=reshape(val,size(lon));
66
67			else;%(use TriScatteredInterp)
68
69			%get interpolant
70			VEC=convert2array(fld); VEC=VEC(mytri.kk);
71			F = TriScatteredInterp(mytri.TRI, VEC);
72			%do interpolation
73			val=F(lon,lat);
74
75			end;
76
77			if doNearFill;%use the nearest neighbor to complement linear interp
78			kk=gcmfaces_bindata(lon(:),lat(:));
79			ARR=convert2array(fld);
80			val2=ARR(kk);
81			val(isnan(val))=val2(isnan(val));
82			end;
83
84