Warning, /utils/matlab/cs_grid/bk_line/grph_CSz.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit 6b4f212f on 2024-06-05 16:55:17 UTC
e8997c8981 Jean*0001 function grph_CSz(var,xcs,ycs,xcg,ycg,c1,c2,shift,cbV,AxBx,kEnv)
0002 % grph_CSz(var,xcs,ycs,xcg,ycg,c1,c2,shift[,cbV,AxBx,kEnv]) :
0003 % produce a flat plot of the cube_sphere "var" (@ C-grid corner position)
0004 % keeping the initial grid (no interpolation, use "surf").
323aa12fe2 Jean*0005 % xcs,ycs,xcg,ycg = center + corner grid point coordinates
0006 % c1 < c2 = min & max for the color graph
0007 % c1 > c2 = scale with min,max of the field, + c1/100 and + c2/100
040160fb38 Andr*0008 % shift=-1 : No coast-line
0009 % shift= 1 : No shift but draw coast-line calling draw_coast
323aa12fe2 Jean*0010 % else : draw coast-line (using m_proj) shifted by "shift" degree.
e8997c8981 Jean*0011 % cbV = 0,1 : horizontal,vertical colorbar; >= 2 : no colorbar;
0012 % kEnv = 0 : standard ; =odd : do not draw the mesh ; >1 : no min,Max written.
323aa12fe2 Jean*0013 % AxBx = do axis(AxBx) to zoom in Box "AxBx" ; only if shift=-1 ;
0014 %-----------------------
6b4f212f67 Jean*0015
4ec37fd829 Jean*0016 % Written by jmc@mit.edu, 2005.
0017
6b4f212f67 Jean*0018 %- for debugging, switch to > 0:
0019 dBug=0;
1c4dce2004 Jean*0020 %- small number (relative to lon,lat in degree)
0021 epsil=1.e-6;
6b4f212f67 Jean*0022 %- big jump in longitude from 2 neighbor points
0023 xJump=240; xJmp2=xJump*0.5;
1c4dce2004 Jean*0024 %- mid-longitude of the grid (truncated @ epsil level):
0025 xMid=mean(xcs(:)); xMid=epsil*round(xMid/epsil);
6b4f212f67 Jean*0026 if dBug > 0, fprintf(' mid longitude of the grid: xMid= %22.16e\n',xMid); end
1c4dce2004 Jean*0027
323aa12fe2 Jean*0028 if nargin < 9, cbV=0 ; end
1c4dce2004 Jean*0029 if nargin < 10, AxBx=[xMid-180 xMid+180 -90 90] ; end
e8997c8981 Jean*0030 if nargin < 11, kEnv=0 ; end
323aa12fe2 Jean*0031
0032 %------------------------------
4ec37fd829 Jean*0033 n1h=size(xcs,1); n2h=size(xcs,2);
0034 if n1h == 6*n2h, nc=n2h;
0035 elseif n1h*6 == n2h, nc=n1h;
0036 else
0037 error([' input xcs size: ',int2str(n1h),' x ',int2str(n2h),' does not fit regular cube !']);
0038 end
0039 ncp=nc+1 ; nPg=nc*nc*6;
0040 nPp2=size(var,1);
0041 %- check input "var" size: assumes it includes the 2 missing corner (nPg+2 lenth) and,
0042 % as a long-vector, always in "compact-format" shape (i.e., 1 face after the other)
0043 if nPp2 ~= nPg+2,
0044 fprintf('Bad dim (input fields): n1h,n2h,nPg,nPp2= %i %i %i %i\n',n1h,n2h,nPg,nPp2);
e8997c8981 Jean*0045 error(' wrong dimensions ');
323aa12fe2 Jean*0046 end
0047 mnV=min(var);
0048 MxV=max(var);
0049 % mx=min(min(var));
0050 % mn=max(max(var));
040160fb38 Andr*0051 %if shift == 1,
323aa12fe2 Jean*0052 % for j=1:nc, for i=1:6*nc,
0053 % if var(i,j) ~= NaN & xcs(i,j) > AxBx(1) & xcs(i,j) < AxBx(2) ...
0054 % & ycs(i,j) > AxBx(3) & ycs(i,j) < AxBx(4) ,
0055 % mn=min(var(i,j),mn); mx=max(var(i,j),mx) ; end
0056 % end ; end ;
0057 %else
0058 % for j=1:nc, for i=1:6*nc,
0059 % if var(i,j) ~= NaN ; mn=min(var(i,j),mn); mx=max(var(i,j),mx) ; end
0060 % end ; end ;
0061 %end
0062 %------------
0063 if c1 >= c2
0064 mb=(MxV-mnV)*0.01;
0065 c1=mnV+mb*c1;
0066 c2=MxV+mb*c2;
1c4dce2004 Jean*0067 if c1*c2 < 0
323aa12fe2 Jean*0068 c2=max(abs([c1 c2]));
0069 c1=-c2;
0070 end
0071 fprintf('min,max %8.3e %8.3e Cmin,max %8.3e %8.3e \n',mnV,MxV,c1,c2)
0072 end
0073 %------------------------------
0074 %figure(1);
040160fb38 Andr*0075 if shift ~= 1 & shift ~= -1,
323aa12fe2 Jean*0076 %pphh=path; ppaa='/u/u2/jmc/MATLAB';
0077 pphh=path; ppaa='/home/jmc/matlab/MATLAB';
0078 llaa=length(ppaa); if pphh(1:llaa) == ppaa, else path(ppaa,path);end
0079 set_axis
0080 fac=rad ;
0081 else
0082 fac=1. ;
0083 end
0084 %---
0085 nbsf = 0 ; ic = 0 ; jc = 0 ;
4ec37fd829 Jean*0086 %- long-vector input "var" always in "compact-format" style (i.e., 1 face after the other)
0087 vv0=reshape(var(1:nPg,1),[nc nc 6]);
0088 [xx2]=split_C_cub(xcs,1);
0089 [yy2]=split_C_cub(ycs,1);
323aa12fe2 Jean*0090 %---
0091 for n=1:6,
1c4dce2004 Jean*0092 %if n > 2 & n < 4,
323aa12fe2 Jean*0093 if n < 7,
0094 %--------------------------------------------------------
0095 vv1=zeros(ncp,ncp) ; xx1=vv1 ; yy1=vv1 ;
1c4dce2004 Jean*0096 vv1(1:nc,1:nc)=vv0(:,:,n);
6b4f212f67 Jean*0097 xx1=xx2(:,:,n); yy1=yy2(:,:,n); xxc=xx1(2:ncp,2:ncp);
862f2fe640 Jean*0098 % if xx1(ncp,1) < xMid-300. ; xx1(ncp,1)=xx1(ncp,1)+360. ; end
0099 % if xx1(1,ncp) < xMid-300. ; xx1(1,ncp)=xx1(1,ncp)+360. ; end
323aa12fe2 Jean*0100 %------------
6b4f212f67 Jean*0101 %- use jump in grid-cell center longitude to decide where to cut 1 face:
0102 cutFace=0;
0103 dxi=xx1(2:ncp,:)-xx1(1:nc,:); dxImx=max(abs(dxi(:)));
0104 dxj=xx1(:,2:ncp)-xx1(:,1:nc); dxJmx=max(abs(dxj(:)));
0105 if dxImx > xJump & dxJmx > xJump, cutFace=3;
0106 elseif dxImx > xJump, cutFace=1;
0107 elseif dxJmx > xJump, cutFace=2; end
0108 if dBug > 0,
0109 fprintf(' face # %i , Max dxI,dxJ = %8.3f , %8.3f',n,dxImx,dxJmx);
0110 if cutFace > 0, fprintf(' ; cutFace= %i',cutFace); end
0111 fprintf('\n');
0112 end
0113 if cutFace == 3,
0114 fprintf(' Jump in both i & j not implemented ==> skip face # %i\n',n);
0115 %------------
0116 elseif cutFace == 2,
0117 [I,J]=find( abs(dxj) > xJump );
0118 if dBug > 1, for l=1:length(I), fprintf(' i,j= %2i, %2i \n',I(l),J(l)); end; end
0119 if min(J) == max(J),
0120 jc =J(1); jp=jc+1;
0121 if dBug > 0, fprintf('--> cut Face @ jc,jp = %3i,%3i\n',jc,jp); end
0122 % cut the face in 2 parts (1: j in [1 jp] ; 2: j in [jc ncp]) and plot separately
0123 % note: duplicate points at the edges to get half grid-cell plotted on both side
0124 for lp=1:2,
0125 if lp == 1,
0126 j1=1; j2=jp; j3=jc-1;
0127 else
0128 j1=jc; j2=ncp; j3=nc;
0129 end
0130 xx1=xx2(:,:,n);
0131 tmp=xxc(:,j1:j3); xLoc=mean(tmp(:));
0132 if dBug > 0, fprintf(' p.%i, %2i -> %2i : %8.3f %8.3f %8.3f\n', ...
0133 lp,j1,j3,min(tmp(:)),xLoc,max(tmp(:))); end
0134 if xLoc > xMid,
0135 %-- case where Xc jump from > 180 to < -180 when j goes from jc to jc+1 :
0136 [I]=find(xx1(:,j1) < xMid-xJmp2); xx1(I,j1)=xx1(I,j1)+360.;
0137 [I]=find(xx1(:,j2) < xMid-xJmp2); xx1(I,j2)=xx1(I,j2)+360.;
0138 else
0139 %-- case where Xc jump from < -180 to > 180 when j goes from jc to jc+1 :
0140 [I]=find(xx1(:,j1) > xMid+xJmp2); xx1(I,j1)=xx1(I,j1)-360.;
0141 [I]=find(xx1(:,j2) > xMid+xJmp2); xx1(I,j2)=xx1(I,j2)-360.;
0142 end
0143 [nbsf,S(nbsf)]=part_surf(nbsf,fac,xx1,yy1,vv1,1,ncp,j1,j2,c1,c2) ;
0144 end
862f2fe640 Jean*0145 % Note: later on, will plot separately N & S poles with the 2 missing corners.
6b4f212f67 Jean*0146 else
0147 fprintf(' Irregular cut not implemented ==> skip face # %i\n',n);
0148 end
0149 %------------
0150 elseif cutFace == 1,
0151 [I,J]=find( abs(dxi) > xJump );
0152 if dBug > 1, for l=1:length(I), fprintf(' i,j= %2i, %2i \n',I(l),J(l)); end; end
0153 if min(I) == max(I),
0154 ic =I(1); ip=ic+1;
0155 if dBug > 0, fprintf('--> cut Face @ ic,ip = %3i,%3i\n',ic,ip); end
0156 % cut the face in 2 parts (1: i in [1 ip] ; 2: i in [ic ncp]) and plot separately
0157 % note: duplicate points at the edges to get half grid-cell plotted on both side
0158 for lp=1:2,
0159 if lp == 1,
0160 i1=1; i2=ip; i3=ic-1;
0161 else
0162 i1=ic; i2=ncp; i3=nc;
0163 end
0164 xx1=xx2(:,:,n);
0165 tmp=xxc(i1:i3,:); xLoc=mean(tmp(:));
0166 if dBug > 0, fprintf(' p.%i, %2i -> %2i : %8.3f %8.3f %8.3f\n', ...
0167 lp,i1,i3,min(tmp(:)),xLoc,max(tmp(:))); end
0168 if xLoc > xMid,
0169 %-- case where X jump from > 180 to < -180 when i goes from ic to ic+1 :
0170 [J]=find(xx1(i1,:) < xMid-xJmp2); xx1(i1,J)=xx1(i1,J)+360.;
0171 [J]=find(xx1(i2,:) < xMid-xJmp2); xx1(i2,J)=xx1(i2,J)+360.;
0172 else
0173 %-- case where X jump from < -180 to > 180 when i goes from ic to ic+1 :
0174 [J]=find(xx1(i1,:) > xMid+xJmp2); xx1(i1,J)=xx1(i1,J)-360.;
0175 [J]=find(xx1(i2,:) > xMid+xJmp2); xx1(i2,J)=xx1(i2,J)-360.;
0176 end
0177 [nbsf,S(nbsf)]=part_surf(nbsf,fac,xx1,yy1,vv1,i1,i2,1,ncp,c1,c2) ;
0178 end
0179 %----
0180 else
0181 fprintf(' Irregular cut not implemented ==> skip face # %i\n',n);
0182 end
0183 %------------
323aa12fe2 Jean*0184 else
862f2fe640 Jean*0185 %-- plot the face in 1 piece :
6b4f212f67 Jean*0186 xLoc=mean(xxc(:));
0187 xx1=rem(xx1-xLoc+180*3,360)+xLoc-180;
323aa12fe2 Jean*0188 [nbsf,S(nbsf)]=part_surf(nbsf,fac,xx1,yy1,vv1,1,ncp,1,ncp,c1,c2) ;
0189 end
0190 %--------------------------------------
1c4dce2004 Jean*0191 end ; end
323aa12fe2 Jean*0192 %--------------
0193 %- add isolated point:
0194 vvI=zeros(2,2); xxI=vvI; yyI=vvI;
0195
4ec37fd829 Jean*0196 %- 1rst missing corner (N.W corner of 1rst face): nPg+1
0197 vvI(1,1)=var(nPg+1);
1c4dce2004 Jean*0198 for l=0:2,
323aa12fe2 Jean*0199 xxI(1+rem(l,2),1+fix(l/2))=xx2(2,ncp,1+2*l);
0200 yyI(1+rem(l,2),1+fix(l/2))=yy2(2,ncp,1+2*l);
0201 end
0202 xxI(2,2)=xxI(1,2); yyI(2,2)=yyI(1,2);
0203 [nbsf,S(nbsf)]=part_surf(nbsf,fac,xxI,yyI,vvI,1,2,1,2,c1,c2) ;
0204
4ec37fd829 Jean*0205 %- 2nd missing corner (S.E corner of 2nd face): nPg+2
0206 vvI(1,1)=var(nPg+2);
1c4dce2004 Jean*0207 for l=0:2,
323aa12fe2 Jean*0208 xxI(1+rem(l,2),1+fix(l/2))=xx2(ncp,2,2+2*l);
0209 yyI(1+rem(l,2),1+fix(l/2))=yy2(ncp,2,2+2*l);
0210 end
0211 xxI(2,2)=xxI(1,2); yyI(2,2)=yyI(1,2);
0212 [nbsf,S(nbsf)]=part_surf(nbsf,fac,xxI,yyI,vvI,1,2,1,2,c1,c2) ;
0213
0214 %- N pole:
4ec37fd829 Jean*0215 vvI(1,1)=var(1+nc/2+nc/2*nc+2*nc*nc);
862f2fe640 Jean*0216 xxI(1,:)=xMid-180; xxI(2,:)=xMid+180;
0217 yyI(:,1)=max(ycs(:)); yyI(:,2)=+90;
323aa12fe2 Jean*0218 [nbsf,S(nbsf)]=part_surf(nbsf,fac,xxI,yyI,vvI,1,2,1,2,c1,c2) ;
0219 %- S pole:
4ec37fd829 Jean*0220 vvI(1,1)=var(1+nc/2+nc/2*nc+5*nc*nc);
862f2fe640 Jean*0221 xxI(1,:)=xMid-180; xxI(2,:)=xMid+180;
0222 yyI(:,2)=min(ycs(:)); yyI(:,1)=-90;
1c4dce2004 Jean*0223 [nbsf,S(nbsf)]=part_surf(nbsf,fac,xxI,yyI,vvI,1,2,1,2,c1,c2) ;
323aa12fe2 Jean*0224 %--------------
04cd582eaa Jean*0225 set(S,'LineStyle','-','LineWidth',0.01);
0226 if rem(kEnv,2) > 0, set(S,'EdgeColor','none'); end
323aa12fe2 Jean*0227 hold off
040160fb38 Andr*0228 if shift == -1,
323aa12fe2 Jean*0229 axis(AxBx); fprintf(' Axis(Box): %i %i %i %i \n',AxBx);
040160fb38 Andr*0230 elseif shift == 1
323aa12fe2 Jean*0231 [L]=draw_coast(fac);
0232 set(L,'color',[1 0 1]);
0233 % set(L,'Color',[0 0 0],'LineWidth',2); % set(L,'LineStyle','-');
0234 axis(AxBx); fprintf(' Axis(Box): %i %i %i %i \n',AxBx);
0235 else
0236 m_proj('Equidistant Cylindrical','lat',90,'lon',[-180+shift 180+shift])
0237 %m_proj('Equidistant Cylindrical','lat',[-0 60],'lon',[-180+shift 180+shift])
0238 m_coast('color',[0 0 0]);
0239 m_grid('box','on')
0240 end
0241
0242 %--
4ec37fd829 Jean*0243 if cbV < 2, bV=fix(cbV); moveHV_colbar([10+bV*2.2 10 7-5*bV 7+2*bV]/10,cbV); end
1c4dce2004 Jean*0244 if mnV < MxV & kEnv < 2,
e8997c8981 Jean*0245 ytxt=min(1,cbV);
040160fb38 Andr*0246 if shift == 1 | shift == -1,
4ec37fd829 Jean*0247 pos=get(gca,'position');
0248 xtxt=mean(AxBx(1:2)) ; ytxt=AxBx(3)-(AxBx(4)-AxBx(3))*(16-pos(4)*7.4)/100;
0249 %fprintf(' pp= %9.6f , ytxt= %7.2f\n',pos(4),ytxt);
0250 else xtxt=0; ytxt=(30*cbV-120)*fac;
323aa12fe2 Jean*0251 end
4ec37fd829 Jean*0252 text(xtxt*fac,ytxt*fac,sprintf('min,Max= %9.5g , %9.5g', mnV, MxV))
323aa12fe2 Jean*0253 else fprintf('Uniform field: min,Max= %g \n',MxV); end
0254 return
0255 %----------------------------------------------------------------------
0256 function [nbsf,S]=part_surf(nbsf,fac,xx,yy,vv,i1,i2,j1,j2,c1,c2)
0257 S=surf(fac*xx(i1:i2,j1:j2),fac*yy(i1:i2,j1:j2), ...
0258 zeros(i2+1-i1,j2+1-j1),vv(i1:i2,j1:j2)) ;
0259 if c1 < c2, caxis([c1 c2]) ; end
0260 %set(S,'LineStyle','-','LineWidth',0.01); %set(S,'EdgeColor','none');
0261 %set(S,'Clipping','off');
0262 %get(S) ;
0263 %nbsf = nbsf+1 ; if nbsf == 1 ; hold on ; view(0,90) ; end
0264 nbsf = nbsf+1 ; if nbsf == 1 ; hold on ; view(0,90) ; end
0265 %axis([-180 180 60 90]) ; % work only without coast-line
0266 %axis([-180 180 -90 -60]) ; % work only without coast-line
0267 %axis([30 60 -45 -25]) ; % work only without coast-line
0268 return
0269 %----------------------------------------------------------------------
0270
