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

;+
;
; @hidden
;
; @file_comments
; recursive function.
; for one given point on an island, check if its neighbors 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
;
  compile_opt idl2, strictarrsubs
;
; flag the point
  msk[indin] = numb
; find its neighbors
  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 neighbor 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-sea mask, give a number to each island
;
; @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
;
  compile_opt idl2, strictarrsubs
;
  time1 = systime(1) ; performance measurement
  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 BEGIN
       ras = report('island number :'+strtrim(islnumb-9, 1))
    ENDIF
    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