source: trunk/Interpolation/extrapolate.pro @ 69

FUNCTION extrapolate, zinput, maskinput, nb_iteration, x_periodic = x_periodic, MINVAL = minval, MAXVAL = maxval
;
  compile_opt strictarr, strictarrsubs 
;
; extrapolate data (zinput) where maskinput eq 0 by filling step by
; step the coastline points with the mean value of the 8 neighbourgs.
;
; check the number of iteration used in the extrapolation.
  IF n_elements(nb_iteration) EQ 0 THEN nb_iteration = 10.E20
  IF nb_iteration EQ 0 THEN return, zinput
  nx = (size(zinput))[1]
  ny = (size(zinput))[2]
; take care of the boundary conditions...
; 
; for the x direction, we put 2 additional columns at the left and
; right side of the array. 
; for the y direction, we put 2 additional lines at the bottom and
; top side of the array. 
; These changes allow us to use shift function without taking care of
; the x and y periodicity.
;
  ztmp = bytarr(nx+2, ny+2)
  ztmp[1:nx, 1:ny] = byte(maskinput)
  msk = temporary(ztmp)
;
  ztmp = replicate(1.e20, nx+2, ny+2)
  ztmp[1:nx, 1:ny] = zinput
  if keyword_set(x_periodic) then begin
    ztmp[0, 1:ny] = zinput[nx-1, *]
    ztmp[nx+1, 1:ny] = zinput[0, *]
  ENDIF
; remove NaN points if there is some...
  nan = where(finite(ztmp) EQ 0, cnt_nan)
  IF cnt_nan NE 0 THEN ztmp[temporary(nan)] = 1.e20
  z = temporary(ztmp)
  nx2 = nx+2
  ny2 = ny+2
;---------------------------------------------------------------
;---------------------------------------------------------------
; extrapolation 
;---------------------------------------------------------------
  sqrtinv = 1./sqrt(2)
;
  cnt = 1
; When we look for the coast line, we don't whant to select the
; borderlines of the array. -> we force the value of the mask for
; those lines.
  msk[0, *] = 1b
  msk[nx+1, *] = 1b
  msk[*, 0] = 1b
  msk[*, ny+1] = 1b
; find the land points
  land = where(msk EQ 0, cnt_land)
;---------------------------------------------------------------
  WHILE cnt LE nb_iteration AND cnt_land NE 0 DO BEGIN
;---------------------------------------------------------------
; find the coast line points...
;---------------------------------------------------------------
; Once the land points list has been found, we change back the the
; mask values for the boundary conditions.
    msk[0, *] = 0b
    msk[nx+1, *] = 0b
    msk[*, 0] = 0b
    msk[*, ny+1] = 0b
    if keyword_set(x_periodic) then begin
      msk[0, *] = msk[nx, *]
      msk[nx+1, *] = msk[1, *]
    endif
;
; we compute the weighted number of sea neighbourgs.
; those 4 neighbours have a weight of 1:
;    *    
;   *+*        
;    *     
;
; those 4 neighbours have a weight of 1/sqrt(2):
;
;    * *
;     +       
;    * *
;
; As we make sure that none of the land points are located on the
; border of the array, we can compute the weight without shift
; (faster).
;
    weight = msk[land+1]+msk[land-1]+msk[land+nx2]+msk[land-nx2] $
             +sqrtinv*(msk[land+nx2+1]+msk[land+nx2-1] $
                       +msk[land-nx2+1]+msk[land-nx2-1])
; list all the points that have sea neighbourgs
    ok = where(weight GT 0)
; the coastline points
    coast = land[ok]
; their weighted number of sea neighbourgs.
    weight = weight[temporary(ok)]
;---------------------------------------------------------------
; fill the coastine points
;---------------------------------------------------------------
    z = temporary(z)*msk
;
    zcoast = z[1+coast]+z[-1+coast]+z[nx2+coast]+z[-nx2+coast] $
             +1./sqrt(2)*(z[nx2+1+coast]+z[nx2-1+coast] $
                          +z[-nx2+1+coast]+z[-nx2-1+coast])
;    
    IF n_elements(minval) NE 0 THEN zcoast = minval > temporary(zcoast)
    IF n_elements(maxval) NE 0 THEN zcoast = temporary(zcoast) < maxval
    z[coast] = temporary(zcoast)/temporary(weight)
; we update the the boundary conditions of z
    if keyword_set(x_periodic) then begin
      z[0, *] = z[nx, *]
      z[nx+1, *] = z[1, *]
    endif
;---------------------------------------------------------------
; we update the mask
;---------------------------------------------------------------
    msk[temporary(coast)] = 1
;
    cnt = cnt + 1
; When we look for the coast line, we don't whant to select the
; borderlines of the array. -> we force the value of the mask for
; those lines.
    msk[0, *] = 1b
    msk[nx+1, *] = 1b
    msk[*, 0] = 1b
    msk[*, ny+1] = 1b
; find the land points
    land = where(msk EQ 0, cnt_land)
;
  ENDWHILE
;---------------------------------------------------------------
; we return the original size of the array
;---------------------------------------------------------------
  
  return, z[1:nx, 1:ny]
END