matlab_class/gcmfaces_calc/gcmfaces_interp_2d.m

function [val]=gcmfaces_interp_2d(fld,lon,lat,method);
%[val]=GCMFACES_INTERP_2D(fld,lon,lat,method);
%    interpolates gcmfaces field (fld) to
%   a set of locations (lon, lat) using one of several methods:
%   'natural', 'linear' (default), 'nearest', 'mix', or 'polygons'.
%
%   This function relies on DelaunayTri except for the 'polygons' 
%   method (still under development). The 'mix' option is 
%   `natural' extended with 'nearest' when the input field 
%   has been land-masked with NaNs.
%
%Example:
%     lon=[-179.9:0.2:179.9]; lat=[-89.9:0.2:89.9];
%     [lat,lon] = meshgrid(lat,lon);
%     fld=mygrid.Depth.*mygrid.mskC(:,:,1);
%     [val]=gcmfaces_interp_2d(fld,lon,lat);
%     figureL; pcolor(lon,lat,val); shading flat;

gcmfaces_global;

%backward compatibility checks:
if ~isfield(myenv,'useDelaunayTri');
    myenv.useDelaunayTri=~isempty(which('DelaunayTri'));
end;

if isempty(whos('method')); method='linear'; end;

if ~ischar(method);
    error('fourth argument to gcmfaces_interp_2d is now ''method''');
end;

if ~myenv.useDelaunayTri;;
    warning('DelaunayTri is missing -> reverting to old tsearch method');

end;

%%use TriScatteredInterp
if strcmp(method,'linear')|strcmp(method,'nearest')|strcmp(method,'natural')|strcmp(method,'mix');

XC=convert2array(mygrid.XC);
YC=convert2array(mygrid.YC);
VEC=convert2array(fld);
val=NaN*lon;
for ii=1:3;
    if ii==1;
        myXC=XC; myXC(myXC<0)=myXC(myXC<0)+360;
        jj=find(lon>90);
    elseif ii==2;
        myXC=XC;
        jj=find(lon>=-90&lon<=90);
    else;
        myXC=XC; myXC(myXC>0)=myXC(myXC>0)-360;
        jj=find(lon<-90);
    end;
    kk=find(~isnan(myXC));
    TRI=DelaunayTri(myXC(kk),YC(kk));
    if strcmp(method,'mix');
        F = TriScatteredInterp(TRI, VEC(kk),'natural');
        tmp1=F(lon(jj),lat(jj));
        F = TriScatteredInterp(TRI, VEC(kk),'nearest');
        tmp2=F(lon(jj),lat(jj));
        tmp1(isnan(tmp1))=tmp2(isnan(tmp1));
        val(jj)=tmp1;
    else;
        F = TriScatteredInterp(TRI, VEC(kk),method);
        val(jj)=F(lon(jj),lat(jj));
    end;
end;

end;

%%old linear interpolation method:
if strcmp(method,'tsearch');

    % Generate triangulation
    global mytri;
    gcmfaces_bindata;

    %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;

    % 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(out)=NaN;
    val=reshape(val,size(lon));

end;%if strcmp(method,'tsearch');


%%use gcmfaces_interp_coeffs (under development):
if strcmp(method,'polygons');
    siz=size(lon);
    lon=lon(:);
    lat=lat(:);

    ni=30; nj=30;
    map_tile=gcmfaces_loc_tile(ni,nj);

    loc_tile=gcmfaces_loc_tile(ni,nj,lon,lat);
    loc_tile=loc_tile.tileNo;
    [loc_interp]=gcmfaces_interp_coeffs(lon,lat,ni,nj);

    fld=exch_T_N(fld);
    [tmp1,tmp2,n3,n4]=size(fld{1});
    arr=NaN*zeros(length(lon),n3,n4);

    for ii_tile=1:max(map_tile);
        ii=find(loc_tile==ii_tile&nansum(loc_interp.w,2));
        if ~isempty(ii);

            %1) determine face of current tile
            tmp1=1*(map_tile==ii_tile);
            tmp11=sum(sum(tmp1,1),2); tmp12=[];
            for kk=1:tmp11.nFaces; tmp12=[tmp12,tmp11{kk}]; end;
            ii_face=find(tmp12);
            %... and its index range within face ...
            tmp1=tmp1{ii_face};
            tmp11=sum(tmp1,2);
            iiMin=min(find(tmp11)); iiMax=max(find(tmp11));
            tmp11=sum(tmp1,1);
            jjMin=min(find(tmp11)); jjMax=max(find(tmp11));

            %2) determine points and weights for next step
            interp_i=1+loc_interp.i(ii,:);
            interp_j=1+loc_interp.j(ii,:);
            interp_k=(interp_j-1)*(ni+2)+interp_i;%index after reshape
            interp_w=loc_interp.w(ii,:);%bininear weights

            %3) get the tile array
            til=fld{ii_face}(iiMin:iiMax+2,jjMin:jjMax+2,:,:);
            msk=1*(~isnan(til));
            til(isnan(til))=0;

            %4) interpolate each month to profile locations
            tmp1=reshape(til,[(ni+2)*(nj+2) n3 n4]);
            tmp2=reshape(msk,[(ni+2)*(nj+2) n3 n4]);
            tmp11=tmp1(interp_k(:),:); tmp11=reshape(tmp11,[size(interp_k) n3 n4]);
            tmp22=tmp2(interp_k(:),:); tmp22=reshape(tmp22,[size(interp_k) n3 n4]);
            tmp1=squeeze(sum(tmp11,2));
            tmp2=squeeze(sum(tmp22,2));
            arr(ii,:,:,:)=tmp1./tmp2;

        end;%if ~isempty(ii);
    end;%for ii_tile=1:117;

    val=squeeze(reshape(arr,[siz n3 n4]));
end;

1	gforget	1.3	function [val]=gcmfaces_interp_2d(fld,lon,lat,method);
2			%[val]=GCMFACES_INTERP_2D(fld,lon,lat,method);
3			% interpolates gcmfaces field (fld) to
4			% a set of locations (lon, lat) using one of several methods:
5			% 'natural', 'linear' (default), 'nearest', 'mix', or 'polygons'.
6	gforget	1.1	%
7	gforget	1.3	% This function relies on DelaunayTri except for the 'polygons'
8			% method (still under development). The 'mix' option is
9			% `natural' extended with 'nearest' when the input field
10			% has been land-masked with NaNs.
11			%
12			%Example:
13			% lon=[-179.9:0.2:179.9]; lat=[-89.9:0.2:89.9];
14			% [lat,lon] = meshgrid(lat,lon);
15			% fld=mygrid.Depth.*mygrid.mskC(:,:,1);
16			% [val]=gcmfaces_interp_2d(fld,lon,lat);
17			% figureL; pcolor(lon,lat,val); shading flat;
18	gforget	1.1
19	gforget	1.3	gcmfaces_global;
20	gforget	1.1
21	gforget	1.3	%backward compatibility checks:
22	gforget	1.1	if ~isfield(myenv,'useDelaunayTri');
23			myenv.useDelaunayTri=~isempty(which('DelaunayTri'));
24			end;
25
26	gforget	1.3	if isempty(whos('method')); method='linear'; end;
27
28			if ~ischar(method);
29			error('fourth argument to gcmfaces_interp_2d is now ''method''');
30	gforget	1.1	end;
31
32	gforget	1.3	if ~myenv.useDelaunayTri;;
33			warning('DelaunayTri is missing -> reverting to old tsearch method');
34
35	gforget	1.1	end;
36
37	gforget	1.3	%%use TriScatteredInterp
38			if strcmp(method,'linear')\|strcmp(method,'nearest')\|strcmp(method,'natural')\|strcmp(method,'mix');
39
40			XC=convert2array(mygrid.XC);
41			YC=convert2array(mygrid.YC);
42			VEC=convert2array(fld);
43			val=NaN*lon;
44			for ii=1:3;
45			if ii==1;
46			myXC=XC; myXC(myXC<0)=myXC(myXC<0)+360;
47			jj=find(lon>90);
48			elseif ii==2;
49			myXC=XC;
50			jj=find(lon>=-90&lon<=90);
51			else;
52			myXC=XC; myXC(myXC>0)=myXC(myXC>0)-360;
53			jj=find(lon<-90);
54			end;
55			kk=find(~isnan(myXC));
56			TRI=DelaunayTri(myXC(kk),YC(kk));
57			if strcmp(method,'mix');
58			F = TriScatteredInterp(TRI, VEC(kk),'natural');
59			tmp1=F(lon(jj),lat(jj));
60			F = TriScatteredInterp(TRI, VEC(kk),'nearest');
61			tmp2=F(lon(jj),lat(jj));
62			tmp1(isnan(tmp1))=tmp2(isnan(tmp1));
63			val(jj)=tmp1;
64			else;
65			F = TriScatteredInterp(TRI, VEC(kk),method);
66			val(jj)=F(lon(jj),lat(jj));
67			end;
68	gforget	1.1	end;
69
70			end;
71
72	gforget	1.3	%%old linear interpolation method:
73			if strcmp(method,'tsearch');
74
75			% Generate triangulation
76			global mytri;
77			gcmfaces_bindata;
78
79			%switch longitude range to -180+180 or 0-360 according to grid
80			if max(mygrid.XC)<0;
81			lon(find(lon>180))=lon(find(lon>180))-360;
82			end;
83			if max(mygrid.XC)>180;
84			lon(find(lon<0))=lon(find(lon<0))+360;
85			end;
86	gforget	1.1
87			% Find the nearest triangle (t)
88			x=convert2array(mygrid.XC); x=x(mytri.kk);
89			y=convert2array(mygrid.YC); y=y(mytri.kk);
90			VEC=convert2array(fld); VEC=VEC(mytri.kk);
91			t = tsearch(x,y,mytri.TRI,lon',lat')';%the order of dims matters!!
92
93			% Only keep the relevant triangles.
94			out = find(isnan(t));
95			if ~isempty(out), t(out) = ones(size(out)); end
96			tri = mytri.TRI(t(:),:);
97
98			% Compute Barycentric coordinates (w). P. 78 in Watson.
99			del = (x(tri(:,2))-x(tri(:,1))) .* (y(tri(:,3))-y(tri(:,1))) - ...
100			(x(tri(:,3))-x(tri(:,1))) .* (y(tri(:,2))-y(tri(:,1)));
101			w(:,3) = ((x(tri(:,1))-lon(:)).*(y(tri(:,2))-lat(:)) - ...
102			(x(tri(:,2))-lon(:)).*(y(tri(:,1))-lat(:))) ./ del;
103			w(:,2) = ((x(tri(:,3))-lon(:)).*(y(tri(:,1))-lat(:)) - ...
104			(x(tri(:,1))-lon(:)).*(y(tri(:,3))-lat(:))) ./ del;
105			w(:,1) = ((x(tri(:,2))-lon(:)).*(y(tri(:,3))-lat(:)) - ...
106			(x(tri(:,3))-lon(:)).*(y(tri(:,2))-lat(:))) ./ del;
107
108			val = sum(VEC(tri) .* w,2);
109	gforget	1.2	val(out)=NaN;
110	gforget	1.1	val=reshape(val,size(lon));
111
112	gforget	1.3	end;%if strcmp(method,'tsearch');
113
114
115			%%use gcmfaces_interp_coeffs (under development):
116			if strcmp(method,'polygons');
117			siz=size(lon);
118			lon=lon(:);
119			lat=lat(:);
120
121			ni=30; nj=30;
122			map_tile=gcmfaces_loc_tile(ni,nj);
123
124			loc_tile=gcmfaces_loc_tile(ni,nj,lon,lat);
125			loc_tile=loc_tile.tileNo;
126			[loc_interp]=gcmfaces_interp_coeffs(lon,lat,ni,nj);
127
128			fld=exch_T_N(fld);
129			[tmp1,tmp2,n3,n4]=size(fld{1});
130			arr=NaN*zeros(length(lon),n3,n4);
131
132			for ii_tile=1:max(map_tile);
133			ii=find(loc_tile==ii_tile&nansum(loc_interp.w,2));
134			if ~isempty(ii);
135
136			%1) determine face of current tile
137			tmp1=1*(map_tile==ii_tile);
138			tmp11=sum(sum(tmp1,1),2); tmp12=[];
139			for kk=1:tmp11.nFaces; tmp12=[tmp12,tmp11{kk}]; end;
140			ii_face=find(tmp12);
141			%... and its index range within face ...
142			tmp1=tmp1{ii_face};
143			tmp11=sum(tmp1,2);
144			iiMin=min(find(tmp11)); iiMax=max(find(tmp11));
145			tmp11=sum(tmp1,1);
146			jjMin=min(find(tmp11)); jjMax=max(find(tmp11));
147
148			%2) determine points and weights for next step
149			interp_i=1+loc_interp.i(ii,:);
150			interp_j=1+loc_interp.j(ii,:);
151			interp_k=(interp_j-1)*(ni+2)+interp_i;%index after reshape
152			interp_w=loc_interp.w(ii,:);%bininear weights
153
154			%3) get the tile array
155			til=fld{ii_face}(iiMin:iiMax+2,jjMin:jjMax+2,:,:);
156			msk=1*(~isnan(til));
157			til(isnan(til))=0;
158
159			%4) interpolate each month to profile locations
160			tmp1=reshape(til,[(ni+2)*(nj+2) n3 n4]);
161			tmp2=reshape(msk,[(ni+2)*(nj+2) n3 n4]);
162			tmp11=tmp1(interp_k(:),:); tmp11=reshape(tmp11,[size(interp_k) n3 n4]);
163			tmp22=tmp2(interp_k(:),:); tmp22=reshape(tmp22,[size(interp_k) n3 n4]);
164			tmp1=squeeze(sum(tmp11,2));
165			tmp2=squeeze(sum(tmp22,2));
166			arr(ii,:,:,:)=tmp1./tmp2;
167
168			end;%if ~isempty(ii);
169			end;%for ii_tile=1:117;
170
171			val=squeeze(reshape(arr,[siz n3 n4]));
172	gforget	1.1	end;
173