Warning, /utils/matlab/cs_grid/bk_line/grph_CS.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
323aa12fe2 Jean*0001 function [fac]=grph_CS(var,xcs,ycs,xcg,ycg,c1,c2,shift,cbV,AxBx,kEnv)
1c4dce2004 Jean*0002 % [fac]=grph_CS(var,xcs,ycs,xcg,ycg,c1,c2,shift[,cbV,AxBx,kEnv]) :
497e5a9d34 Jean*0003 % produce a flat plot of the cube_sphere field "var",
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 % 1 : No shift but draw coast-line calling draw_coast
323aa12fe2 Jean*0010 % else : draw coast-line (using m_proj) shifted by "shift" degree.
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.
040160fb38 Andr*0013 % AxBx = do axis(AxBx) to zoom in Box "AxBx" ; only if shift=-1,1 ;
323aa12fe2 Jean*0014 %-----------------------
497e5a9d34 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
323aa12fe2 Jean*0030 if nargin < 11, kEnv=0 ; end
0031
0032 %------------------------------
4ec37fd829 Jean*0033 n1h=size(var,1); n2h=size(var,2);
0034 if n1h == 6*n2h, nc=n2h;
0035 elseif n1h*6 == n2h, nc=n1h;
0036 else
0037 error([' input var size: ',int2str(n1h),' x ',int2str(n2h),' does not fit regular cube !']);
0038 end
0039 ncp=nc+1 ; nPg=nc*nc*6;
323aa12fe2 Jean*0040 MxV=min(min(var));
0041 mnV=max(max(var));
451cd38a73 Jean*0042 if shift == 1 | shift == -1,
c8c5756e04 Jean*0043 xxt=reshape(xcs,[nPg 1]); yyt=reshape(ycs,[nPg 1]);
0044 tmp=(xxt-AxBx(1)).*(AxBx(2)-xxt); [I]=find(tmp > 0 );
0045 tmp=(yyt(I)-AxBx(3)).*(AxBx(4)-yyt(I)); [J]=find(tmp > 0 );
0046 IJ=I(J); tmp=var(IJ); mnV=min(tmp); MxV=max(tmp);
323aa12fe2 Jean*0047 end
0048 %------------
0049 if c1 >= c2
0050 mb=(MxV-mnV)*0.01;
0051 c1=mnV+mb*c1;
0052 c2=MxV+mb*c2;
1c4dce2004 Jean*0053 if c1*c2 < 0
323aa12fe2 Jean*0054 c2=max(abs([c1 c2]));
0055 c1=-c2;
0056 end
0057 fprintf('min,max %8.3e %8.3e Cmin,max %8.3e %8.3e \n',mnV,MxV,c1,c2)
0058 end
0059 %------------------------------
0060 %figure(1);
040160fb38 Andr*0061 if shift ~= 1 & shift ~= -1,
0062 fac=pi/180.;
323aa12fe2 Jean*0063 else
0064 fac=1. ;
0065 end
0066 %---
0067 nbsf = 0 ; ic = 0 ; jc = 0 ;
497e5a9d34 Jean*0068 nPx=prod(size(xcg)); nPy=prod(size(ycg));
0069 if nPx == nPg & nPy == nPg,
4ec37fd829 Jean*0070 %- when stored in long-vector, use "compact" convention (i.e., 1 face after the other)
0071 if n2h == nc,
0072 xcg=permute(reshape(xcg,[nc 6 nc]),[1 3 2]);
0073 ycg=permute(reshape(ycg,[nc 6 nc]),[1 3 2]);
0074 end
497e5a9d34 Jean*0075 xcg=reshape(xcg,[nPg 1]); ycg=reshape(ycg,[nPg 1]);
323aa12fe2 Jean*0076 %- add the 2 missing corners:
0077 %fprintf(' Local version of grph_CS : ---------------------------------- \n');
4ec37fd829 Jean*0078 xcg(nPg+1)=xcg(1); ycg(nPg+1)=ycg(1+2*nc*nc);
0079 xcg(nPg+2)=xcg(1+3*nc*nc); ycg(nPg+2)=ycg(1);
497e5a9d34 Jean*0080 elseif nPx ~= nPg+2 | nPy ~= nPg+2,
0081 error([' wrong xcg,ycg dimensions : ',int2str(nPx),' , ',int2str(nPy)]);
0082 end
6b4f212f67 Jean*0083 [x6c]=split_C_cub(xcs,0);
323aa12fe2 Jean*0084 [xx2]=split_Z_cub(xcg);
0085 [yy2]=split_Z_cub(ycg);
0086 %---
4ec37fd829 Jean*0087 if n2h == nc,
0088 var=permute(reshape(var,[nc 6 nc]),[1 3 2]);
0089 else
0090 var=reshape(var,[nc nc 6]);
0091 end
323aa12fe2 Jean*0092 for n=1:6,
0093 %if n < 5 & n > 2,
0094 if n < 7,
0095 %--------------------------------------------------------
6b4f212f67 Jean*0096 vv1=zeros(ncp,ncp) ; xx1=vv1 ; yy1=vv1 ;
323aa12fe2 Jean*0097 vv1(1:nc,1:nc)=var(:,:,n);
6b4f212f67 Jean*0098 xx1=xx2(:,:,n); yy1=yy2(:,:,n); xxc=x6c(:,:,n);
862f2fe640 Jean*0099 % if xx1(ncp,1) < xMid-300. ; xx1(ncp,1)=xx1(ncp,1)+360. ; end
0100 % if xx1(1,ncp) < xMid-300. ; xx1(1,ncp)=xx1(1,ncp)+360. ; end
323aa12fe2 Jean*0101 %------------
6b4f212f67 Jean*0102 %- use jump in grid-cell center longitude to decide where to cut 1 face:
0103 cutFace=0;
0104 dxi=xxc(2:nc,:)-xxc(1:nc-1,:); dxImx=max(abs(dxi(:)));
0105 dxj=xxc(:,2:nc)-xxc(:,1:nc-1); dxJmx=max(abs(dxj(:)));
0106 if dxImx > xJump & dxJmx > xJump, cutFace=3;
0107 elseif dxImx > xJump, cutFace=1;
0108 elseif dxJmx > xJump, cutFace=2; end
0109 if dBug > 0,
0110 fprintf(' face # %i , Max dxI,dxJ = %8.3f , %8.3f',n,dxImx,dxJmx);
0111 if cutFace > 0, fprintf(' ; cutFace= %i',cutFace); end
0112 fprintf('\n');
0113 end
0114 if cutFace == 3,
0115 fprintf(' Jump in both i & j not implemented ==> skip face # %i\n',n);
0116 %------------
0117 elseif cutFace == 2,
0118 [I,J]=find( abs(dxj) > xJump );
0119 if dBug > 1, for l=1:length(I), fprintf(' i,j= %2i, %2i \n',I(l),J(l)); end; end
0120 if min(J) == max(J),
0121 jc =J(1); jp=jc+1;
0122 if dBug > 0, fprintf('--> cut Face @ jc,jp = %3i,%3i\n',jc,jp); end
0123 % cut the face in 2 parts (1: j in [1 jp] ; 2: j in [jp ncp]) and plot separately
0124 for lp=1:2,
0125 if lp == 1,
0126 j1=1; j2=jp;
0127 else
0128 j1=jp; j2=ncp;
0129 end
0130 tmp=xxc(:,j1:j2-1); xLoc=mean(tmp(:));
0131 if dBug > 0, fprintf(' p.%i, %2i -> %2i : %8.3f %8.3f %8.3f\n', ...
0132 lp,j1,j2-1,min(tmp(:)),xLoc,max(tmp(:))); end
0133 if xLoc > xMid,
0134 %-- case where Xc jump from > 180 to < -180 when j goes from jc to jc+1 :
0135 [I]=find(xx1(:,jp) < xMid-xJmp2); xx1(I,jp)=xx1(I,jp)+360.;
0136 % Pole longitude is arbitrary: set to +90 to get a nicer plot:
0137 xx1(find( abs(yy1) > 90-epsil ))=xMid+90;
0138 else
0139 %-- case where Xc jump from < -180 to > 180 when j goes from jc to jc+1 :
0140 [I]=find(xx1(:,jp) > xMid+xJmp2); xx1(I,jp)=xx1(I,jp)-360.;
0141 % Pole longitude is arbitrary: set to -90 to get a nicer plot:
0142 xx1(find( abs(yy1) > 90-epsil ))=xMid-90;
0143 end
0144 [nbsf,S(nbsf)]=part_surf(nbsf,fac,xx1,yy1,vv1,1,ncp,j1,j2,c1,c2) ;
0145 end
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 [ip ncp]) and plot separately
0157 for lp=1:2,
0158 if lp == 1,
0159 i1=1; i2=ip;
0160 else
0161 i1=ip; i2=ncp;
0162 end
0163 tmp=xxc(i1:i2-1,:); xLoc=mean(tmp(:));
0164 if dBug > 0, fprintf(' p.%i, %2i -> %2i : %8.3f %8.3f %8.3f\n', ...
0165 lp,i1,i2-1,min(tmp(:)),xLoc,max(tmp(:))); end
0166 if xLoc > xMid,
0167 %-- case where X jump from > 180 to < -180 when i goes from ic to ic+1 :
0168 [J]=find(xx1(ip,:) < xMid-xJmp2); xx1(ip,J)=xx1(ip,J)+360.;
0169 % Pole longitude is arbitrary: set to +90 to get a nicer plot:
0170 xx1(find( abs(yy1) > 90-epsil ))=xMid+90;
0171 else
0172 %-- case where X jump from < -180 to > 180 when i goes from ic to ic+1 :
0173 [J]=find(xx1(ip,:) > xMid+xJmp2); xx1(ip,J)=xx1(ip,J)-360.;
0174 % Pole longitude is arbitrary: set to -90 to get a nicer plot:
0175 xx1(find( abs(yy1) > 90-epsil ))=xMid-90;
0176 end
0177 [nbsf,S(nbsf)]=part_surf(nbsf,fac,xx1,yy1,vv1,i1,i2,1,ncp,c1,c2) ;
862f2fe640 Jean*0178 %-
6b4f212f67 Jean*0179 end
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
6b4f212f67 Jean*0190 %--------------------------------------------------------
1c4dce2004 Jean*0191 end ; end
0192 set(S,'LineStyle','-','LineWidth',0.01);
323aa12fe2 Jean*0193 if rem(kEnv,2) > 0, set(S,'EdgeColor','none'); end
0194 hold off
040160fb38 Andr*0195 if shift == -1,
1c4dce2004 Jean*0196 axis(AxBx); fprintf(' Axis(Box): %i %i %i %i \n',AxBx);
0197 elseif shift == 1,
323aa12fe2 Jean*0198 [L]=draw_coast(fac);
0199 set(L,'color',[1 0 1]);
0200 % set(L,'Color',[0 0 0],'LineWidth',2); % set(L,'LineStyle','-');
1c4dce2004 Jean*0201 axis(AxBx); fprintf(' Axis(Box): %i %i %i %i \n',AxBx);
323aa12fe2 Jean*0202 else
1c4dce2004 Jean*0203 m_proj('Equidistant Cylindrical','lat',90,'lon',[-180+shift 180+shift])
0204 %m_proj('Equidistant Cylindrical','lat',[-0 60],'lon',[-180+shift 180+shift])
323aa12fe2 Jean*0205 m_coast('color',[0 0 0]);
0206 m_grid('box','on')
0207 end
0208
0209 %--
c8c5756e04 Jean*0210 if cbV < 2, bV=fix(cbV); moveHV_colbar([10+bV*2.2 10 7-5*bV 7+2*bV]/10,cbV); end
323aa12fe2 Jean*0211 if mnV < MxV & kEnv < 2,
c8c5756e04 Jean*0212 ytxt=min(1,fix(cbV));
040160fb38 Andr*0213 if shift == 1 | shift == -1,
4ec37fd829 Jean*0214 pos=get(gca,'position');
0215 xtxt=mean(AxBx(1:2)) ; ytxt=AxBx(3)-(AxBx(4)-AxBx(3))*(16-pos(4)*7.4)/100;
0216 %fprintf(' pp= %9.6f , ytxt= %7.2f\n',pos(4),ytxt);
1c4dce2004 Jean*0217 else
323aa12fe2 Jean*0218 xtxt=60 ; ytxt=30*ytxt-145 ;
0219 %xtxt= 0 ; ytxt=30*ytxt-120 ;
0220 end
1c4dce2004 Jean*0221 %fprintf('min,Max= %9.5g , %9.5g (xtxt,ytxt= %f %f)\n',mnV,MxV,xtxt,ytxt);
0222 text(xtxt*fac,ytxt*fac,sprintf('min,Max= %9.5g , %9.5g', mnV, MxV))
323aa12fe2 Jean*0223 %set(gca,'position',[-.1 0.2 0.8 0.75]); % xmin,ymin,xmax,ymax in [0,1]
1c4dce2004 Jean*0224 elseif mnV < MxV, fprintf('min,Max= %9.5g , %9.5g \n',mnV,MxV);
323aa12fe2 Jean*0225 else fprintf('Uniform field: min,Max= %g \n',MxV); end
0226 return
0227 %----------------------------------------------------------------------
0228 function [nbsf,S]=part_surf(nbsf,fac,xx,yy,vv,i1,i2,j1,j2,c1,c2)
0229 S=surf(fac*xx(i1:i2,j1:j2),fac*yy(i1:i2,j1:j2), ...
0230 zeros(i2+1-i1,j2+1-j1),vv(i1:i2,j1:j2)) ;
0231 %shading flat ; %- or faceted (=default) or interp
0232 if c1 < c2, caxis([c1 c2]) ; end
0233 %set(S,'LineStyle','-','LineWidth',0.01); set(S,'EdgeColor','none');
0234 %set(S,'Clipping','off');
0235 %get(S) ;
0236 %nbsf = nbsf+1 ; if nbsf == 1 ; hold on ; view(0,90) ; end
0237 nbsf = nbsf+1 ; if nbsf == 1 ; hold on ; view(0,90) ; end
0238 %axis([-180 180 60 90]) ; % work only without coast-line
0239 %axis([-180 180 -90 -60]) ; % work only without coast-line
0240 %axis([30 60 -45 -25]) ; % work only without coast-line
0241 return
0242 %----------------------------------------------------------------------
0243
