source: trunk/SRC/Grid/numbisland.pro @ 216

;+
; @hidden
;
; @file_comments
; recursive function.
; for one given point on an island, check if its neighbourgs are on the same island
;
; @param msk {in}{required}{type=2D array of 0 and 1}
; land-sea mask with 0 on the land and 1 on the ocean
;
; @param nx {in}{required}{type=scalar}
; x size of the mask
;
; @param ny {in}{required}{type=scalar}
; y size of the mask
;
; @param indin {in}{required}{type=1D array}
; index listing the point of the mask which are on the island
;
; @param numb {in}{required}{type=scalar}
; number of the island
;
;-
PRO mskneig, msk, nx, ny, indin, numb
; flag the point
  msk[indin] = numb
; find its neighbourgs
  indx = (indin MOD nx) + [-1L, 0L, 1L] 
  indx = (temporary(indx) + nx) MOD nx ; x periodicity
  indy = (indin/nx) + [-1L, 0L, 1L] 
  indy = 0 > temporary(indy) < (ny-1L) ; y periodicity
; build 1d index
  ind = temporary(indx)#replicate(1L, 3) + replicate(nx, 3)#temporary(indy)
; for each neighbourgs on the same island, call again mskneig
  FOR i = 0, n_elements(ind)-1 DO BEGIN  
    IF msk[ind[i]] EQ 0 THEN mskneig, msk, nx, ny, ind[i], numb
  ENDFOR 

  return
END

;+
; @file_comments
; given a 2D land-se mask, give a number to each islands
;
; @categories
; Grid
;
; @param mskin {in}{required}{type=2D array of 0 and 1}
; land-sea mask with 0 on the land and 1 on the ocean
;
; @returns
; A 2d array with 0 on the ocean and a different number for each island
;
; @restrictions
;    - mask must be periodic along the x direction.
;    - Periodicity along y direction is not taken into account
;    - must have less than ~32700 islands;
;
; @examples
;   IDL> a = numbisland(tmask[*,*,0])
;
; @history
;  Jan 2006: sebastien masson (smasson\@locean-ipsl.upmc.fr)
;
; @version
; $Id$
;-
;
FUNCTION numbisland, mskin
  time1 = systime(1) ; performance measurment
  szmsk = size(reform(mskin))
  IF szmsk[0] NE 2 THEN stop
  nx = szmsk[1]
  ny = szmsk[2]
  msk = fix(mskin)

  islnumb = 10 ; default value

  land = (where(msk EQ 0, count))[0]
  WHILE count NE 0 DO BEGIN
    IF (islnumb-9) MOD 10 EQ 0 THEN print,  'island number :'+strtrim(islnumb-9, 1)

    mskneig, msk, nx, ny, land, islnumb 
    
    land = (where(msk EQ 0, count))[0]
    islnumb = islnumb + 1
  ENDWHILE

  msk = msk-9
  msk[where(msk EQ -8)] = 0

  print, 'time:', systime(1)-time1 

  RETURN, msk
END