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