matlab_class/gcmfaces_misc/gcmfaces_interp_coeffs.m

function [prof_interp,tile_corners]=gcmfaces_interp_coeffs(prof_lon,prof_lat,varargin);
%[prof_interp]=gcmfaces_interp_coeffs(prof_lon,prof_lat);
%object:    compute bilinear interpolation weights for prof_lon, prof_lat
%inputs:    prof_lon, prof_lat are column vectors
%(optional) ni,nj is the MITgcm tile size (2 numbers total)
%outputs:   prof_interp contains face and tile numbers,
%           indices within tiles (within 0:ni+1,0:nj+1)
%           and interpolation weights (between 0 and 1)
%           of the four neighboring grid points.
%(optional) tile_corners contains XC11,XCNINJ,YC11,YCNINJ (for MITprof)
%
%note: pathological cases (e.g. at edges) remain to be treated.
%example:
%prof=MITprof_load('ctd_feb2013.nc');
%[prof_interp]=gcmfaces_interp_coeffs(prof.prof_lon,prof.prof_lat);

gcmfaces_global;

doDisplay=0;
doVerbose=0; if myenv.verbose>=1; doVerbose=myenv.verbose; end;
%set doVerbose to display points no associated with a polygon

%set-up tile information (domain decomposition to ni,nj blocs)
if nargin<=2;
    ni=30; nj=30;
else;
    ni=varargin{1};
    nj=varargin{2};
end;

map_tile=gcmfaces_loc_tile(ni,nj);
loc_tile=gcmfaces_loc_tile(ni,nj,prof_lon,prof_lat);
prof_tile=loc_tile.tileNo;
list_tile=unique(prof_tile);

%initialize output:
prof_interp.face=NaN*prof_lon;
prof_interp.tile=NaN*prof_lon;
prof_interp.i=NaN*repmat(prof_lon,[1 4]);
prof_interp.j=NaN*repmat(prof_lon,[1 4]);
prof_interp.w=NaN*repmat(prof_lon,[1 4]);
%
tile_corners.XC11=NaN*prof_lon;
tile_corners.YC11=NaN*prof_lon;
tile_corners.XCNINJ=NaN*prof_lon;
tile_corners.YCNINJ=NaN*prof_lon;

%loop over tiles
for ii=1:length(list_tile);
    %1) determine face of current tile ...
    tmp1=1*(map_tile==list_tile(ii));
    tmp11=sum(sum(tmp1,1),2); tmp12=[];
    for ff=1:tmp11.nFaces; tmp12=[tmp12,tmp11{ff}]; end;
    iiFace=find(tmp12);
    %... and its index range within face ...
    tmp1=tmp1{iiFace};
    tmp11=sum(tmp1,2);
    iiMin=min(find(tmp11)); iiMax=max(find(tmp11));
    tmp11=sum(tmp1,1);
    jjMin=min(find(tmp11)); jjMax=max(find(tmp11));
    %... as well as the list of profiles in tile
    ii_prof=find(prof_tile==list_tile(ii));
    %tile corners
    XC11=mygrid.XC{iiFace}(iiMin,jjMin);
    YC11=mygrid.YC{iiFace}(iiMin,jjMin);
    XCNINJ=mygrid.XC{iiFace}(iiMax,jjMax);
    YCNINJ=mygrid.YC{iiFace}(iiMax,jjMax);

    clear tmp*;
    
    %2) stereographic projection to current tile center:
    ii0=floor((iiMin+iiMax)/2);
    jj0=floor((jjMin+jjMax)/2);
    XC0=mygrid.XC{iiFace}(ii0,jj0);
    YC0=mygrid.YC{iiFace}(ii0,jj0);
    %for grid locations:
    [xx,yy]=gcmfaces_stereoproj(XC0,YC0);
    %for profile locations
    [prof_x,prof_y]=gcmfaces_stereoproj(XC0,YC0,prof_lon,prof_lat);
    
    % nrm=sqrt(prof_x.^2+prof_y.^2);
    %ii_prof=find(nrm<tan(pi/4/2));%points inside of pi/4 cone
    
    %3) form array of grid cell quadrilaterals
    xxx=exch_T_N(xx); yyy=exch_T_N(yy);
    x_quad=[]; y_quad=[]; i_quad=[]; j_quad=[];
    for pp=1:4;
        switch pp;
            case 1; di=0; dj=0;
            case 2; di=1; dj=0;
            case 3; di=1; dj=1;
            case 4; di=0; dj=1;
        end;
        %note the shift in indices due to exchange above
        tmpx=xxx{iiFace}(iiMin+di:iiMax+1+di,jjMin+dj:jjMax+1+dj);
        tmpx=tmpx(:); x_quad=[x_quad tmpx];
        tmpy=yyy{iiFace}(iiMin+di:iiMax+1+di,jjMin+dj:jjMax+1+dj);
        tmpy=tmpy(:); y_quad=[y_quad tmpy];
        %
        tmpi=[0+di:iiMax-iiMin+1+di]'*ones(1,jjMax-jjMin+2);
        tmpi=tmpi(:); i_quad=[i_quad tmpi];
        tmpj=ones(jjMax-jjMin+2,1)*[0+dj:jjMax-jjMin+1+dj];
        tmpj=tmpj(:); j_quad=[j_quad tmpj];
    end;
    
    %4) associate profile locations with quadrilaterals
    [angsum]=gcmfaces_polygonangle(x_quad,y_quad,prof_x(ii_prof),prof_y(ii_prof));
    [II,JJ]=find(abs(angsum)>179);%+-360 for an interior point (+-180 for an edge point)
    if length(unique(JJ))~=length(JJ)&doVerbose;
        n0=num2str(length(JJ)-length(unique(JJ)));
        warning(['multiple polygons (' n0 ')']);
        %store indices of such instances for display
        [a,b]=hist(JJ,unique(JJ));
        KK=find(a>1);
        kk_prof=ii_prof(KK);
        kk_quad={};
        for kk=1:length(KK);
          kk_quad{kk}=II(find(JJ==KK(kk)));
        end
    else;
        kk_prof=[];
    end;
    if length(unique(JJ))<length(ii_prof)&doVerbose;
        n0=num2str(length(ii_prof)-length(unique(JJ)));
        n1=num2str(length(ii_prof));
        warning(['no polygon for ' n0 ' / ' n1]);
        %the following will then remove the corresponding profiles form ii_prof
    end;
    [C,IA,IC] = unique(JJ);
    %
    ii_prof0=ii_prof;
    ii_prof=ii_prof(C);%treated profiles
    jj_prof=setdiff(ii_prof0,ii_prof);%un-treated profiles
    %
    ii_quad=II(IA);

    if length(kk_prof)>0&doVerbose>=3;
      for kk=1:length(kk_prof);
        figureL;
        tmpx=x_quad(kk_quad{kk},[1:4 1])';
        tmpy=y_quad(kk_quad{kk},[1:4 1])';
        plot(tmpx,tmpy,'k.-','MarkerSize',36); hold on;
        plot(prof_x(kk_prof(kk)),prof_y(kk_prof(kk)),'r.','MarkerSize',36)
        aa=axis;
        aa(1:2)=aa(1:2)+abs(diff(aa(1:2)))*[-0.1 0.1];
        aa(3:4)=aa(3:4)+abs(diff(aa(3:4)))*[-0.1 0.1];
        axis(aa);
        keyboard;
      end;
    end;

    if length(jj_prof)>0&doVerbose>=2;
        figureL;
        %
        subplot(2,1,1);
        plot(prof_x(ii_prof),prof_y(ii_prof),'.');
        hold on; plot(x_quad(:),y_quad(:),'r.');
        plot(prof_x(jj_prof),prof_y(jj_prof),'k.','MarkerSize',36);
        %
        subplot(2,1,2);
        tmpx=convert2vector(mygrid.XC);
        tmpy=convert2vector(mygrid.YC);
        tmpi=convert2vector(map_tile);
        tmpi=find(tmpi==ii);
        plot(tmpx(:),tmpy(:),'r.'); hold on; 
        plot(tmpx(tmpi),tmpy(tmpi),'c.');
        plot(prof_lon(ii_prof),prof_lat(ii_prof),'.');
        plot(prof_lon(jj_prof),prof_lat(jj_prof),'k.','MarkerSize',36);
        %
        tmp1=prof_lon([ii_prof;jj_prof]);
        tmp2=prof_lat([ii_prof;jj_prof]);
        axis([min(tmp1) max(tmp1) min(tmp2) max(tmp2)]);
        %
        keyboard;
    end;
    
    if doDisplay;
        figureL;
        xx_tile=xxx{iiFace}(iiMin:iiMax+2,jjMin:jjMax+2);
        yy_tile=yyy{iiFace}(iiMin:iiMax+2,jjMin:jjMax+2);
        pcolor(xx_tile,yy_tile,sqrt(xx_tile.^2+yy_tile.^2)); hold on;
        cc=caxis; cc(2)=cc(2)*2; caxis(cc);
        plot(prof_x(ii_prof),prof_y(ii_prof),'r.','MarkerSize',20);
        plot(prof_x(jj_prof),prof_y(jj_prof),'k.','MarkerSize',60);
        axis([-0.6 0.6 -0.6 0.6]/6);
    end;
    
    if ~isempty(ii_prof);
        %5) determine bi-linear interpolation weights:
        px=x_quad(ii_quad,:);
        py=y_quad(ii_quad,:);
        ox=prof_x(ii_prof);
        oy=prof_y(ii_prof);
        [ow]=gcmfaces_quadmap(px,py,ox,oy);
        
        %to double check interpolation
        % pw=squeeze(ow);
        % oxInterp=sum(pw.*px,2);
        % oyInterp=sum(pw.*py,2);

        %round up coefficient to 4th digit (also to avoid slight negatives)
        test1=~isempty(find(ow(:)<-1e-5));
        if test1; error('interp weight < 0 -- something went wrong'); end;
        test1=~isempty(find(ow(:)>1+1e-5));
        if test1; error('interp weight >1 -- something went wrong'); end;
        %
        ow=1e-4*round(ow*1e4);
        sumw=repmat(sum(ow,3),[1 1 4]);
        ow=ow./sumw;
        
        %6) output interpolation specs:
        prof_interp.face(ii_prof,1)=iiFace*(1+0*ii_quad);
        prof_interp.tile(ii_prof,1)=list_tile(ii)*(1+0*ii_quad);
        prof_interp.i(ii_prof,:)=i_quad(ii_quad,:);
        prof_interp.j(ii_prof,:)=j_quad(ii_quad,:);
        prof_interp.w(ii_prof,:)=squeeze(ow);
        %
        tile_corners.XC11(ii_prof)=XC11;
        tile_corners.YC11(ii_prof)=YC11;
        tile_corners.XCNINJ(ii_prof)=XCNINJ;
        tile_corners.YCNINJ(ii_prof)=YCNINJ;
    end;
    
end;%for ii=1:length(list_tile);

% if doVerbose;
n1=sum(~isnan(prof_interp.face));
n2=sum(isnan(prof_interp.face));
fprintf(['interpolated points: ' num2str(n1) '\n']);
fprintf(['un-treated points:   ' num2str(n2) '\n']);
% end;

1	gforget	1.2	function [prof_interp,tile_corners]=gcmfaces_interp_coeffs(prof_lon,prof_lat,varargin);
2	gforget	1.1	%[prof_interp]=gcmfaces_interp_coeffs(prof_lon,prof_lat);
3			%object: compute bilinear interpolation weights for prof_lon, prof_lat
4			%inputs: prof_lon, prof_lat are column vectors
5			%(optional) ni,nj is the MITgcm tile size (2 numbers total)
6			%outputs: prof_interp contains face and tile numbers,
7			% indices within tiles (within 0:ni+1,0:nj+1)
8			% and interpolation weights (between 0 and 1)
9			% of the four neighboring grid points.
10	gforget	1.2	%(optional) tile_corners contains XC11,XCNINJ,YC11,YCNINJ (for MITprof)
11	gforget	1.1	%
12			%note: pathological cases (e.g. at edges) remain to be treated.
13			%example:
14			%prof=MITprof_load('ctd_feb2013.nc');
15			%[prof_interp]=gcmfaces_interp_coeffs(prof.prof_lon,prof.prof_lat);
16
17			gcmfaces_global;
18
19			doDisplay=0;
20	gforget	1.4	doVerbose=0; if myenv.verbose>=1; doVerbose=myenv.verbose; end;
21			%set doVerbose to display points no associated with a polygon
22	gforget	1.1
23			%set-up tile information (domain decomposition to ni,nj blocs)
24			if nargin<=2;
25			ni=30; nj=30;
26			else;
27			ni=varargin{1};
28			nj=varargin{2};
29			end;
30
31			map_tile=gcmfaces_loc_tile(ni,nj);
32			loc_tile=gcmfaces_loc_tile(ni,nj,prof_lon,prof_lat);
33			prof_tile=loc_tile.tileNo;
34			list_tile=unique(prof_tile);
35
36			%initialize output:
37			prof_interp.face=NaN*prof_lon;
38			prof_interp.tile=NaN*prof_lon;
39			prof_interp.i=NaN*repmat(prof_lon,[1 4]);
40			prof_interp.j=NaN*repmat(prof_lon,[1 4]);
41			prof_interp.w=NaN*repmat(prof_lon,[1 4]);
42	gforget	1.2	%
43			tile_corners.XC11=NaN*prof_lon;
44			tile_corners.YC11=NaN*prof_lon;
45			tile_corners.XCNINJ=NaN*prof_lon;
46			tile_corners.YCNINJ=NaN*prof_lon;
47	gforget	1.1
48			%loop over tiles
49			for ii=1:length(list_tile);
50			%1) determine face of current tile ...
51			tmp1=1*(map_tile==list_tile(ii));
52			tmp11=sum(sum(tmp1,1),2); tmp12=[];
53			for ff=1:tmp11.nFaces; tmp12=[tmp12,tmp11{ff}]; end;
54			iiFace=find(tmp12);
55			%... and its index range within face ...
56			tmp1=tmp1{iiFace};
57			tmp11=sum(tmp1,2);
58			iiMin=min(find(tmp11)); iiMax=max(find(tmp11));
59			tmp11=sum(tmp1,1);
60			jjMin=min(find(tmp11)); jjMax=max(find(tmp11));
61			%... as well as the list of profiles in tile
62			ii_prof=find(prof_tile==list_tile(ii));
63	gforget	1.2	%tile corners
64			XC11=mygrid.XC{iiFace}(iiMin,jjMin);
65			YC11=mygrid.YC{iiFace}(iiMin,jjMin);
66			XCNINJ=mygrid.XC{iiFace}(iiMax,jjMax);
67			YCNINJ=mygrid.YC{iiFace}(iiMax,jjMax);
68
69	gforget	1.1	clear tmp*;
70
71			%2) stereographic projection to current tile center:
72			ii0=floor((iiMin+iiMax)/2);
73			jj0=floor((jjMin+jjMax)/2);
74			XC0=mygrid.XC{iiFace}(ii0,jj0);
75			YC0=mygrid.YC{iiFace}(ii0,jj0);
76			%for grid locations:
77			[xx,yy]=gcmfaces_stereoproj(XC0,YC0);
78			%for profile locations
79			[prof_x,prof_y]=gcmfaces_stereoproj(XC0,YC0,prof_lon,prof_lat);
80
81			% nrm=sqrt(prof_x.^2+prof_y.^2);
82			%ii_prof=find(nrm<tan(pi/4/2));%points inside of pi/4 cone
83
84			%3) form array of grid cell quadrilaterals
85			xxx=exch_T_N(xx); yyy=exch_T_N(yy);
86			x_quad=[]; y_quad=[]; i_quad=[]; j_quad=[];
87			for pp=1:4;
88			switch pp;
89			case 1; di=0; dj=0;
90			case 2; di=1; dj=0;
91			case 3; di=1; dj=1;
92			case 4; di=0; dj=1;
93			end;
94			%note the shift in indices due to exchange above
95			tmpx=xxx{iiFace}(iiMin+di:iiMax+1+di,jjMin+dj:jjMax+1+dj);
96			tmpx=tmpx(:); x_quad=[x_quad tmpx];
97			tmpy=yyy{iiFace}(iiMin+di:iiMax+1+di,jjMin+dj:jjMax+1+dj);
98			tmpy=tmpy(:); y_quad=[y_quad tmpy];
99			%
100			tmpi=[0+di:iiMax-iiMin+1+di]'*ones(1,jjMax-jjMin+2);
101			tmpi=tmpi(:); i_quad=[i_quad tmpi];
102			tmpj=ones(jjMax-jjMin+2,1)*[0+dj:jjMax-jjMin+1+dj];
103			tmpj=tmpj(:); j_quad=[j_quad tmpj];
104			end;
105
106			%4) associate profile locations with quadrilaterals
107			[angsum]=gcmfaces_polygonangle(x_quad,y_quad,prof_x(ii_prof),prof_y(ii_prof));
108			[II,JJ]=find(abs(angsum)>179);%+-360 for an interior point (+-180 for an edge point)
109			if length(unique(JJ))~=length(JJ)&doVerbose;
110			n0=num2str(length(JJ)-length(unique(JJ)));
111			warning(['multiple polygons (' n0 ')']);
112	gforget	1.4	%store indices of such instances for display
113			[a,b]=hist(JJ,unique(JJ));
114			KK=find(a>1);
115			kk_prof=ii_prof(KK);
116			kk_quad={};
117			for kk=1:length(KK);
118			kk_quad{kk}=II(find(JJ==KK(kk)));
119			end
120			else;
121			kk_prof=[];
122	gforget	1.1	end;
123			if length(unique(JJ))<length(ii_prof)&doVerbose;
124			n0=num2str(length(ii_prof)-length(unique(JJ)));
125			n1=num2str(length(ii_prof));
126			warning(['no polygon for ' n0 ' / ' n1]);
127			%the following will then remove the corresponding profiles form ii_prof
128			end;
129			[C,IA,IC] = unique(JJ);
130			%
131			ii_prof0=ii_prof;
132			ii_prof=ii_prof(C);%treated profiles
133			jj_prof=setdiff(ii_prof0,ii_prof);%un-treated profiles
134			%
135			ii_quad=II(IA);
136	gforget	1.4
137			if length(kk_prof)>0&doVerbose>=3;
138			for kk=1:length(kk_prof);
139			figureL;
140			tmpx=x_quad(kk_quad{kk},[1:4 1])';
141			tmpy=y_quad(kk_quad{kk},[1:4 1])';
142			plot(tmpx,tmpy,'k.-','MarkerSize',36); hold on;
143			plot(prof_x(kk_prof(kk)),prof_y(kk_prof(kk)),'r.','MarkerSize',36)
144			aa=axis;
145			aa(1:2)=aa(1:2)+abs(diff(aa(1:2)))*[-0.1 0.1];
146			aa(3:4)=aa(3:4)+abs(diff(aa(3:4)))*[-0.1 0.1];
147			axis(aa);
148			keyboard;
149			end;
150			end;
151
152			if length(jj_prof)>0&doVerbose>=2;
153			figureL;
154			%
155			subplot(2,1,1);
156			plot(prof_x(ii_prof),prof_y(ii_prof),'.');
157			hold on; plot(x_quad(:),y_quad(:),'r.');
158			plot(prof_x(jj_prof),prof_y(jj_prof),'k.','MarkerSize',36);
159			%
160			subplot(2,1,2);
161			tmpx=convert2vector(mygrid.XC);
162			tmpy=convert2vector(mygrid.YC);
163			tmpi=convert2vector(map_tile);
164			tmpi=find(tmpi==ii);
165			plot(tmpx(:),tmpy(:),'r.'); hold on;
166			plot(tmpx(tmpi),tmpy(tmpi),'c.');
167			plot(prof_lon(ii_prof),prof_lat(ii_prof),'.');
168			plot(prof_lon(jj_prof),prof_lat(jj_prof),'k.','MarkerSize',36);
169			%
170			tmp1=prof_lon([ii_prof;jj_prof]);
171			tmp2=prof_lat([ii_prof;jj_prof]);
172			axis([min(tmp1) max(tmp1) min(tmp2) max(tmp2)]);
173			%
174			keyboard;
175			end;
176	gforget	1.1
177			if doDisplay;
178			figureL;
179			xx_tile=xxx{iiFace}(iiMin:iiMax+2,jjMin:jjMax+2);
180			yy_tile=yyy{iiFace}(iiMin:iiMax+2,jjMin:jjMax+2);
181			pcolor(xx_tile,yy_tile,sqrt(xx_tile.^2+yy_tile.^2)); hold on;
182			cc=caxis; cc(2)=cc(2)*2; caxis(cc);
183			plot(prof_x(ii_prof),prof_y(ii_prof),'r.','MarkerSize',20);
184			plot(prof_x(jj_prof),prof_y(jj_prof),'k.','MarkerSize',60);
185			axis([-0.6 0.6 -0.6 0.6]/6);
186			end;
187
188			if ~isempty(ii_prof);
189			%5) determine bi-linear interpolation weights:
190			px=x_quad(ii_quad,:);
191			py=y_quad(ii_quad,:);
192			ox=prof_x(ii_prof);
193			oy=prof_y(ii_prof);
194			[ow]=gcmfaces_quadmap(px,py,ox,oy);
195
196			%to double check interpolation
197			% pw=squeeze(ow);
198			% oxInterp=sum(pw.*px,2);
199			% oyInterp=sum(pw.*py,2);
200	gforget	1.3
201			%round up coefficient to 4th digit (also to avoid slight negatives)
202			test1=~isempty(find(ow(:)<-1e-5));
203			if test1; error('interp weight < 0 -- something went wrong'); end;
204			test1=~isempty(find(ow(:)>1+1e-5));
205			if test1; error('interp weight >1 -- something went wrong'); end;
206			%
207			ow=1e-4round(ow1e4);
208			sumw=repmat(sum(ow,3),[1 1 4]);
209			ow=ow./sumw;
210	gforget	1.1
211			%6) output interpolation specs:
212			prof_interp.face(ii_prof,1)=iiFace(1+0ii_quad);
213			prof_interp.tile(ii_prof,1)=list_tile(ii)(1+0ii_quad);
214			prof_interp.i(ii_prof,:)=i_quad(ii_quad,:);
215			prof_interp.j(ii_prof,:)=j_quad(ii_quad,:);
216			prof_interp.w(ii_prof,:)=squeeze(ow);
217	gforget	1.2	%
218			tile_corners.XC11(ii_prof)=XC11;
219			tile_corners.YC11(ii_prof)=YC11;
220			tile_corners.XCNINJ(ii_prof)=XCNINJ;
221			tile_corners.YCNINJ(ii_prof)=YCNINJ;
222	gforget	1.1	end;
223
224			end;%for ii=1:length(list_tile);
225
226			% if doVerbose;
227			n1=sum(~isnan(prof_interp.face));
228			n2=sum(isnan(prof_interp.face));
229			fprintf(['interpolated points: ' num2str(n1) '\n']);
230			fprintf(['un-treated points: ' num2str(n2) '\n']);
231			% end;
232