Issuing rollback() due to DESTROY without explicit disconnect() of DBD::mysql::db handle dbname=MITgcm at /usr/local/share/lxr/lib/LXR/Common.pm line 1224.
Last-Modified: Fri, 28 Nov 2024 06:11:57 GMT
Content-Type: text/html; charset=utf-8
MITgcm/MITgcm/utils/matlab/cs_grid/latloncap/compact2map.m
c391107ae5 Jean*0001 function b=compact2map(a,flag,rgbDim);
0002 % function b=compact2map(a,[flag,rgbDim]);
0003 % convert between mdsio compact format and a sensible 2d map for llc fields;
0004 % unless specified explicitly (optional argument rgbDim), try to guess
0005 % (assuming regular Lat-Lon cap) what the 3 dimensions of faces are.
176918a0f7 Mart*0006 % a: input field (2d, 3d, or 4d)
0007 % b: output field
0008 % flag == 0 (default): convert from compact format to 2d map
0009 % flag == 1 : convert from 2d map to compact format
00100011 % author: Martin Losch (Martin.Losch@awi.de)
0012 if nargin < 2
0013 flag = 0;
0014 end
0015 if flag
0016 % convert from "nice" to compact format
0017 [nx,ny,nz,nt]=size(a);
0018 if nx/4~=round(nx/4);
0019 error('not a llc field');
0020 end
00210022 n4=nx/4;
0023 m=(ny-n4)/n4;
0024 m4=m*n4;
0025 for kt = 1:nt
0026 for kz = 1:nz
0027 for k=1:4
0028 sides{k} = a([1:n4]+(k-1)*n4,1:m4,kz,kt);
0029 end
0030 % cap = a([1:n4]+1*n4,m4+1:end,kz,kt);
0031 cap = rot90(a([1:n4],m4+1:end,kz,kt),-1);
0032 % reformat
0033 btmp = [sides{1} ...
0034 sides{2} ...
0035 cap ...
0036 reshape(rot90(sides{3},1),[n4 m4]) ...
0037 reshape(rot90(sides{4},1),[n4 m4]) ...
0038 cap*0];
0039 b(:,:,kz,kt)=btmp;
0040 end
0041 end
0042 else
0043 % convert from compact to "nice" format
c391107ae5 Jean*0044 dims=size(a); ndim=length(dims);
0045 if ndim == 2, dims=[dims 1]; end; n3d=prod(dims(3:end));
0046 nPts=prod(dims(1:2)); a=reshape(a,[nPts n3d]);
0047 if nargin < 3,
0048 nx=dims(1); m=dims(2);
0049 if m/nx~=round(m/nx);
0050 error('not a compact format llc field');
0051 end
0052 ntiles = m/nx;
0053 if mod(ntiles,2); % odd number of tiles
0054 tilesPerFace = (ntiles-1)/4;
0055 else
0056 tilesPerFace = (ntiles-2)/4;
0057 end
0058 ng=tilesPerFace*nx;
0059 nr=nx; nb=nx; Nfaces=5;
0060 elseif length(rgbDim) ~= 3,
0061 error('3 elements array expected as 3rd argument (rgbDim)');
176918a0f7 Mart*0062 else
c391107ae5 Jean*0063 nr=rgbDim(1); ng=rgbDim(2); nb=rgbDim(3);
176918a0f7 Mart*0064 end
c391107ae5 Jean*0065 %-- set dimensions of each of the 6 faces:
0066 nf=ones(6,2);
0067 nf(1,:)=[nr ng]; nf(2,:)=[nb ng]; nf(3,:)=[nb nr];
0068 nf(4,:)=[ng nr]; nf(5,:)=[ng nb]; nf(6,:)=[nr nb];
0069 fdim=prod(nf,2); fd2= cumsum(fdim); fd1=fd2-fdim+1;
0070 if nargin == 3,
0071 [Nfaces]=find(fd2==nPts);
0072 if length(Nfaces) ~= 1,
0073 error('array size does not match 5 or 6 faces compact array dim');
176918a0f7 Mart*0074 end
0075 end
c391107ae5 Jean*0076 %-- set output 2-D map dimensions:
0077 n2=nr*2+nb*2;
0078 j1=0; if Nfaces==6, j1=max(nr,nb); end
0079 j2=j1+ng; m2=j2+max(nr,nb);
0080 b=zeros(n2,m2,n3d);
00810082 for kz = 1:n3d
0083 i1=1; i2=0;
0084 for n=1:Nfaces,
0085 i1=i2+1;
00860087 %-- extract face "n" data from compact array "a", level kz:
0088 vf1=reshape(a(fd1(n):fd2(n),kz),[nf(n,:)]);
0089 %-- put face data into a "nice" 2-D map:
0090 if n < 3,
0091 %- face 1 & 2 : strait copy
0092 i2=i2+nf(n,1);
0093 b(i1:i2,j1+1:j2,kz)=vf1;
0094 elseif n == 3,
0095 %- face 3 is copied 4 times (in continuation of each equatorial faces)
0096 ii1=1; ii2=nr;
0097 for k=1:4,
0098 if rem(k,2) == 0, mj=nr; else mj=nb ; end
0099 if k == 2,
0100 b(ii1:ii2,j2+1:j2+mj,kz)=vf1;
0101 else
0102 b(ii1:ii2,j2+1:j2+mj,kz)=rot90(vf1,2-k);
0103 end
0104 ii1=ii2+1;
0105 ii2=ii2+mj;
0106 end
0107 elseif n == 6,
0108 %- face 6 is copied 4 times (in continuation of each equatorial faces)
0109 ii1=1; ii2=nr;
0110 for k=1:4,
0111 if rem(k,2) == 0, mj=nr; else mj=nb ; end
0112 if k == 1,
0113 b(ii1:ii2,j1-mj+1:j1,kz)=vf1;
0114 else
0115 b(ii1:ii2,j1-mj+1:j1,kz)=rot90(vf1,k-1);
0116 end
0117 ii1=ii2+1;
0118 ii2=ii2+mj;
0119 end
0120 else
0121 %- face 4 & 5 : rotate & copy
0122 i2=i2+nf(n,2);
0123 b(i1:i2,j1+1:j2,kz)=rot90(vf1,-1);
0124 end
01250126 end
0127 end
0128 if ndim > 3, b=reshape(b,[n2 m2 dims(3:end)]);
0129 elseif ndim < 3, b=squeeze(b); end
0130176918a0f7 Mart*0131 end
01320133 return
