| 1 |
function [z] = cubeZ2latlon(x,y,c,xi,yi) |
| 2 |
% z=cubeZ2latlon(x,y,c,xi,yi); |
| 3 |
% |
| 4 |
% Re-grids model output on expanded spherical cube to lat-lon grid. |
| 5 |
% x,y are 1-D arrays (reshaped Horiz. grid) of the cell-corner coordinates |
| 6 |
% c is a 1-D or 2-D scalar field |
| 7 |
% xi,yi are vectors of the new regular lat-lon grid to interpolate to. |
| 8 |
% z is the interpolated data with dimensions of size(xi) by size(yi). |
| 9 |
% |
| 10 |
% e.g. |
| 11 |
% >> xg=rdmds('XG'); nPts=prod(size(xg)); x=rehsape(xg,nPts,1); |
| 12 |
% >> yg=rdmds('YG'); y=rehsape(yg,nPts,1); |
| 13 |
% >> PsiH=rdmds('psiH.0000513360'); |
| 14 |
% >> xi=-179:2:180;yi=-89:2:90; |
| 15 |
% >> psiI=cubeZ2latlon(x,y,psiH,xi,yi); |
| 16 |
% can also add the 2 missing corner : |
| 17 |
% >> nc=sqrt(nPts/6); |
| 18 |
% >> x(end+1)=mean(xg([1:2*nc:6*nc],nc)); y(end+1)=mean(yg([1:2*nc:6*nc],nc)); |
| 19 |
% >> x(end+1)=mean(xg([2*nc:2*nc:6*nc],1)); y(end+1)=mean(yg([2*nc:2*nc:6*nc],1)); |
| 20 |
% |
| 21 |
% $Header: /u/gcmpack/MITgcm/utils/cs_grid/cubeZ2latlon.m,v 1.2 2005/08/05 23:08:09 jmc Exp $ |
| 22 |
|
| 23 |
NN=size(c); |
| 24 |
[nPt2 nz]=size(c); |
| 25 |
nc=fix(sqrt(nPt2/6)); nPts=6*nc*nc; |
| 26 |
|
| 27 |
for k=1:nz; |
| 28 |
X=x(1:nPts);X=reshape(X,6*nc,nc); |
| 29 |
Y=y(1:nPts);Y=reshape(Y,6*nc,nc); |
| 30 |
C=c(1:nPts,k);C=reshape(C,6*nc,nc); |
| 31 |
|
| 32 |
i=3*nc+(1:nc);j=floor(nc/2)+1; |
| 33 |
X(end+1,:)=X(i,j)'-360; Y(end+1,:)=Y(i,j)'; C(end+1,:)=C(i,j)'; |
| 34 |
%i=3*nc+(1:nc);j=floor(nc/2)+1; |
| 35 |
%X(end+1,:)=X(i,j)'+360; Y(end+1,:)=Y(i,j)'; C(end+1,:)=C(i,j)'; |
| 36 |
i=5*nc+1+floor(nc/2);j=1:floor(nc/2); |
| 37 |
X(end+1,j)=X(i,j)-360; |
| 38 |
Y(end+1,j)=Y(i,j); |
| 39 |
C(end+1,j)=C(i,j); |
| 40 |
%-- |
| 41 |
j=1+floor(nc/2);i=2*nc+j; if Y(i,j)==90, X(i,j)=180; end |
| 42 |
i=2*nc+(floor(nc/2)+1:nc);j=1+floor(nc/2); |
| 43 |
X(end,i-2*nc)=X(i,j)'-360; |
| 44 |
Y(end,i-2*nc)=Y(i,j)'; |
| 45 |
C(end,i-2*nc)=C(i,j)'; |
| 46 |
%-- |
| 47 |
%i=5*nc+floor(nc/2)+1;j=1:floor(nc/2); |
| 48 |
%X(end,j+nc/2)=X(i,j)-360; |
| 49 |
%Y(end,j+nc/2)=Y(i,j); |
| 50 |
%C(end,j+nc/2)=C(i,j); |
| 51 |
%i=2*nc+(nc/2+1:nc);j=floor(nc/2); |
| 52 |
%X(end+1,1:nc/2)=X(i,j)'-360; |
| 53 |
%Y(end+1,1:nc/2)=Y(i,j)'; |
| 54 |
%C(end+1,1:nc/2)=C(i,j)'; |
| 55 |
%i=2*nc+(nc/2+1:nc);j=floor(nc/2)+1; |
| 56 |
%X(end,nc/2+1:nc)=X(i,j)'+360; |
| 57 |
%Y(end,nc/2+1:nc)=Y(i,j)'; |
| 58 |
%C(end,nc/2+1:nc)=C(i,j)'; |
| 59 |
|
| 60 |
j=1+floor(nc/2);i=5*nc+j; ij=i+(j-1)*nc*6; |
| 61 |
if Y(i,j)==-90, |
| 62 |
% fprintf('South pole: %i %i %f %f\n',i,j,X(i,j),Y(i,j)); |
| 63 |
X(i,j)=180; |
| 64 |
end |
| 65 |
|
| 66 |
X=reshape(X,prod(size(X)),1); |
| 67 |
Y=reshape(Y,prod(size(Y)),1); |
| 68 |
C=reshape(C,prod(size(C)),1); |
| 69 |
|
| 70 |
[I]=find(Y==-90); |
| 71 |
if length(I)==1, |
| 72 |
% fprintf('South pole: %i %f %f\n',I,X(I),Y(I)); |
| 73 |
X(end+1)=X(I)-360; |
| 74 |
Y(end+1)=Y(I); |
| 75 |
C(end+1)=C(I); |
| 76 |
end |
| 77 |
|
| 78 |
if nPt2 > nPts, |
| 79 |
X(end+1)=x(nPts+1); |
| 80 |
Y(end+1)=y(nPts+1); |
| 81 |
C(end+1)=c(nPts+1,k); |
| 82 |
end |
| 83 |
if nPt2 == nPts+2, |
| 84 |
X(end+1)=x(nPt2); |
| 85 |
Y(end+1)=y(nPt2); |
| 86 |
C(end+1)=c(nPt2,k); |
| 87 |
end |
| 88 |
|
| 89 |
z(:,:,k)=griddata(Y,X,C,yi,xi'); |
| 90 |
end % k |
| 91 |
|
| 92 |
if length(NN)>1 |
| 93 |
z=reshape(z,[size(z,1) size(z,2) NN(2:end)]); |
| 94 |
end |
Parent Directory
| 
Revision Log
| 
Revision Graph