Warning, /utils/matlab/cs_grid/cubeZ2latlon.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
896b20e91e Jean*0001 function [z] = cubeZ2latlon(x,y,c,xi,yi)
0002 % z=cubeZ2latlon(x,y,c,xi,yi);
0003 %
0004 % Re-grids model output on expanded spherical cube to lat-lon grid.
0005 % x,y are 1-D arrays (reshaped Horiz. grid) of the cell-corner coordinates
0006 % c is a 1-D or 2-D scalar field
0007 % xi,yi are vectors of the new regular lat-lon grid to interpolate to.
0008 % z is the interpolated data with dimensions of size(xi) by size(yi).
0009 %
0010 % e.g.
0011 % >> xg=rdmds('XG'); nPts=prod(size(xg)); x=rehsape(xg,nPts,1);
0012 % >> yg=rdmds('YG'); y=rehsape(yg,nPts,1);
0013 % >> PsiH=rdmds('psiH.0000513360');
0014 % >> xi=-179:2:180;yi=-89:2:90;
0015 % >> psiI=cubeZ2latlon(x,y,psiH,xi,yi);
0016 % can also add the 2 missing corner :
0017 % >> nc=sqrt(nPts/6);
0018 % >> x(end+1)=mean(xg([1:2*nc:6*nc],nc)); y(end+1)=mean(yg([1:2*nc:6*nc],nc));
0019 % >> x(end+1)=mean(xg([2*nc:2*nc:6*nc],1)); y(end+1)=mean(yg([2*nc:2*nc:6*nc],1));
0020 %
0021 % Written by jmc@ocean.mit.edu, 2005.
0022 NN=size(c);
0023 [nPt2 nz]=size(c);
0024 nc=fix(sqrt(nPt2/6)); nPts=6*nc*nc;
0025
0026 for k=1:nz;
0027 X=x(1:nPts);X=reshape(X,6*nc,nc);
0028 Y=y(1:nPts);Y=reshape(Y,6*nc,nc);
0029 C=c(1:nPts,k);C=reshape(C,6*nc,nc);
0030
0031 i=3*nc+(1:nc);j=floor(nc/2)+1;
0032 X(end+1,:)=X(i,j)'-360; Y(end+1,:)=Y(i,j)'; C(end+1,:)=C(i,j)';
0033 %i=3*nc+(1:nc);j=floor(nc/2)+1;
0034 %X(end+1,:)=X(i,j)'+360; Y(end+1,:)=Y(i,j)'; C(end+1,:)=C(i,j)';
0035 i=5*nc+1+floor(nc/2);j=1:floor(nc/2);
0036 X(end+1,j)=X(i,j)-360;
0037 Y(end+1,j)=Y(i,j);
0038 C(end+1,j)=C(i,j);
0039 %--
0040 j=1+floor(nc/2);i=2*nc+j; if Y(i,j)==90, X(i,j)=180; end
0041 i=2*nc+(floor(nc/2)+1:nc);j=1+floor(nc/2);
0042 X(end,i-2*nc)=X(i,j)'-360;
0043 Y(end,i-2*nc)=Y(i,j)';
0044 C(end,i-2*nc)=C(i,j)';
0045 %--
0046 %i=5*nc+floor(nc/2)+1;j=1:floor(nc/2);
0047 %X(end,j+nc/2)=X(i,j)-360;
0048 %Y(end,j+nc/2)=Y(i,j);
0049 %C(end,j+nc/2)=C(i,j);
0050 %i=2*nc+(nc/2+1:nc);j=floor(nc/2);
0051 %X(end+1,1:nc/2)=X(i,j)'-360;
0052 %Y(end+1,1:nc/2)=Y(i,j)';
0053 %C(end+1,1:nc/2)=C(i,j)';
0054 %i=2*nc+(nc/2+1:nc);j=floor(nc/2)+1;
0055 %X(end,nc/2+1:nc)=X(i,j)'+360;
0056 %Y(end,nc/2+1:nc)=Y(i,j)';
0057 %C(end,nc/2+1:nc)=C(i,j)';
0058
0059 j=1+floor(nc/2);i=5*nc+j; ij=i+(j-1)*nc*6;
0060 if Y(i,j)==-90,
0061 % fprintf('South pole: %i %i %f %f\n',i,j,X(i,j),Y(i,j));
0062 X(i,j)=180;
0063 end
0064
0065 X=reshape(X,prod(size(X)),1);
0066 Y=reshape(Y,prod(size(Y)),1);
0067 C=reshape(C,prod(size(C)),1);
0068
0069 [I]=find(Y==-90);
0070 if length(I)==1,
0071 % fprintf('South pole: %i %f %f\n',I,X(I),Y(I));
0072 X(end+1)=X(I)-360;
0073 Y(end+1)=Y(I);
0074 C(end+1)=C(I);
0075 end
0076
0077 if nPt2 > nPts,
0078 X(end+1)=x(nPts+1);
0079 Y(end+1)=y(nPts+1);
0080 C(end+1)=c(nPts+1,k);
0081 end
0082 if nPt2 == nPts+2,
0083 X(end+1)=x(nPt2);
0084 Y(end+1)=y(nPt2);
0085 C(end+1)=c(nPt2,k);
0086 end
0087
0088 z(:,:,k)=griddata(Y,X,C,yi,xi');
0089 end % k
0090
0091 if length(NN)>1
0092 z=reshape(z,[size(z,1) size(z,2) NN(2:end)]);
0093 end
