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