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