Warning, /utils/matlab/cs_grid/bk_line/def_API_msk.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
323aa12fe2 Jean*0001
0002 krd=1; kpr=2; kgr=0; kwr=1;
0003 %krd=0; kpr=0; kgr=1; kwr=0; % <- execute a 2nd time & draw some plot
0004 %set_axis
d3e0534631 Jean*0005 %krd=1; kpr=3; kgr=3; kwr=1; % <- generate Continent Mask
323aa12fe2 Jean*0006
0007 if krd == 1,
0c7ba451ba Jean*0008 %- set ncdf=1 to load MNC (NetCDF) grid-files ;
0009 % or ncdf=0 to load MDS (binary) grid-files :
0010 ncdf=0;
0011 %rac='/home/jmc/grid_cs32/';
323aa12fe2 Jean*0012 rac='grid_files/';
0c7ba451ba Jean*0013 G=load_grid(rac,ncdf);
0014 xcs=G.xC; ycs=G.yC; xcg=G.xG; ycg=G.yG; arc=G.rAc;
0015 mskC=G.hFacC; mskW=G.hFacW; mskS=G.hFacS;
323aa12fe2 Jean*0016 mskC=min(1,mskC); mskC=ceil(mskC);
0017 mskW=min(1,mskW); mskW=ceil(mskW);
0018 mskS=min(1,mskS); mskS=ceil(mskS);
0019 ncx=size(mskC,1); nc=size(mskC,2); ncp=nc+1;
0020 xc1=reshape(xcs,1,ncx*nc);
0021 yc1=reshape(ycs,1,ncx*nc);
0022 xg1=reshape(xcg,1,ncx*nc);
0023 yg1=reshape(ycg,1,ncx*nc);
0024 mk1=reshape(mskC(:,:,1),1,ncx*nc);
0025 %xg6=split_Z_cub(xcg);
0026 %yg6=split_Z_cub(ycg);
0027 xc6=split_C_cub(xcs,1);
0028 yc6=split_C_cub(ycs,1);
0029 end
0030
d3e0534631 Jean*0031 if kpr == 1 | kpr == 2,
323aa12fe2 Jean*0032 %- define basin mask for Tracer pt: =1: Atlantic ; =2: Indian ; =3: Pacific
0033 [xcb,xPA,yPA,xAI,yAI,xIP,yIP]=line_sep(yc1);
0034 mbs=3*ones(1,ncx*nc);
0035 mbs(find( xc1>xcb(1,:) & xc1<xcb(2,:)) )=1;
0036 mbs(find( xc1>xcb(2,:) & xc1<xcb(3,:)) )=2;
0037 bas=reshape(mbs,ncx,nc);
0038
0039 %- define basin mask for U,V pt: Atlantic=1 ; Indian=2 ; Pacific=3
0040 mskBasU=zeros(ncx,nc,3); mskBasV=zeros(ncx,nc,3); mskBasC=zeros(ncx,nc,3);
0041 varU=zeros(nc,nc);varV=zeros(nc,nc);
0042 msk6=split_C_cub(bas,1);
0043 for b=1:3,
0044 mskBasC(:,:,b)=abs(b-bas);
0045 mskB=abs(b-msk6); mskB=1-min(1,mskB);
0046 for k=1:6,
0047 i1=1+nc*(k-1);i2=nc*k;
0048 var=mskB(1:ncp,2:ncp,k);
0049 varU(1:nc,:)=max(var(1:nc,:),var(2:ncp,:));
0050 mskBasU(i1:i2,:,b)=varU(1:nc,:);
0051 var=mskB(2:ncp,1:ncp,k);
0052 varV(:,1:nc)=max(var(:,1:nc),var(:,2:ncp));
0053 mskBasV(i1:i2,:,b)=varV(1:nc,:);
0054 end
d3e0534631 Jean*0055 end
323aa12fe2 Jean*0056 mskBasC=1-min(1,mskBasC);
0057
0058 end
0059
0060 if kpr==2,
0061
0062 %- find list of separation points :
0063 mskU=zeros(ncx,nc); mskV=zeros(ncx,nc);
0064 msk6=msk6+max(2,msk6)-2;
0065 for k=1:6,
0066 i1=1+nc*(k-1);i2=nc*k;
0067 var=msk6(1:ncp,2:ncp,k);
0068 varU(1:nc,:)=var(2:ncp,:)-var(1:nc,:);
0069 mskU(i1:i2,:)=varU(1:nc,:);
0070 var=msk6(2:ncp,1:ncp,k);
0071 varV(:,1:nc)=var(:,2:ncp)-var(:,1:nc);
0072 mskV(i1:i2,:)=varV(1:nc,:);
0073 end
0074 %mskU=mskU.*mskW(:,:,1); mskV=mskV.*mskS(:,:,1);
0075 mskU=reshape(mskU,1,ncx*nc); mskV=reshape(mskV,1,ncx*nc);
f3103ba969 Jean*0076 %- Separation bs=1 : Atl/Ind ; bs=2 : Ind/Pac ; bs=3 : Pac/Atl.
0077 % Defining msk6 as above (msk6 =1: Atl, =2: Ind, =4: Pac)
0078 % allows to get separation number as gradient (=mskU,mskV) of msk6
0079 for bs=1:3,
0080 [I]=find(abs(mskU)==bs);
0081 Nu(bs)=length(I); Iu(bs,1:Nu(bs))=I;
323aa12fe2 Jean*0082 i0=rem(I-1,ncx);j1=1+fix((I-1)/ncx);
0083 i1=1+rem(i0,nc); k1=1+fix(i0/nc);
f3103ba969 Jean*0084 xx=zeros(2,Nu(bs));
0085 for i=1:Nu(bs),
0086 xx(1,i)=xc6(1+i1(i),1+j1(i),k1(i));
0087 xx(2,i)=xc6( i1(i), 1+j1(i),k1(i));
0088 ySepU(bs,i)=(yc6(1+i1(i),1+j1(i),k1(i))+yc6(i1(i),1+j1(i),k1(i)))/2;
0089 end
0090 var=reshape(xcg,[1 ncx*nc]);
0091 xx(1,:)=var(I)-180+rem(xx(1,:)-var(I)+3*180,360);
0092 xx(2,:)=var(I)-180+rem(xx(2,:)-var(I)+3*180,360);
0093 xSepU(bs,[1:Nu(bs)])=(xx(1,:)+xx(2,:))/2;
0094 [J]=find(abs(mskV)==bs);
0095 Nv(bs)=length(J); Iv(bs,1:Nv(bs))=J;
323aa12fe2 Jean*0096 i0=rem(J-1,ncx);j1=1+fix((J-1)/ncx);
0097 i1=1+rem(i0,nc); k1=1+fix(i0/nc);
f3103ba969 Jean*0098 xx=zeros(2,Nv(bs));
0099 for i=1:Nv(bs),
0100 xx(1,i)=xc6(1+i1(i),1+j1(i),k1(i));
0101 xx(2,i)=xc6(1+i1(i), j1(i), k1(i));
0102 ySepV(bs,i)=(yc6(1+i1(i),1+j1(i),k1(i))+yc6(1+i1(i),j1(i),k1(i)))/2;
323aa12fe2 Jean*0103 end
f3103ba969 Jean*0104 xx(1,:)=var(J)-180+rem(xx(1,:)-var(J)+3*180,360);
0105 xx(2,:)=var(J)-180+rem(xx(2,:)-var(J)+3*180,360);
0106 xSepV(bs,[1:Nv(bs)])=(xx(1,:)+xx(2,:))/2;
323aa12fe2 Jean*0107 %----
f3103ba969 Jean*0108 fprintf('Sep. Line %i , Open Pts : %i %i \n',bs,Nu(bs),Nv(bs));
323aa12fe2 Jean*0109 end
0110 end
0111
d3e0534631 Jean*0112 if kpr == 3,
0113 %-- create a mask for each set of connected continents:
0114 % 1 = Asia+Europe+Africa ; 2 = Americas ; 2 = Antarctica ; 4 = Australia
0115 %- start from 1 point:
0116 x0cont=[60 -70 140 140]; nbCont=length(x0cont);
0117 y0cont=[40 0 -85 -25];
0118 ContinC=zeros(ncx,nc); locMskZ=zeros(ncp+1,ncp+1,6);
0119 for k=1:nbCont,
0120 locMskC=1.-mskC(:,:,1);
0121 fprintf('Cont. k= %i , Starting Point: %7.3f , %7.3f\n',k,x0cont(k),y0cont(k));
0122 dd=abs(xcs-x0cont(k))+abs(ycs-y0cont(k));
0123 [I J]=find(dd==min(dd(:)));
0124 for i=1:length(I),
0125 if mskC(I(i),J(i),1)==0, locMskC(I(i),J(i))=2; end
0126 end
0127 if max(locMskC(:)) < 2,
0128 error 'Starting point in water => need to pick an other point'
0129 end
0130 it=0; nSum=0; mSum=sum(locMskC(:));
0131 fprintf(' it=%i , mSum= %i\n',it,mSum);
0132 % for it=1:1,
0133 while mSum > nSum,
0134 it=it+1; nSum=mSum;
0135 msk2=split_C_cub(locMskC,2);
0136 msk2(1,1,:)=0; msk2(ncp+1,ncp+1,:)=0; msk2(ncp+1,1,:)=0; msk2(1,ncp+1,:)=0;
0137 for n=1:6,
0138 tmp=min(msk2(:,:,n),1);
0139 var=msk2(:,:,n)-1; var=max(var,0);
0140 var(2:ncp,2:ncp)=var(2:ncp,2:ncp) ...
0141 +var(1:nc,2:ncp)+var(3:ncp+1,2:ncp) ...
0142 +var(2:ncp,1:nc)+var(2:ncp,3:ncp+1) ...
0143 +var(1:nc,1:nc)+var(3:ncp+1,3:ncp+1) ...
0144 +var(1:nc,3:ncp+1)+var(3:ncp+1,1:nc) ...
0145 ; var=var.*tmp;
0146 var=min(var,1);
0147 msk2(:,:,n)=msk2(:,:,n)+var;
0148 end
0149 msk2=min(msk2,2);
0150 var=msk2(2:ncp,2:ncp,:);locMskC=reshape(permute(var,[1 3 2]),[ncx nc]);
0151 mSum=sum(locMskC(:));
0152 fprintf(' it=%i , mSum= %i\n',it,mSum);
0153 %- exit infinite loop if something is wrong,
0154 if it > ncx,
0155 mSum=nSum; fprintf('Exit loop (k=%i) at it= %i\n',k,it);
0156 end
0157 end
0158 locMskC=max(locMskC,1)-1;
0159 msk2=split_C_cub(locMskC,2);
0160 msk2(1,1,:)=0; msk2(ncp+1,ncp+1,:)=0; msk2(ncp+1,1,:)=0; msk2(1,ncp+1,:)=0;
0161 for n=1:6,
0162 var=msk2(:,:,n);
0163 var(2:ncp+1,2:ncp+1)=var(2:ncp+1,2:ncp+1)+var(1:ncp,1:ncp) ...
0164 +var(1:ncp,2:ncp+1)+var(2:ncp+1,1:ncp);
0165 var=min(var,1);
0166 msk2(:,:,n)=var;
0167 end
0168 ContinC=ContinC+k*locMskC;
0169 locMskZ=locMskZ+k*msk2;
0170 end
0171 ContinZ=zeros(ncx*nc+2,1);
0172 var=locMskZ(2:ncp,2:ncp,:);
0173 ContinZ(1:ncx*nc)=reshape(permute(var,[1 3 2]),[ncx*nc 1]);
0174 ContinZ(ncx*nc+1)=locMskZ(2,ncp+1,1);
0175 ContinZ(ncx*nc+2)=locMskZ(ncp+1,2,2);
0176 end
0177
323aa12fe2 Jean*0178 if kwr==1,
d3e0534631 Jean*0179 if kpr==1 | kpr==2,
0180 namf='maskC_bas.bin';
0181 fid=fopen(namf,'w','b'); fwrite(fid,mskBasC,'real*4'); fclose(fid);
0182 fprintf(['Write Atl-Ind-Pac Bass. mask on file:',namf,' : O.K. \n']);
0183 namf='maskW_bas.bin';
0184 fid=fopen(namf,'w','b'); fwrite(fid,mskBasU,'real*4'); fclose(fid);
0185 fprintf(['Write Atl-Ind-Pac Bass. mask on file:',namf,' : O.K. \n']);
0186 namf='maskS_bas.bin';
0187 fid=fopen(namf,'w','b'); fwrite(fid,mskBasV,'real*4'); fclose(fid);
0188 fprintf(['Write Atl-Ind-Pac Bass. mask on file:',namf,' : O.K. \n']);
0189 end
0190 if kpr==2,
0191 namf='open_basins_section';
0192 save(namf,'Nu','Nv','Iu','Iv');
0193 fprintf(['write Basin connection on file:',namf,'.mat : O.K. \n']);
0194 end
0195 if kpr==3,
0196 namf='mask_Cont.bin';
0197 fid=fopen(namf,'w','b');
0198 fwrite(fid,ContinZ,'real*4');
0199 fwrite(fid,ContinC,'real*4');
0200 fclose(fid);
0201 fprintf('Write %i Continents mask on file: %s : O.K. \n',nbCont,namf);
0202 end
323aa12fe2 Jean*0203 end
0204
d3e0534631 Jean*0205 if kgr==1,
323aa12fe2 Jean*0206 %yy1=[-89.5:1:89.5]; [x3b]=line_sep(yy1);
0207 figure(1); clf;
f3103ba969 Jean*0208 shift=-1; cbV=1; ccB=[0 0]; AxBx=[-182 182 -91 91];
0209 var=bas ; ccB=[-1.5 5.5]; var(find(mskC(:,:,1)==0))=NaN;
323aa12fe2 Jean*0210 grph_CS(var,xcs,ycs,xcg,ycg,ccB(1),ccB(2),shift,cbV,AxBx);
0211 hold on ;
0212 %plot(x3b(1,:),yy1,'r-');
0213 %plot(x3b(2,:),yy1,'r-');
0214 %plot(x3b(3,:),yy1,'r-');
0215 plot(xPA,yPA,'*b-');
0216 plot(xAI,yAI,'*b-');
0217 plot(xIP,yIP,'*b-');
f3103ba969 Jean*0218 for bs=1:3,
0219 plot(xSepU(bs,1:Nu(bs)),ySepU(bs,1:Nu(bs)),'r*');
0220 plot(xSepV(bs,1:Nv(bs)),ySepV(bs,1:Nv(bs)),'ro');
323aa12fe2 Jean*0221 end
0222 hold off
f3103ba969 Jean*0223 %figure(2); clf;
0224 %for b=1:3,
0225 % subplot(310+b);
0226 % var=mskBasC(:,:,b); titv='mskBasC';
0227 % %var=mskBasU(:,:,b); titv='mskBasU';
0228 % %var=mskBasV(:,:,b); titv='mskBasV';
0229 % grph_CS(var,xcs,ycs,xcg,ycg,ccB(1),ccB(2),shift,cbV,AxBx);
0230 % title([titv,' : nBas= ',int2str(b)]);
0231 %end
0232 figure(3); clf; bs=3;
0233 shift=-1; cbV=1; ccB=[0 0]; AxBx=[-182 182 -91 91]; ccB=[-1 1]*5;
0234 subplot(211);
0235 var=reshape(mskU,[ncx nc]) ;
0236 grph_CS(var,xcs,ycs,xcg,ycg,ccB(1),ccB(2),shift,cbV,AxBx);
0237 hold on; plot(xSepU(bs,1:Nu(bs)),ySepU(bs,1:Nu(bs)),'m*'); hold off
0238 title(['mskU and SepU for separation bs= ',int2str(bs)]);
0239 subplot(212);
0240 var=reshape(mskV,[ncx nc]) ;
0241 grph_CS(var,xcs,ycs,xcg,ycg,ccB(1),ccB(2),shift,cbV,AxBx);
0242 hold on; plot(xSepV(bs,1:Nv(bs)),ySepV(bs,1:Nv(bs)),'mo'); hold off
0243 title(['mskV and SepV for separation bs= ',int2str(bs)]);
0244
323aa12fe2 Jean*0245 end
0246
d3e0534631 Jean*0247 %----
0248 if kgr==3,
0249 figure(2); clf;
0250 shift=-1; cbV=1; ccB=[0 0]; AxBx=[-180 180 -90 90];
0251 ccB=[-1 nbCont+1];
0252 subplot(211)
0253 var=ContinC;
0254 var(find(mskC(:,:,1)==1))=NaN;
0255 grph_CS(var,xcs,ycs,xcg,ycg,ccB(1),ccB(2),shift,cbV,AxBx);
0256 subplot(212)
0257 var=ContinZ; ccB=[-1 nbCont+1];
0258 grph_CSz(var,xcs,ycs,xcg,ycg,ccB(1),ccB(2),shift,cbV,AxBx);
0259 end
0260
323aa12fe2 Jean*0261 return
