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&~strcmp(method,'polygons');
    warning('DelaunayTri is missing -> reverting to old tsearch method');
    if ~isa(fld,'gcmfaces');
      error('old tsearch method only treats gcmfaces inputs (fld)');
    end;
end;

if ~isa(fld,'gcmfaces')&strcmp(method,'polygons');
    error('polygons method only treats gcmfaces inputs (fld)');
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;
    %ni=51; nj=51;
    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.7	if ~myenv.useDelaunayTri&~strcmp(method,'polygons');
44	gforget	1.3	warning('DelaunayTri is missing -> reverting to old tsearch method');
45	gforget	1.4	if ~isa(fld,'gcmfaces');
46	gforget	1.7	error('old tsearch method only treats gcmfaces inputs (fld)');
47	gforget	1.4	end;
48	gforget	1.1	end;
49
50	gforget	1.7	if ~isa(fld,'gcmfaces')&strcmp(method,'polygons');
51			error('polygons method only treats gcmfaces inputs (fld)');
52			end;
53
54	gforget	1.3	%%use TriScatteredInterp
55			if strcmp(method,'linear')\|strcmp(method,'nearest')\|strcmp(method,'natural')\|strcmp(method,'mix');
56
57	gforget	1.4	if isa(fld,'gcmfaces');
58			XC=convert2array(mygrid.XC);
59			YC=convert2array(mygrid.YC);
60			VEC=convert2array(fld);
61			else;
62			XC=lon;
63			YC=lat;
64			VEC=fld;
65			lon=convert2gcmfaces(mygrid.XC);
66			lat=convert2gcmfaces(mygrid.YC);
67			end;
68
69	gforget	1.3	val=NaN*lon;
70			for ii=1:3;
71			if ii==1;
72			myXC=XC; myXC(myXC<0)=myXC(myXC<0)+360;
73			jj=find(lon>90);
74			elseif ii==2;
75			myXC=XC;
76			jj=find(lon>=-90&lon<=90);
77			else;
78			myXC=XC; myXC(myXC>0)=myXC(myXC>0)-360;
79			jj=find(lon<-90);
80			end;
81			kk=find(~isnan(myXC));
82			TRI=DelaunayTri(myXC(kk),YC(kk));
83			if strcmp(method,'mix');
84			F = TriScatteredInterp(TRI, VEC(kk),'natural');
85			tmp1=F(lon(jj),lat(jj));
86			F = TriScatteredInterp(TRI, VEC(kk),'nearest');
87			tmp2=F(lon(jj),lat(jj));
88			tmp1(isnan(tmp1))=tmp2(isnan(tmp1));
89			val(jj)=tmp1;
90			else;
91			F = TriScatteredInterp(TRI, VEC(kk),method);
92			val(jj)=F(lon(jj),lat(jj));
93			end;
94	gforget	1.1	end;
95
96	gforget	1.4	if ~isa(fld,'gcmfaces');
97			val=convert2gcmfaces(val);
98			end;
99
100	gforget	1.1	end;
101
102	gforget	1.3	%%old linear interpolation method:
103			if strcmp(method,'tsearch');
104
105			% Generate triangulation
106			global mytri;
107			gcmfaces_bindata;
108
109			%switch longitude range to -180+180 or 0-360 according to grid
110			if max(mygrid.XC)<0;
111			lon(find(lon>180))=lon(find(lon>180))-360;
112			end;
113			if max(mygrid.XC)>180;
114			lon(find(lon<0))=lon(find(lon<0))+360;
115			end;
116	gforget	1.1
117			% Find the nearest triangle (t)
118			x=convert2array(mygrid.XC); x=x(mytri.kk);
119			y=convert2array(mygrid.YC); y=y(mytri.kk);
120			VEC=convert2array(fld); VEC=VEC(mytri.kk);
121			t = tsearch(x,y,mytri.TRI,lon',lat')';%the order of dims matters!!
122
123			% Only keep the relevant triangles.
124			out = find(isnan(t));
125			if ~isempty(out), t(out) = ones(size(out)); end
126			tri = mytri.TRI(t(:),:);
127
128			% Compute Barycentric coordinates (w). P. 78 in Watson.
129			del = (x(tri(:,2))-x(tri(:,1))) .* (y(tri(:,3))-y(tri(:,1))) - ...
130			(x(tri(:,3))-x(tri(:,1))) .* (y(tri(:,2))-y(tri(:,1)));
131			w(:,3) = ((x(tri(:,1))-lon(:)).*(y(tri(:,2))-lat(:)) - ...
132			(x(tri(:,2))-lon(:)).*(y(tri(:,1))-lat(:))) ./ del;
133			w(:,2) = ((x(tri(:,3))-lon(:)).*(y(tri(:,1))-lat(:)) - ...
134			(x(tri(:,1))-lon(:)).*(y(tri(:,3))-lat(:))) ./ del;
135			w(:,1) = ((x(tri(:,2))-lon(:)).*(y(tri(:,3))-lat(:)) - ...
136			(x(tri(:,3))-lon(:)).*(y(tri(:,2))-lat(:))) ./ del;
137
138			val = sum(VEC(tri) .* w,2);
139	gforget	1.2	val(out)=NaN;
140	gforget	1.1	val=reshape(val,size(lon));
141
142	gforget	1.3	end;%if strcmp(method,'tsearch');
143
144
145			%%use gcmfaces_interp_coeffs (under development):
146			if strcmp(method,'polygons');
147			siz=size(lon);
148			lon=lon(:);
149			lat=lat(:);
150
151			ni=30; nj=30;
152	gforget	1.7	%ni=51; nj=51;
153	gforget	1.3	map_tile=gcmfaces_loc_tile(ni,nj);
154
155			loc_tile=gcmfaces_loc_tile(ni,nj,lon,lat);
156			loc_tile=loc_tile.tileNo;
157			[loc_interp]=gcmfaces_interp_coeffs(lon,lat,ni,nj);
158
159			fld=exch_T_N(fld);
160			[tmp1,tmp2,n3,n4]=size(fld{1});
161			arr=NaN*zeros(length(lon),n3,n4);
162
163			for ii_tile=1:max(map_tile);
164			ii=find(loc_tile==ii_tile&nansum(loc_interp.w,2));
165			if ~isempty(ii);
166			%1) determine face of current tile
167			tmp1=1*(map_tile==ii_tile);
168			tmp11=sum(sum(tmp1,1),2); tmp12=[];
169			for kk=1:tmp11.nFaces; tmp12=[tmp12,tmp11{kk}]; end;
170			ii_face=find(tmp12);
171			%... and its index range within face ...
172			tmp1=tmp1{ii_face};
173			tmp11=sum(tmp1,2);
174			iiMin=min(find(tmp11)); iiMax=max(find(tmp11));
175			tmp11=sum(tmp1,1);
176			jjMin=min(find(tmp11)); jjMax=max(find(tmp11));
177
178			%2) determine points and weights for next step
179			interp_i=1+loc_interp.i(ii,:);
180			interp_j=1+loc_interp.j(ii,:);
181			interp_k=(interp_j-1)*(ni+2)+interp_i;%index after reshape
182			interp_w=loc_interp.w(ii,:);%bininear weights
183
184			%3) get the tile array
185			til=fld{ii_face}(iiMin:iiMax+2,jjMin:jjMax+2,:,:);
186			msk=1*(~isnan(til));
187			til(isnan(til))=0;
188
189			%4) interpolate each month to profile locations
190	gforget	1.5	tmp0=repmat(interp_w,[1 1 n3 n4]);
191	gforget	1.3	tmp1=reshape(til,[(ni+2)*(nj+2) n3 n4]);
192			tmp2=reshape(msk,[(ni+2)*(nj+2) n3 n4]);
193			tmp11=tmp1(interp_k(:),:); tmp11=reshape(tmp11,[size(interp_k) n3 n4]);
194			tmp22=tmp2(interp_k(:),:); tmp22=reshape(tmp22,[size(interp_k) n3 n4]);
195	gforget	1.5	tmp1=squeeze(sum(tmp0.*tmp11,2));
196			tmp2=squeeze(sum(tmp0.*tmp22,2));
197	gforget	1.3	arr(ii,:,:,:)=tmp1./tmp2;
198
199			end;%if ~isempty(ii);
200			end;%for ii_tile=1:117;
201
202	gforget	1.6	if sum(siz~=1)<=1;
203			val=reshape(arr,[prod(siz) n3 n4]);
204			else;
205			val=reshape(arr,[siz n3 n4]);
206			end;
207			% val=squeeze(reshape(arr,[siz n3 n4]));
208			% if size(val,1)~=size(lon,1); val=val'; end;
209	gforget	1.1	end;
210