| 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.8 |
doVerbose=0; %if myenv.verbose>=1; doVerbose=myenv.verbose; end; |
| 21 |
gforget |
1.5 |
%set doVerbose to display points that could not be triangulated (if any) |
| 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 |
gforget |
1.7 |
%point indices in vector format |
| 37 |
|
|
map_vpi=convert2vector(mygrid.XC,'new'); |
| 38 |
|
|
map_vpi=convert2vector([1:length(map_vpi)]','new'); |
| 39 |
|
|
|
| 40 |
gforget |
1.1 |
%initialize output: |
| 41 |
gforget |
1.6 |
prof_interp.point=loc_tile.point; |
| 42 |
gforget |
1.1 |
prof_interp.face=NaN*prof_lon; |
| 43 |
|
|
prof_interp.tile=NaN*prof_lon; |
| 44 |
|
|
prof_interp.i=NaN*repmat(prof_lon,[1 4]); |
| 45 |
|
|
prof_interp.j=NaN*repmat(prof_lon,[1 4]); |
| 46 |
|
|
prof_interp.w=NaN*repmat(prof_lon,[1 4]); |
| 47 |
gforget |
1.6 |
prof_interp.XC=NaN*repmat(prof_lon,[1 4]); |
| 48 |
|
|
prof_interp.YC=NaN*repmat(prof_lon,[1 4]); |
| 49 |
gforget |
1.7 |
prof_interp.vpi=NaN*repmat(prof_lon,[1 4]); |
| 50 |
gforget |
1.2 |
% |
| 51 |
|
|
tile_corners.XC11=NaN*prof_lon; |
| 52 |
|
|
tile_corners.YC11=NaN*prof_lon; |
| 53 |
|
|
tile_corners.XCNINJ=NaN*prof_lon; |
| 54 |
|
|
tile_corners.YCNINJ=NaN*prof_lon; |
| 55 |
gforget |
1.1 |
|
| 56 |
|
|
%loop over tiles |
| 57 |
|
|
for ii=1:length(list_tile); |
| 58 |
|
|
%1) determine face of current tile ... |
| 59 |
|
|
tmp1=1*(map_tile==list_tile(ii)); |
| 60 |
|
|
tmp11=sum(sum(tmp1,1),2); tmp12=[]; |
| 61 |
|
|
for ff=1:tmp11.nFaces; tmp12=[tmp12,tmp11{ff}]; end; |
| 62 |
|
|
iiFace=find(tmp12); |
| 63 |
|
|
%... and its index range within face ... |
| 64 |
|
|
tmp1=tmp1{iiFace}; |
| 65 |
|
|
tmp11=sum(tmp1,2); |
| 66 |
|
|
iiMin=min(find(tmp11)); iiMax=max(find(tmp11)); |
| 67 |
|
|
tmp11=sum(tmp1,1); |
| 68 |
|
|
jjMin=min(find(tmp11)); jjMax=max(find(tmp11)); |
| 69 |
|
|
%... as well as the list of profiles in tile |
| 70 |
|
|
ii_prof=find(prof_tile==list_tile(ii)); |
| 71 |
gforget |
1.2 |
%tile corners |
| 72 |
|
|
XC11=mygrid.XC{iiFace}(iiMin,jjMin); |
| 73 |
|
|
YC11=mygrid.YC{iiFace}(iiMin,jjMin); |
| 74 |
|
|
XCNINJ=mygrid.XC{iiFace}(iiMax,jjMax); |
| 75 |
|
|
YCNINJ=mygrid.YC{iiFace}(iiMax,jjMax); |
| 76 |
|
|
|
| 77 |
gforget |
1.1 |
clear tmp*; |
| 78 |
|
|
|
| 79 |
|
|
%2) stereographic projection to current tile center: |
| 80 |
|
|
ii0=floor((iiMin+iiMax)/2); |
| 81 |
|
|
jj0=floor((jjMin+jjMax)/2); |
| 82 |
|
|
XC0=mygrid.XC{iiFace}(ii0,jj0); |
| 83 |
|
|
YC0=mygrid.YC{iiFace}(ii0,jj0); |
| 84 |
|
|
%for grid locations: |
| 85 |
|
|
[xx,yy]=gcmfaces_stereoproj(XC0,YC0); |
| 86 |
|
|
%for profile locations |
| 87 |
|
|
[prof_x,prof_y]=gcmfaces_stereoproj(XC0,YC0,prof_lon,prof_lat); |
| 88 |
|
|
|
| 89 |
|
|
% nrm=sqrt(prof_x.^2+prof_y.^2); |
| 90 |
|
|
%ii_prof=find(nrm<tan(pi/4/2));%points inside of pi/4 cone |
| 91 |
|
|
|
| 92 |
|
|
%3) form array of grid cell quadrilaterals |
| 93 |
|
|
xxx=exch_T_N(xx); yyy=exch_T_N(yy); |
| 94 |
gforget |
1.7 |
xc=exch_T_N(mygrid.XC); yc=exch_T_N(mygrid.YC); vpi=exch_T_N(map_vpi); |
| 95 |
|
|
x_quad=[]; y_quad=[]; xc_quad=[]; yc_quad=[]; vpi_quad=[]; i_quad=[]; j_quad=[]; |
| 96 |
gforget |
1.1 |
for pp=1:4; |
| 97 |
|
|
switch pp; |
| 98 |
|
|
case 1; di=0; dj=0; |
| 99 |
|
|
case 2; di=1; dj=0; |
| 100 |
|
|
case 3; di=1; dj=1; |
| 101 |
|
|
case 4; di=0; dj=1; |
| 102 |
|
|
end; |
| 103 |
|
|
%note the shift in indices due to exchange above |
| 104 |
|
|
tmpx=xxx{iiFace}(iiMin+di:iiMax+1+di,jjMin+dj:jjMax+1+dj); |
| 105 |
|
|
tmpx=tmpx(:); x_quad=[x_quad tmpx]; |
| 106 |
|
|
tmpy=yyy{iiFace}(iiMin+di:iiMax+1+di,jjMin+dj:jjMax+1+dj); |
| 107 |
|
|
tmpy=tmpy(:); y_quad=[y_quad tmpy]; |
| 108 |
|
|
% |
| 109 |
gforget |
1.6 |
tmpx=xc{iiFace}(iiMin+di:iiMax+1+di,jjMin+dj:jjMax+1+dj); |
| 110 |
|
|
tmpx=tmpx(:); xc_quad=[xc_quad tmpx]; |
| 111 |
|
|
tmpy=yc{iiFace}(iiMin+di:iiMax+1+di,jjMin+dj:jjMax+1+dj); |
| 112 |
|
|
tmpy=tmpy(:); yc_quad=[yc_quad tmpy]; |
| 113 |
gforget |
1.7 |
tmpvpi=vpi{iiFace}(iiMin+di:iiMax+1+di,jjMin+dj:jjMax+1+dj); |
| 114 |
|
|
tmpvpi=tmpvpi(:); vpi_quad=[vpi_quad tmpvpi]; |
| 115 |
gforget |
1.6 |
% |
| 116 |
gforget |
1.1 |
tmpi=[0+di:iiMax-iiMin+1+di]'*ones(1,jjMax-jjMin+2); |
| 117 |
|
|
tmpi=tmpi(:); i_quad=[i_quad tmpi]; |
| 118 |
|
|
tmpj=ones(jjMax-jjMin+2,1)*[0+dj:jjMax-jjMin+1+dj]; |
| 119 |
|
|
tmpj=tmpj(:); j_quad=[j_quad tmpj]; |
| 120 |
|
|
end; |
| 121 |
|
|
|
| 122 |
|
|
%4) associate profile locations with quadrilaterals |
| 123 |
|
|
[angsum]=gcmfaces_polygonangle(x_quad,y_quad,prof_x(ii_prof),prof_y(ii_prof)); |
| 124 |
|
|
[II,JJ]=find(abs(angsum)>179);%+-360 for an interior point (+-180 for an edge point) |
| 125 |
|
|
if length(unique(JJ))~=length(JJ)&doVerbose; |
| 126 |
|
|
n0=num2str(length(JJ)-length(unique(JJ))); |
| 127 |
|
|
warning(['multiple polygons (' n0 ')']); |
| 128 |
gforget |
1.4 |
%store indices of such instances for display |
| 129 |
|
|
[a,b]=hist(JJ,unique(JJ)); |
| 130 |
|
|
KK=find(a>1); |
| 131 |
|
|
kk_prof=ii_prof(KK); |
| 132 |
|
|
kk_quad={}; |
| 133 |
|
|
for kk=1:length(KK); |
| 134 |
|
|
kk_quad{kk}=II(find(JJ==KK(kk))); |
| 135 |
|
|
end |
| 136 |
|
|
else; |
| 137 |
|
|
kk_prof=[]; |
| 138 |
gforget |
1.1 |
end; |
| 139 |
|
|
if length(unique(JJ))<length(ii_prof)&doVerbose; |
| 140 |
|
|
n0=num2str(length(ii_prof)-length(unique(JJ))); |
| 141 |
|
|
n1=num2str(length(ii_prof)); |
| 142 |
|
|
warning(['no polygon for ' n0 ' / ' n1]); |
| 143 |
|
|
%the following will then remove the corresponding profiles form ii_prof |
| 144 |
|
|
end; |
| 145 |
|
|
[C,IA,IC] = unique(JJ); |
| 146 |
|
|
% |
| 147 |
|
|
ii_prof0=ii_prof; |
| 148 |
|
|
ii_prof=ii_prof(C);%treated profiles |
| 149 |
|
|
jj_prof=setdiff(ii_prof0,ii_prof);%un-treated profiles |
| 150 |
|
|
% |
| 151 |
|
|
ii_quad=II(IA); |
| 152 |
gforget |
1.4 |
|
| 153 |
|
|
if length(kk_prof)>0&doVerbose>=3; |
| 154 |
|
|
for kk=1:length(kk_prof); |
| 155 |
|
|
figureL; |
| 156 |
|
|
tmpx=x_quad(kk_quad{kk},[1:4 1])'; |
| 157 |
|
|
tmpy=y_quad(kk_quad{kk},[1:4 1])'; |
| 158 |
|
|
plot(tmpx,tmpy,'k.-','MarkerSize',36); hold on; |
| 159 |
|
|
plot(prof_x(kk_prof(kk)),prof_y(kk_prof(kk)),'r.','MarkerSize',36) |
| 160 |
|
|
aa=axis; |
| 161 |
|
|
aa(1:2)=aa(1:2)+abs(diff(aa(1:2)))*[-0.1 0.1]; |
| 162 |
|
|
aa(3:4)=aa(3:4)+abs(diff(aa(3:4)))*[-0.1 0.1]; |
| 163 |
|
|
axis(aa); |
| 164 |
|
|
keyboard; |
| 165 |
|
|
end; |
| 166 |
|
|
end; |
| 167 |
|
|
|
| 168 |
|
|
if length(jj_prof)>0&doVerbose>=2; |
| 169 |
|
|
figureL; |
| 170 |
|
|
% |
| 171 |
|
|
subplot(2,1,1); |
| 172 |
|
|
plot(prof_x(ii_prof),prof_y(ii_prof),'.'); |
| 173 |
|
|
hold on; plot(x_quad(:),y_quad(:),'r.'); |
| 174 |
|
|
plot(prof_x(jj_prof),prof_y(jj_prof),'k.','MarkerSize',36); |
| 175 |
|
|
% |
| 176 |
|
|
subplot(2,1,2); |
| 177 |
|
|
tmpx=convert2vector(mygrid.XC); |
| 178 |
|
|
tmpy=convert2vector(mygrid.YC); |
| 179 |
|
|
tmpi=convert2vector(map_tile); |
| 180 |
|
|
tmpi=find(tmpi==ii); |
| 181 |
|
|
plot(tmpx(:),tmpy(:),'r.'); hold on; |
| 182 |
|
|
plot(tmpx(tmpi),tmpy(tmpi),'c.'); |
| 183 |
|
|
plot(prof_lon(ii_prof),prof_lat(ii_prof),'.'); |
| 184 |
|
|
plot(prof_lon(jj_prof),prof_lat(jj_prof),'k.','MarkerSize',36); |
| 185 |
|
|
% |
| 186 |
|
|
tmp1=prof_lon([ii_prof;jj_prof]); |
| 187 |
|
|
tmp2=prof_lat([ii_prof;jj_prof]); |
| 188 |
|
|
axis([min(tmp1) max(tmp1) min(tmp2) max(tmp2)]); |
| 189 |
|
|
% |
| 190 |
|
|
keyboard; |
| 191 |
|
|
end; |
| 192 |
gforget |
1.1 |
|
| 193 |
|
|
if doDisplay; |
| 194 |
|
|
figureL; |
| 195 |
|
|
xx_tile=xxx{iiFace}(iiMin:iiMax+2,jjMin:jjMax+2); |
| 196 |
|
|
yy_tile=yyy{iiFace}(iiMin:iiMax+2,jjMin:jjMax+2); |
| 197 |
|
|
pcolor(xx_tile,yy_tile,sqrt(xx_tile.^2+yy_tile.^2)); hold on; |
| 198 |
|
|
cc=caxis; cc(2)=cc(2)*2; caxis(cc); |
| 199 |
|
|
plot(prof_x(ii_prof),prof_y(ii_prof),'r.','MarkerSize',20); |
| 200 |
|
|
plot(prof_x(jj_prof),prof_y(jj_prof),'k.','MarkerSize',60); |
| 201 |
|
|
axis([-0.6 0.6 -0.6 0.6]/6); |
| 202 |
|
|
end; |
| 203 |
|
|
|
| 204 |
|
|
if ~isempty(ii_prof); |
| 205 |
|
|
%5) determine bi-linear interpolation weights: |
| 206 |
|
|
px=x_quad(ii_quad,:); |
| 207 |
|
|
py=y_quad(ii_quad,:); |
| 208 |
|
|
ox=prof_x(ii_prof); |
| 209 |
|
|
oy=prof_y(ii_prof); |
| 210 |
|
|
[ow]=gcmfaces_quadmap(px,py,ox,oy); |
| 211 |
|
|
|
| 212 |
|
|
%to double check interpolation |
| 213 |
|
|
% pw=squeeze(ow); |
| 214 |
|
|
% oxInterp=sum(pw.*px,2); |
| 215 |
|
|
% oyInterp=sum(pw.*py,2); |
| 216 |
gforget |
1.3 |
|
| 217 |
|
|
%round up coefficient to 4th digit (also to avoid slight negatives) |
| 218 |
|
|
test1=~isempty(find(ow(:)<-1e-5)); |
| 219 |
|
|
if test1; error('interp weight < 0 -- something went wrong'); end; |
| 220 |
|
|
test1=~isempty(find(ow(:)>1+1e-5)); |
| 221 |
|
|
if test1; error('interp weight >1 -- something went wrong'); end; |
| 222 |
|
|
% |
| 223 |
|
|
ow=1e-4*round(ow*1e4); |
| 224 |
|
|
sumw=repmat(sum(ow,3),[1 1 4]); |
| 225 |
|
|
ow=ow./sumw; |
| 226 |
gforget |
1.1 |
|
| 227 |
|
|
%6) output interpolation specs: |
| 228 |
|
|
prof_interp.face(ii_prof,1)=iiFace*(1+0*ii_quad); |
| 229 |
|
|
prof_interp.tile(ii_prof,1)=list_tile(ii)*(1+0*ii_quad); |
| 230 |
|
|
prof_interp.i(ii_prof,:)=i_quad(ii_quad,:); |
| 231 |
|
|
prof_interp.j(ii_prof,:)=j_quad(ii_quad,:); |
| 232 |
|
|
prof_interp.w(ii_prof,:)=squeeze(ow); |
| 233 |
gforget |
1.6 |
prof_interp.XC(ii_prof,:)=xc_quad(ii_quad,:); |
| 234 |
|
|
prof_interp.YC(ii_prof,:)=yc_quad(ii_quad,:); |
| 235 |
gforget |
1.7 |
prof_interp.vpi(ii_prof,:)=vpi_quad(ii_quad,:); |
| 236 |
gforget |
1.2 |
% |
| 237 |
|
|
tile_corners.XC11(ii_prof)=XC11; |
| 238 |
|
|
tile_corners.YC11(ii_prof)=YC11; |
| 239 |
|
|
tile_corners.XCNINJ(ii_prof)=XCNINJ; |
| 240 |
|
|
tile_corners.YCNINJ(ii_prof)=YCNINJ; |
| 241 |
gforget |
1.1 |
end; |
| 242 |
gforget |
1.7 |
|
| 243 |
gforget |
1.1 |
end;%for ii=1:length(list_tile); |
| 244 |
|
|
|
| 245 |
gforget |
1.7 |
%now format interpolation as sparse matrix and delete vpi: |
| 246 |
|
|
sizin=convert2vector(mygrid.XC,'new'); sizin=length(sizin); |
| 247 |
|
|
sizout=size(prof_interp.vpi); |
| 248 |
|
|
tmpi=[1:sizout(1)]'*ones(1,sizout(2)); |
| 249 |
|
|
tmpi=tmpi(:); tmpj=prof_interp.vpi(:); tmpw=prof_interp.w(:); |
| 250 |
|
|
kk=find(~isnan(tmpj)); |
| 251 |
|
|
tmpi=tmpi(kk); tmpj=tmpj(kk); tmpw=tmpw(kk); |
| 252 |
|
|
prof_interp.SPM=sparse(tmpi,tmpj,tmpw,sizout(1),sizin,prod(sizout)); |
| 253 |
|
|
prof_interp=rmfield(prof_interp,'vpi'); |
| 254 |
|
|
|
| 255 |
gforget |
1.5 |
if doVerbose; |
| 256 |
|
|
n1=sum(~isnan(prof_interp.face)); |
| 257 |
|
|
n2=sum(isnan(prof_interp.face)); |
| 258 |
|
|
fprintf(['interpolated points: ' num2str(n1) '\n']); |
| 259 |
|
|
fprintf(['un-treated points: ' num2str(n2) '\n']); |
| 260 |
|
|
end; |
| 261 |
gforget |
1.1 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph