matlab_class/gcmfaces_calc/gcmfaces_close_basin.m

function [msk]=gcmfaces_close_basin(basinEdge,landMsk,doLargerBasin);
%object : compute closed basin mask from a specificed edge and a "land" mask
%         By design basinEdge is required to split the global ocean into two and 
%         only two basins. Successive calls allow for numerous basin definitions.
%
%input :  - basinEdge is an index (standard cases; defined by great circle arcs
%           between pairs of locations; see below) or a 2D map (interactive case)
%           - test cases : basinEdge = 0 to -3 (included hereafter)
%           - user cases : basinEdge > 0 (for user to add hereafter)
%           - interactive case : basinEdge is a 2D map of 1s and NaNs 
%             where NaNs delineate the "wet" egde of the basin.
%         - landMsk is the "land" mask (1s & NaNs; mygrid.mskC(:,:,1) default)
%           basinEdge values where landMsk is NaN are of no consequence. 
%         - doLargerBasin is 0 (default) or 1. If set to 1 then output 
%           larger basin mask, otherwise output smaller basin mask. The ocean 
%           edge points are always counted as part of the output basin mask. 
%
%output : msk is the closed basin mask, with 1s over the basin of interest, 
%         0s over the rest of the ocean, and NaNs consistent with landMsk.
%
%notes : - unless used in interactive mode the ocean basin edge is based on 
%          great circle arcs. That is usually perfectly adequate. The most 
%          notable exception is when one wants to delimit a basin by latitude lines, 
%          but those are easy to do in interactive mode. The NOT implemented general 
%          approach would use small circle arcs rather than great circle arcs.
%        - the test cases below demonstrate the algorithm using great circles (0), 
%          the algorithm error cases (-1, -2) and the limitation of great circles (-3)

gcmfaces_global;

if isempty(whos('basinEdge')); basinEdge=0; end;
if isempty(whos('landMsk')); landMsk=mygrid.mskC(:,:,1); end;
if isempty(whos('doLargerBasin')); doLargerBasin=0; end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%BEGINING OF EDGE DEFINITION%%%%%%%%%%%%%%%%%%%%%%%%%%%

if isnumeric(basinEdge)&length(basinEdge(:))==1;
%standard cases: to be extended by user

if basinEdge==0;
lonPairs=[[  10   19];[  19   60];[  60  100];[ 100  110]];
latPairs=[[  60   80];[  80   81];[  81   80];[  80   60]];
%a case that matches v4_basin_one
elseif basinEdge==-1;
lonPairs=[[ -70   20];[  20  110];[ 110  200];[ 200  -70]];
latPairs=[[ -30  -30];[ -30  -30];[ -30  -30];[ -30  -30]];
%a case that fails because it delimits three basins
elseif basinEdge==-2;
lonPairs=[[ -70   20];[  20  145]];
latPairs=[[ -35  -35];[ -35  -35]];
%a case that fails because it does not close basins
elseif basinEdge==-3;
lonPairs=[[ -70   20];[  20  145];[ 145  -70]];
latPairs=[[ -35  -35];[ -35  -35];[ -35  -35]];
%a case that does not fail, but illustrates that long great 
%circle arcs are not always the most appopriate basin edges 
%(a line of constant latitude would be better in this case)
else;
error('unknown basin');
end;

%close basin:
msk=landMsk;
msk0=msk;
for iy=1:size(lonPairs,1);
    %define great circle line by end points:
    lonPair=lonPairs(iy,:);
    latPair=latPairs(iy,:);
    %get line:
    line_cur=gcmfaces_lines_transp(lonPair,latPair,{'tmp'});
    %close line:
    msk(line_cur.mskCedge==1)=NaN;
end;


else;
%interactive cases: based on field provided by user

msk=landMsk;
msk0=msk;
msk(isnan(basinEdge))=NaN;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%END OF EDGE DEFINITION%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%compute the closed basin labels for each grid face:
msk1=msk;
msk=1*~isnan(msk);
for iF=1:msk.nFaces; 
  msk{iF}=bwlabel(msk{iF},4); 
end;

%remove potential aliases / account for domain periodicity:
msk2=msk;
test0=1;%do at least one faces sweep
while test0;
   test0=0;%unless something changes in the upcoming sweep we will stop
   tmpmsk=exch_T_N(msk);
   %list aliases
   list0=[];
   for iF=1:msk.nFaces;
     for iS=1:4;
       tmp1=msk{iF};
       if iS==1;     tmp2=tmpmsk{iF}(1:end-2,2:end-1);
       elseif iS==2; tmp2=tmpmsk{iF}(3:end,2:end-1);
       elseif iS==3; tmp2=tmpmsk{iF}(2:end-1,1:end-2); 
       elseif iS==4; tmp2=tmpmsk{iF}(2:end-1,3:end); 
       end;
       tmp3=find( (tmp1~=tmp2)&(tmp1>0)&(tmp2~=0));
       %note : tmp1~=tmp2 is an alias if the two other conditions apply
       %note : "land" (label 0) is not an alias (tmp1>0&tmp2~=0 test)
       %note : basins cannot be relabeld multiple times (tmp1>0 test) but
       %       relabeld basins can be tested against (tmp2~=0 rather than tmp2>0)
       list0=[list0;[tmp1(tmp3) tmp2(tmp3) iF+0*tmp3]];
     end;
   end;
   %not needed : list0=unique(list0,'rows');
   if ~isempty(list0); 
     %start with previously relabeld basin if any (needed to avoid conflicts)
     ii=find(list0(:,2)==min(list0(:,2))); ii=ii(1);
     iF=list0(ii,3);
     list1=list0(ii,1:2);
     tmp1=msk{iF};
     %since we will have relabeld at least one basin in this 
     %faces sweep, then another faces sweep will be needed :
     test0=1;
     %apply unique negative label
     if list1(2)<0;
       %carry over unique label from previously relabeld basin
       tmp4=list1(2);
     else;
       %define new unique negative label
       tmp4=min(min(msk)-1,-1);
     end;
     tmp1(tmp1==list1(1))=tmp4;
     msk{iF}=tmp1;
   end;
end;

%final list should have 2 elements (1 per basin)
list2=convert2array(msk);
list2=unique(list2(:));
list2=list2(find(~isnan(list2)&list2~=0));
if length(list2)<2; msk=msk1; fprintf('failure : basin was not closed properly\n'); return; end;
if length(list2)>2; fprintf('failure : more than two basins\n'); return; end;

%attribute edge and prepare output:
tmp1=sum(mygrid.RAC.*(msk==list2(1)));
tmp2=sum(mygrid.RAC.*(msk==list2(2)));
if     doLargerBasin&tmp1>=tmp2; excl=list2(2);
elseif doLargerBasin           ; excl=list2(1);
elseif               tmp1<=tmp2; excl=list2(2);
else                           ; excl=list2(1);
end;
msk3=msk;
msk=(msk~=excl).*mygrid.mskC(:,:,1);


1	function [msk]=gcmfaces_close_basin(basinEdge,landMsk,doLargerBasin);
2	%object : compute closed basin mask from a specificed edge and a "land" mask
3	% By design basinEdge is required to split the global ocean into two and
4	% only two basins. Successive calls allow for numerous basin definitions.
5	%
6	%input : - basinEdge is an index (standard cases; defined by great circle arcs
7	% between pairs of locations; see below) or a 2D map (interactive case)
8	% - test cases : basinEdge = 0 to -3 (included hereafter)
9	% - user cases : basinEdge > 0 (for user to add hereafter)
10	% - interactive case : basinEdge is a 2D map of 1s and NaNs
11	% where NaNs delineate the "wet" egde of the basin.
12	% - landMsk is the "land" mask (1s & NaNs; mygrid.mskC(:,:,1) default)
13	% basinEdge values where landMsk is NaN are of no consequence.
14	% - doLargerBasin is 0 (default) or 1. If set to 1 then output
15	% larger basin mask, otherwise output smaller basin mask. The ocean
16	% edge points are always counted as part of the output basin mask.
17	%
18	%output : msk is the closed basin mask, with 1s over the basin of interest,
19	% 0s over the rest of the ocean, and NaNs consistent with landMsk.
20	%
21	%notes : - unless used in interactive mode the ocean basin edge is based on
22	% great circle arcs. That is usually perfectly adequate. The most
23	% notable exception is when one wants to delimit a basin by latitude lines,
24	% but those are easy to do in interactive mode. The NOT implemented general
25	% approach would use small circle arcs rather than great circle arcs.
26	% - the test cases below demonstrate the algorithm using great circles (0),
27	% the algorithm error cases (-1, -2) and the limitation of great circles (-3)
28
29	gcmfaces_global;
30
31	if isempty(whos('basinEdge')); basinEdge=0; end;
32	if isempty(whos('landMsk')); landMsk=mygrid.mskC(:,:,1); end;
33	if isempty(whos('doLargerBasin')); doLargerBasin=0; end;
34
35	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%BEGINING OF EDGE DEFINITION%%%%%%%%%%%%%%%%%%%%%%%%%%%
36
37	if isnumeric(basinEdge)&length(basinEdge(:))==1;
38	%standard cases: to be extended by user
39
40	if basinEdge==0;
41	lonPairs=[[ 10 19];[ 19 60];[ 60 100];[ 100 110]];
42	latPairs=[[ 60 80];[ 80 81];[ 81 80];[ 80 60]];
43	%a case that matches v4_basin_one
44	elseif basinEdge==-1;
45	lonPairs=[[ -70 20];[ 20 110];[ 110 200];[ 200 -70]];
46	latPairs=[[ -30 -30];[ -30 -30];[ -30 -30];[ -30 -30]];
47	%a case that fails because it delimits three basins
48	elseif basinEdge==-2;
49	lonPairs=[[ -70 20];[ 20 145]];
50	latPairs=[[ -35 -35];[ -35 -35]];
51	%a case that fails because it does not close basins
52	elseif basinEdge==-3;
53	lonPairs=[[ -70 20];[ 20 145];[ 145 -70]];
54	latPairs=[[ -35 -35];[ -35 -35];[ -35 -35]];
55	%a case that does not fail, but illustrates that long great
56	%circle arcs are not always the most appopriate basin edges
57	%(a line of constant latitude would be better in this case)
58	else;
59	error('unknown basin');
60	end;
61
62	%close basin:
63	msk=landMsk;
64	msk0=msk;
65	for iy=1:size(lonPairs,1);
66	%define great circle line by end points:
67	lonPair=lonPairs(iy,:);
68	latPair=latPairs(iy,:);
69	%get line:
70	line_cur=gcmfaces_lines_transp(lonPair,latPair,{'tmp'});
71	%close line:
72	msk(line_cur.mskCedge==1)=NaN;
73	end;
74
75
76	else;
77	%interactive cases: based on field provided by user
78
79	msk=landMsk;
80	msk0=msk;
81	msk(isnan(basinEdge))=NaN;
82	end;
83
84	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%END OF EDGE DEFINITION%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
85
86	%compute the closed basin labels for each grid face:
87	msk1=msk;
88	msk=1*~isnan(msk);
89	for iF=1:msk.nFaces;
90	msk{iF}=bwlabel(msk{iF},4);
91	end;
92
93	%remove potential aliases / account for domain periodicity:
94	msk2=msk;
95	test0=1;%do at least one faces sweep
96	while test0;
97	test0=0;%unless something changes in the upcoming sweep we will stop
98	tmpmsk=exch_T_N(msk);
99	%list aliases
100	list0=[];
101	for iF=1:msk.nFaces;
102	for iS=1:4;
103	tmp1=msk{iF};
104	if iS==1; tmp2=tmpmsk{iF}(1:end-2,2:end-1);
105	elseif iS==2; tmp2=tmpmsk{iF}(3:end,2:end-1);
106	elseif iS==3; tmp2=tmpmsk{iF}(2:end-1,1:end-2);
107	elseif iS==4; tmp2=tmpmsk{iF}(2:end-1,3:end);
108	end;
109	tmp3=find( (tmp1~=tmp2)&(tmp1>0)&(tmp2~=0));
110	%note : tmp1~=tmp2 is an alias if the two other conditions apply
111	%note : "land" (label 0) is not an alias (tmp1>0&tmp2~=0 test)
112	%note : basins cannot be relabeld multiple times (tmp1>0 test) but
113	% relabeld basins can be tested against (tmp2~=0 rather than tmp2>0)
114	list0=[list0;[tmp1(tmp3) tmp2(tmp3) iF+0*tmp3]];
115	end;
116	end;
117	%not needed : list0=unique(list0,'rows');
118	if ~isempty(list0);
119	%start with previously relabeld basin if any (needed to avoid conflicts)
120	ii=find(list0(:,2)==min(list0(:,2))); ii=ii(1);
121	iF=list0(ii,3);
122	list1=list0(ii,1:2);
123	tmp1=msk{iF};
124	%since we will have relabeld at least one basin in this
125	%faces sweep, then another faces sweep will be needed :
126	test0=1;
127	%apply unique negative label
128	if list1(2)<0;
129	%carry over unique label from previously relabeld basin
130	tmp4=list1(2);
131	else;
132	%define new unique negative label
133	tmp4=min(min(msk)-1,-1);
134	end;
135	tmp1(tmp1==list1(1))=tmp4;
136	msk{iF}=tmp1;
137	end;
138	end;
139
140	%final list should have 2 elements (1 per basin)
141	list2=convert2array(msk);
142	list2=unique(list2(:));
143	list2=list2(find(~isnan(list2)&list2~=0));
144	if length(list2)<2; msk=msk1; fprintf('failure : basin was not closed properly\n'); return; end;
145	if length(list2)>2; fprintf('failure : more than two basins\n'); return; end;
146
147	%attribute edge and prepare output:
148	tmp1=sum(mygrid.RAC.*(msk==list2(1)));
149	tmp2=sum(mygrid.RAC.*(msk==list2(2)));
150	if doLargerBasin&tmp1>=tmp2; excl=list2(2);
151	elseif doLargerBasin ; excl=list2(1);
152	elseif tmp1<=tmp2; excl=list2(2);
153	else ; excl=list2(1);
154	end;
155	msk3=msk;
156	msk=(msk~=excl).*mygrid.mskC(:,:,1);
157
158