source: trunk/SRC/Interpolation/fromreg.pro

;+
;
; @file_comments
; interpolate data from a "regular/rectangular grid" to any grid.
;
; 2 methods available: bilinear and imoms3
;
; A "regular/rectangular grid" is defined as a grid for which
;
; Each longitude lines have the same latitude and each latitude columns
; have the same longitude.
;
; @categories
; Interpolation
;
; @param method {in}{required}{type=string}
; the interpolation method.
; must be 'bilinear' or 'imoms3'
;
; @param datain {in}{required}{type=2d array}
; the input data to interpolate
;
; @param lonin {in}{required}{type=1d or 2d array}
; the longitude of the input data
;
; @param latin {in}{required}{type=1d or 2d array}
; the latitude of the input data
;
; @param lonout {in}{required}{type=1d or 2d array}
; the longitude of the output data
;
; @param latout {in}{required}{type=1d or 2d array}
; the latitude of the output data
;
; @keyword WEIG {type=2d array or variable name}
; (see ADDR)
;
; @keyword ADDR {type=2d array or variable name}
; 1) at the first call of fromreg:
; This keyword can be set to a named variable (undefined or equal to 0) into which the
; addresses used to perform the interpolation will be copied when the current routine exits.
; 2) Next, once this keyword is set to a defined 2d array, it is used to bypass the computation
; of the weights and addresses used to perform the interpolation. In this case, fromreg simply
; compute the interpolated field as:
;          dataout = total(weig*datain[addr], 1)
;          dataout = reform(dataout, jpio, jpjo, /over)
;
; In that case, method, lonin, latin, are not used (but are necessary).
; lonout, latout are used only to know the output domain size
;
; @keyword NONORTHERNLINE
; activate if you don't want to take into account the northern line
; of the input data when performing the interpolation.
;
; @keyword NOSOUTHERNLINE
; activate if you don't want to take into account the southern line
; of the input data when performing the interpolation.
;
; @keyword
; _EXTRA to be able to call fromreg with _extra keyword
;
; @returns
; 2D array the interpolated data
;
; @restrictions
; We supposed the data are located on a sphere, with a periodicity along the
; longitude.
;
; @examples
;
;  To interpolate 1 field:
;
;   IDL> topa = fromreg('bilinear', tncep, xncep, yncep, glamt, gphit)
;
;  or if you have several fields to interpolate from the same source and target grids
;
; 1) get back the weights and addresses in variables a and b
;   (that must be undefined or equal to 0 before calling fromreg
;
;   IDL> t1opa = fromreg('bilinear', t1ncep, xncep, yncep, glamt, gphit, WEIG = a, ADDR = b)
;   IDL> help, a, b
;
; 2) use a and b that are now defined to bypass the computation of the weights and addresses
; and speed-up the computation!
;
;   IDL> t2opa = fromreg('bilinear', t2ncep, xncep, yncep, glamt, gphit, WEIG = a, ADDR = b)
;
; @history
;  November 2005: Sebastien Masson (smasson\@lodyc.jussieu.fr)
;
; @version
; $Id$
;
;-
FUNCTION fromreg, method, datain, lonin, latin, lonout, latout $
                  , WEIG=weig, ADDR=addr $
                  , NONORTHERNLINE=nonorthernline $
                  , NOSOUTHERNLINE=nosouthernline, _EXTRA=ex
;
  compile_opt idl2, strictarrsubs
;
;---------------
; atmospheric grid parameters
;---------------
    alon = lonin
    alat = latin
    get_gridparams, alon, alat, jpia, jpja, 1, /double
;---------------
; Oceanic grid parameters
;---------------
    olon = lonout
    olat = latout
    get_gridparams, olon, olat, jpio, jpjo, 2, /double
;---------------
; Compute weight and address
;---------------
  IF NOT (keyword_set(weig) AND keyword_set(addr)) THEN BEGIN
    CASE method OF
      'bilinear':compute_fromreg_bilinear_weigaddr, alon, alat, olon, olat, weig, addr, NONORTHERNLINE = nonorthernline, NOSOUTHERNLINE = nosouthernline
      'imoms3':  compute_fromreg_imoms3_weigaddr,   alon, alat, olon, olat, weig, addr, NONORTHERNLINE = nonorthernline, NOSOUTHERNLINE = nosouthernline
      ELSE:BEGIN
        ras = report(' unknown interpolation method... we stop')
        stop
      ENDELSE
    ENDCASE
  ENDIF
;
;---------------
; do the interpolation
;---------------
  intype = size(datain, /type)
  if intype LE 3 THEN dataout = total(weig*float(datain[addr]), 1) $
  ELSE                dataout = total(weig*      datain[addr] , 1)
  dataout = reform(dataout, jpio, jpjo, /over)
  IF intype LE 3 THEN dataout = round(temporary(dataout))
;
  RETURN, dataout
END