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 a gcmfaces field (fld) to a set of locations (lon, lat) 
%   using one of several methods: 'polygons' (default), 'natural', 'linear',
%   'nearest', or 'mix'.
%
%   If instead fld is is an array (of the same size as lon and lat)
%   then the field is interpolated to the gcmfaces grid (in mygrid).
%
%   The 'polygons' method is a particular implementation of bilinear 
%   interpolation on the sphere. The other methods applies on DelaunayTri
%   in longitude/latitude coordinates. In particular, 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;

% note: using gcmfaces_bindata to determine nearest neighbors on the sphere
% can be problematic. To illustrate the problem: uncomment the old
% method (using gcmfaces_bindata) in gcmfaces_loc_tile.m, set myenv.verbose
% to 2 and run the example as explained above. Points at the edge
% of tile 59 provide an example -- see display generated by
% gcmfaces_interp_coeffs.m

gcmfaces_global;

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

if isempty(whos('method')); method='polygons'; 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');
    if ~isa(fld,'gcmfaces');
      error('old tsearch method only works with a gcmfaces input (fld)');
    end;
end;

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

if isa(fld,'gcmfaces');
  XC=convert2array(mygrid.XC);
  YC=convert2array(mygrid.YC);
  VEC=convert2array(fld);
else;
  XC=lon;
  YC=lat;
  VEC=fld;
  lon=convert2gcmfaces(mygrid.XC);
  lat=convert2gcmfaces(mygrid.YC);
end;

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;

if ~isa(fld,'gcmfaces');
  val=convert2gcmfaces(val);
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
            tmp0=repmat(interp_w,[1 1 n3 n4]);
            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(tmp0.*tmp11,2));
            tmp2=squeeze(sum(tmp0.*tmp22,2));
            arr(ii,:,:,:)=tmp1./tmp2;

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

    if sum(siz~=1)<=1;
      val=reshape(arr,[prod(siz) n3 n4]);
    else;
      val=reshape(arr,[siz n3 n4]);
    end;
%   val=squeeze(reshape(arr,[siz n3 n4]));
%   if size(val,1)~=size(lon,1); val=val'; end;
end;

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