| 1 |
gforget |
1.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=v4_read_bin('/net/nares/raid11/ecco-shared/ecco-version-4/input/input_rest/basin_masks_ext_ecco_v4.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