1 | ;+ |
---|
2 | ; |
---|
3 | ; @hidden |
---|
4 | ; |
---|
5 | ; @file_comments |
---|
6 | ; recursive function. |
---|
7 | ; for one given point on an island, check if its neighbourgs are on the same island |
---|
8 | ; |
---|
9 | ; @param msk {in}{required}{type=2D array of 0 and 1} |
---|
10 | ; land-sea mask with 0 on the land and 1 on the ocean |
---|
11 | ; |
---|
12 | ; @param nx {in}{required}{type=scalar} |
---|
13 | ; x size of the mask |
---|
14 | ; |
---|
15 | ; @param ny {in}{required}{type=scalar} |
---|
16 | ; y size of the mask |
---|
17 | ; |
---|
18 | ; @param indin {in}{required}{type=1D array} |
---|
19 | ; index listing the point of the mask which are on the island |
---|
20 | ; |
---|
21 | ; @param numb {in}{required}{type=scalar} |
---|
22 | ; number of the island |
---|
23 | ; |
---|
24 | ;- |
---|
25 | ; |
---|
26 | PRO mskneig, msk, nx, ny, indin, numb |
---|
27 | ; |
---|
28 | compile_opt idl2, strictarrsubs |
---|
29 | ; |
---|
30 | ; flag the point |
---|
31 | msk[indin] = numb |
---|
32 | ; find its neighbourgs |
---|
33 | indx = (indin MOD nx) + [-1L, 0L, 1L] |
---|
34 | indx = (temporary(indx) + nx) MOD nx ; x periodicity |
---|
35 | indy = (indin/nx) + [-1L, 0L, 1L] |
---|
36 | indy = 0 > temporary(indy) < (ny-1L) ; y periodicity |
---|
37 | ; build 1d index |
---|
38 | ind = temporary(indx)#replicate(1L, 3) + replicate(nx, 3)#temporary(indy) |
---|
39 | ; for each neighbourg on the same island, call again mskneig |
---|
40 | FOR i = 0, n_elements(ind)-1 DO BEGIN |
---|
41 | IF msk[ind[i]] EQ 0 THEN mskneig, msk, nx, ny, ind[i], numb |
---|
42 | ENDFOR |
---|
43 | |
---|
44 | return |
---|
45 | END |
---|
46 | ; |
---|
47 | ;+ |
---|
48 | ; |
---|
49 | ; @file_comments |
---|
50 | ; given a 2D land-sea mask, give a number to each island |
---|
51 | ; |
---|
52 | ; @categories |
---|
53 | ; Grid |
---|
54 | ; |
---|
55 | ; @param mskin {in}{required}{type=2D array of 0 and 1} |
---|
56 | ; land-sea mask with 0 on the land and 1 on the ocean |
---|
57 | ; |
---|
58 | ; @returns |
---|
59 | ; A 2d array with 0 on the ocean and a different number for each island |
---|
60 | ; |
---|
61 | ; @restrictions |
---|
62 | ; - mask must be periodic along the x direction. |
---|
63 | ; - Periodicity along y direction is not taken into account |
---|
64 | ; - must have less than ~32700 islands; |
---|
65 | ; |
---|
66 | ; @examples |
---|
67 | ; IDL> a = numbisland(tmask[*,*,0]) |
---|
68 | ; |
---|
69 | ; @history |
---|
70 | ; Jan 2006: sebastien masson (smasson\@locean-ipsl.upmc.fr) |
---|
71 | ; |
---|
72 | ; @version |
---|
73 | ; $Id$ |
---|
74 | ; |
---|
75 | ;- |
---|
76 | ; |
---|
77 | FUNCTION numbisland, mskin |
---|
78 | ; |
---|
79 | compile_opt idl2, strictarrsubs |
---|
80 | ; |
---|
81 | time1 = systime(1) ; performance measurement |
---|
82 | szmsk = size(reform(mskin)) |
---|
83 | IF szmsk[0] NE 2 THEN stop |
---|
84 | nx = szmsk[1] |
---|
85 | ny = szmsk[2] |
---|
86 | msk = fix(mskin) |
---|
87 | |
---|
88 | islnumb = 10 ; default value |
---|
89 | |
---|
90 | land = (where(msk EQ 0, count))[0] |
---|
91 | WHILE count NE 0 DO BEGIN |
---|
92 | IF (islnumb-9) MOD 10 EQ 0 THEN BEGIN |
---|
93 | ras = report('island number :'+strtrim(islnumb-9, 1)) |
---|
94 | ENDIF |
---|
95 | mskneig, msk, nx, ny, land, islnumb |
---|
96 | |
---|
97 | land = (where(msk EQ 0, count))[0] |
---|
98 | islnumb = islnumb + 1 |
---|
99 | ENDWHILE |
---|
100 | |
---|
101 | msk = msk-9 |
---|
102 | msk[where(msk EQ -8)] = 0 |
---|
103 | |
---|
104 | print, 'time:', systime(1)-time1 |
---|
105 | |
---|
106 | RETURN, msk |
---|
107 | END |
---|