| 1 |
function [msk]=v4_basin_one(numBasin); |
| 2 |
%object : compute mask for a basin defined by a vector of points |
| 3 |
% |
| 4 |
%notes : - algorithm is based on great circles |
| 5 |
% - individual basin needs to be convex |
| 6 |
% |
| 7 |
%example of usage : |
| 8 |
%fld=read_bin('v4_basin.bin',0,1); |
| 9 |
%fld(fld>17)=18;%reset indices to proper basin |
| 10 |
%fld(fld>17&mygrid.XC>-40)=19;%reset indices to proper basin |
| 11 |
%msk1=v4_basin_one(1); msk2=v4_basin_one(2); |
| 12 |
%fld(fld<10&(msk1==1|msk2==1))=20;%add new basin index |
| 13 |
|
| 14 |
gcmfaces_global; |
| 15 |
|
| 16 |
if numBasin==1; |
| 17 |
lonList=[10 19 60 60 ]; |
| 18 |
latList=[60 80 81 60]; |
| 19 |
elseif numBasin==2; |
| 20 |
lonList=[60 100 110 60]; |
| 21 |
latList=[81 80 60 60]; |
| 22 |
else; |
| 23 |
error('unknown basin'); |
| 24 |
end; |
| 25 |
|
| 26 |
lonList=[lonList lonList(1)]; |
| 27 |
latList=[latList latList(1)]; |
| 28 |
|
| 29 |
msk=mygrid.mskC(:,:,1); |
| 30 |
for iy=1:length(lonList)-1; |
| 31 |
|
| 32 |
lonPair=lonList(iy:iy+1); |
| 33 |
latPair=latList(iy:iy+1); |
| 34 |
|
| 35 |
lon=mygrid.XC; lat=mygrid.YC; |
| 36 |
x=cos(lat*pi/180).*cos(lon*pi/180); |
| 37 |
y=cos(lat*pi/180).*sin(lon*pi/180); |
| 38 |
z=sin(lat*pi/180); |
| 39 |
|
| 40 |
x0=cos(latPair*pi/180).*cos(lonPair*pi/180); |
| 41 |
y0=cos(latPair*pi/180).*sin(lonPair*pi/180); |
| 42 |
z0=sin(latPair*pi/180); |
| 43 |
|
| 44 |
%get the rotation matrix: |
| 45 |
%1) rotate around x axis to put first point at z=0 |
| 46 |
theta=atan2(-z0(1),y0(1)); |
| 47 |
R1=[[1;0;0] [0;cos(theta);sin(theta)] [0;-sin(theta);cos(theta)]]; |
| 48 |
tmp0=[x0;y0;z0]; tmp1=R1*tmp0; x1=tmp1(1,:); y1=tmp1(2,:); z1=tmp1(3,:); |
| 49 |
x0=x1; y0=y1; z0=z1; |
| 50 |
%2) rotate around z axis to put first point at y=0 |
| 51 |
theta=atan2(x0(1),y0(1)); |
| 52 |
R2=[[cos(theta);sin(theta);0] [-sin(theta);cos(theta);0] [0;0;1]]; |
| 53 |
tmp0=[x0;y0;z0]; tmp1=R2*tmp0; x1=tmp1(1,:); y1=tmp1(2,:); z1=tmp1(3,:); |
| 54 |
x0=x1; y0=y1; z0=z1; |
| 55 |
%3) rotate around y axis to put second point at z=0 |
| 56 |
theta=atan2(-z0(2),-x0(2)); |
| 57 |
R3=[[cos(theta);0;-sin(theta)] [0;1;0] [sin(theta);0;cos(theta)]]; |
| 58 |
tmp0=[x0;y0;z0]; tmp1=R3*tmp0; x1=tmp1(1,:); y1=tmp1(2,:); z1=tmp1(3,:); |
| 59 |
x0=x1; y0=y1; z0=z1; |
| 60 |
|
| 61 |
%apply rotation to grid: |
| 62 |
tmpx=convert2array(x); tmpy=convert2array(y); tmpz=convert2array(z); |
| 63 |
tmp1=find(~isnan(tmpx)); |
| 64 |
tmpx2=tmpx(tmp1); tmpy2=tmpy(tmp1); tmpz2=tmpz(tmp1); |
| 65 |
tmp2=[tmpx2';tmpy2';tmpz2']; |
| 66 |
tmp3=R3*R2*R1*tmp2; |
| 67 |
tmpx2=tmp3(1,:); tmpy2=tmp3(2,:); tmpz2=tmp3(3,:); |
| 68 |
tmpx(tmp1)=tmpx2; tmpy(tmp1)=tmpy2; tmpz(tmp1)=tmpz2; |
| 69 |
x=convert2array(tmpx); y=convert2array(tmpy); z=convert2array(tmpz); |
| 70 |
|
| 71 |
%compute the great circle mask: |
| 72 |
mskCint=1*(z>0); |
| 73 |
|
| 74 |
%increment msk: |
| 75 |
msk(mskCint>0)=0; |
| 76 |
end; |
| 77 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph