;+
;
; @file_comments
; similar to extrapolate but could to the job in a better way
; because the extrapolated values are smoothed...
; takes more time than extrapolate.
; extrapolate data where mskin is equal 0 by filling
; step by step the coastline points with the mean value of the 8 neighbors.
;
; @categories
; Interpolation
;
; @param in {in}{required}{type=2d array}
; data to be extrapolate
;
; @param mskin {in}{required}{type=2d array or -1}
; a 2D array, the land-sea mask of the output data (1 on ocean, 0 on land)
; put -1 if input data are not masked
;
; @keyword MINVAL {type=scalar}{default=not used}
; to specify a minimum value to the extrapolated values
;
; @keyword MAXVAL {type=scalar}{default=not used}
; to specify a maximum value to the extrapolated values
;
; @keyword GE0 {type=scalar 0 or 1}{default=0}
; put 1 to force the extrapolated values to be larger than 0, same as using minval=0.
;
; @keyword
; _EXTRA to be able to call extrapsmooth with _extra keyword
;
; @returns
; the extrapolated array with no more masked values
;
; @restrictions
; You cannot specify the number of iterations done in the extrapolation process
;
; @examples
; IDL> a=extrapsmooth(dist(jpi,jpj),tmask[*,*,0],/x_periodic)
; IDL> tvplus, a
; compare to extrapolate result:
; IDL> b=extrapolate(dist(jpi,jpj),tmask[*,*,0],/x_periodic)
; IDL> tvplus, b, window = 1
;
; @history
; January 2007: Sebastien Masson (smasson\@lodyc.jussieu.fr)
;
; @version
; $Id$
;-
FUNCTION extrapsmooth, in, mskin, X_PERIODIC = x_periodic, MINVAL = minval, MAXVAL = maxval, GE0 = ge0, _EXTRA = ex
;
compile_opt strictarr, strictarrsubs
;
sz = size(reform(in))
IF sz[0] NE 2 THEN BEGIN
ras = report('Input arrays must have 2 dimensions')
return, -1
ENDIF
nx = sz[1]
ny = sz[2]
IF n_elements(mskin) EQ 1 AND mskin[0] EQ -1 THEN mskin = replicate(1b, nx, ny)
IF n_elements(mskin) NE nx*ny THEN BEGIN
ras = report('input grid mask do not have the good size')
return, -1
ENDIF
;
out = reform(in)
whmsk = where(mskin EQ 0, nbr)
IF nbr NE 0 THEN out[temporary(whmsk)] = !values.f_nan
; add values on each side of the array to avoid boundary effects
nx2 = nx/2
ny2 = ny/2
add = replicate(!values.f_nan, nx, ny2)
out = [[add], [temporary(out)], [add]]
IF keyword_set(x_periodic) THEN BEGIN
add1 = out[0:nx2, *]
add2 = out[nx-nx2:nx-1, *]
out = [add2, [temporary(out)], add1]
ENDIF ELSE BEGIN
add = replicate(!values.f_nan, nx2, ny+2*ny2)
out = [add, [temporary(out)], add]
ENDELSE
;
msk0 = where(finite(out) EQ 0)
nnan = total(finite(out, /nan))
i = 1
WHILE nnan NE 0 DO BEGIN
tmp = smooth(out, 2*i + 1, /nan)
; find only the changed values that where on land
new = where(finite(out) EQ 0 AND finite(tmp) EQ 1, nnew)
IF nnew EQ 0 then nnan = 0 ELSE BEGIN
IF keyword_set(ge0) THEN tmp = 0. > temporary(tmp)
IF n_elements(minval) NE 0 THEN tmp = minval > temporary(tmp)
IF n_elements(maxval) NE 0 THEN tmp = temporary(tmp) < maxval
out[new] = tmp[new]
i = i+1
nnan = total(finite(out, /nan))
ENDELSE
ENDWHILE
; a final smooth is needed
out = smooth(temporary(out), 5, /nan)
IF keyword_set(ge0) THEN out = 0. > temporary(out)
IF n_elements(minval) NE 0 THEN out = minval > temporary(out)
IF n_elements(maxval) NE 0 THEN out = temporary(out) < maxval
; get back to the original domain
out = (temporary(out))[nx2:nx+nx2-1, ny2:ny+ny2-1]
; put back the non-maskqed values
whmsk = where(mskin NE 0)
out[whmsk] = in[whmsk]
;
return, out
END