| 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 |
|
|
doVerbose=1; if myenv.verbose>=1; doVerbose=1; end; |
| 21 |
|
|
|
| 22 |
|
|
%set-up tile information (domain decomposition to ni,nj blocs) |
| 23 |
|
|
if nargin<=2; |
| 24 |
|
|
ni=30; nj=30; |
| 25 |
|
|
else; |
| 26 |
|
|
ni=varargin{1}; |
| 27 |
|
|
nj=varargin{2}; |
| 28 |
|
|
end; |
| 29 |
|
|
|
| 30 |
|
|
map_tile=gcmfaces_loc_tile(ni,nj); |
| 31 |
|
|
loc_tile=gcmfaces_loc_tile(ni,nj,prof_lon,prof_lat); |
| 32 |
|
|
prof_tile=loc_tile.tileNo; |
| 33 |
|
|
list_tile=unique(prof_tile); |
| 34 |
|
|
|
| 35 |
|
|
%initialize output: |
| 36 |
|
|
prof_interp.face=NaN*prof_lon; |
| 37 |
|
|
prof_interp.tile=NaN*prof_lon; |
| 38 |
|
|
prof_interp.i=NaN*repmat(prof_lon,[1 4]); |
| 39 |
|
|
prof_interp.j=NaN*repmat(prof_lon,[1 4]); |
| 40 |
|
|
prof_interp.w=NaN*repmat(prof_lon,[1 4]); |
| 41 |
gforget |
1.2 |
% |
| 42 |
|
|
tile_corners.XC11=NaN*prof_lon; |
| 43 |
|
|
tile_corners.YC11=NaN*prof_lon; |
| 44 |
|
|
tile_corners.XCNINJ=NaN*prof_lon; |
| 45 |
|
|
tile_corners.YCNINJ=NaN*prof_lon; |
| 46 |
gforget |
1.1 |
|
| 47 |
|
|
%loop over tiles |
| 48 |
|
|
for ii=1:length(list_tile); |
| 49 |
|
|
|
| 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 |
|
|
end; |
| 113 |
|
|
if length(unique(JJ))<length(ii_prof)&doVerbose; |
| 114 |
|
|
n0=num2str(length(ii_prof)-length(unique(JJ))); |
| 115 |
|
|
n1=num2str(length(ii_prof)); |
| 116 |
|
|
warning(['no polygon for ' n0 ' / ' n1]); |
| 117 |
|
|
%the following will then remove the corresponding profiles form ii_prof |
| 118 |
|
|
end; |
| 119 |
|
|
[C,IA,IC] = unique(JJ); |
| 120 |
|
|
% |
| 121 |
|
|
ii_prof0=ii_prof; |
| 122 |
|
|
ii_prof=ii_prof(C);%treated profiles |
| 123 |
|
|
jj_prof=setdiff(ii_prof0,ii_prof);%un-treated profiles |
| 124 |
|
|
% |
| 125 |
|
|
ii_quad=II(IA); |
| 126 |
|
|
|
| 127 |
|
|
if doDisplay; |
| 128 |
|
|
figureL; |
| 129 |
|
|
xx_tile=xxx{iiFace}(iiMin:iiMax+2,jjMin:jjMax+2); |
| 130 |
|
|
yy_tile=yyy{iiFace}(iiMin:iiMax+2,jjMin:jjMax+2); |
| 131 |
|
|
pcolor(xx_tile,yy_tile,sqrt(xx_tile.^2+yy_tile.^2)); hold on; |
| 132 |
|
|
cc=caxis; cc(2)=cc(2)*2; caxis(cc); |
| 133 |
|
|
plot(prof_x(ii_prof),prof_y(ii_prof),'r.','MarkerSize',20); |
| 134 |
|
|
plot(prof_x(jj_prof),prof_y(jj_prof),'k.','MarkerSize',60); |
| 135 |
|
|
axis([-0.6 0.6 -0.6 0.6]/6); |
| 136 |
|
|
end; |
| 137 |
|
|
|
| 138 |
|
|
if ~isempty(ii_prof); |
| 139 |
|
|
%5) determine bi-linear interpolation weights: |
| 140 |
|
|
px=x_quad(ii_quad,:); |
| 141 |
|
|
py=y_quad(ii_quad,:); |
| 142 |
|
|
ox=prof_x(ii_prof); |
| 143 |
|
|
oy=prof_y(ii_prof); |
| 144 |
|
|
[ow]=gcmfaces_quadmap(px,py,ox,oy); |
| 145 |
|
|
|
| 146 |
|
|
%to double check interpolation |
| 147 |
|
|
% pw=squeeze(ow); |
| 148 |
|
|
% oxInterp=sum(pw.*px,2); |
| 149 |
|
|
% oyInterp=sum(pw.*py,2); |
| 150 |
gforget |
1.3 |
|
| 151 |
|
|
%round up coefficient to 4th digit (also to avoid slight negatives) |
| 152 |
|
|
test1=~isempty(find(ow(:)<-1e-5)); |
| 153 |
|
|
if test1; error('interp weight < 0 -- something went wrong'); end; |
| 154 |
|
|
test1=~isempty(find(ow(:)>1+1e-5)); |
| 155 |
|
|
if test1; error('interp weight >1 -- something went wrong'); end; |
| 156 |
|
|
% |
| 157 |
|
|
ow=1e-4*round(ow*1e4); |
| 158 |
|
|
sumw=repmat(sum(ow,3),[1 1 4]); |
| 159 |
|
|
ow=ow./sumw; |
| 160 |
gforget |
1.1 |
|
| 161 |
|
|
%6) output interpolation specs: |
| 162 |
|
|
prof_interp.face(ii_prof,1)=iiFace*(1+0*ii_quad); |
| 163 |
|
|
prof_interp.tile(ii_prof,1)=list_tile(ii)*(1+0*ii_quad); |
| 164 |
|
|
prof_interp.i(ii_prof,:)=i_quad(ii_quad,:); |
| 165 |
|
|
prof_interp.j(ii_prof,:)=j_quad(ii_quad,:); |
| 166 |
|
|
prof_interp.w(ii_prof,:)=squeeze(ow); |
| 167 |
gforget |
1.2 |
% |
| 168 |
|
|
tile_corners.XC11(ii_prof)=XC11; |
| 169 |
|
|
tile_corners.YC11(ii_prof)=YC11; |
| 170 |
|
|
tile_corners.XCNINJ(ii_prof)=XCNINJ; |
| 171 |
|
|
tile_corners.YCNINJ(ii_prof)=YCNINJ; |
| 172 |
gforget |
1.1 |
end; |
| 173 |
|
|
|
| 174 |
|
|
end;%for ii=1:length(list_tile); |
| 175 |
|
|
|
| 176 |
|
|
% if doVerbose; |
| 177 |
|
|
n1=sum(~isnan(prof_interp.face)); |
| 178 |
|
|
n2=sum(isnan(prof_interp.face)); |
| 179 |
|
|
fprintf(['interpolated points: ' num2str(n1) '\n']); |
| 180 |
|
|
fprintf(['un-treated points: ' num2str(n2) '\n']); |
| 181 |
|
|
% end; |
| 182 |
|
|
|
Parent Directory
| 
Revision Log
| 
Revision Graph