| 1 |
gforget |
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 |
|
|
|
Parent Directory
| 
Revision Log
| 
Revision Graph