Changeset 433

Timestamp:

06/14/10 12:31:14 (14 years ago)

Author:

smasson

Message:

add box mean interpolation in file_interp

Location:

trunk/SRC/Interpolation

Files:

• file_interp.pro (modified) (14 diffs)
• fromirr.pro (modified) (1 diff)
• fromreg.pro (modified) (1 diff)

trunk/SRC/Interpolation/file_interp.pro

@@ +1 @@
+;+
+;
+; @hidden
+;
+;-
+FUNCTION boxmean_interp, data, divx, divy
+ ;
+  compile_opt idl2, strictarrsubs
+;
+  sz = size(data, /dimensions )
+  jpiout = sz[0] / divx
+  jpjout = sz[1] / divy
+;
+  data = reform(data, divx, jpiout, divy, jpjout, /overwrite)
+  data = total(temporary(data), 1)   ; ave along divx
+  data = total(temporary(data), 2)   ; ave along divy
+  data = temporary(data) / float(divx*divy)
+
+  RETURN, data
+END
 ;+
 ;
@@ -7 +27 @@
                         , INIRR = inirr, METHOD = method, SMOOTH = smooth $
                         , WEIG = weig, ADDR = addr, MISSING_VALUE = missing_value $
+                        , DIVX = divx, DIVY = divy $
                         , GETHAN = gethan, LETHAN = lethan, _EXTRA = ex
 ;
@@ -35 +56 @@
     IF mask[0] NE -1 THEN mask = temporary(missmask) * mask ELSE mask = temporary(missmask)
   ENDIF
+
 ; extrapolation
   IF keyword_set(smooth) THEN data = extrapsmooth(temporary(data), mask, /x_periodic, _extra = ex) $
   ELSE data = extrapolate(temporary(data), mask, /x_periodic, _extra = ex)
 ; interpolation
-  IF NOT keyword_set(inirr) THEN BEGIN
-    data = fromreg(method, temporary(data), inlon, inlat, outlon, outlat, WEIG = weig, ADDR = addr, _extra = ex)
-  ENDIF ELSE BEGIN
-    data = fromirr(method, temporary(data), inlon, inlat, -1, outlon, outlat, -1, WEIG = weig, ADDR = addr)
-  ENDELSE
+  IF method EQ 'boxmean' THEN BEGIN
+    data = boxmean_interp(temporary(data), divx, divy)
+  ENDIF ELSE BEGIN 
+    IF NOT keyword_set(inirr) THEN BEGIN
+      data = fromreg(method, temporary(data), inlon, inlat, outlon, outlat, WEIG = weig, ADDR = addr, _extra = ex)
+    ENDIF ELSE BEGIN
+      data = fromirr(method, temporary(data), inlon, inlat, -1, outlon, outlat, -1, WEIG = weig, ADDR = addr)
+    ENDELSE
+  ENDELSE 
 
   IF n_elements(gethan) EQ 1 THEN data = gethan > temporary(data)
@@ -66 +92 @@
 ; output file name (will be overwritten if already exist)
 ;
-; @param gridout {in}{type=scalar string}
-; output grid file name (must exist and must contain the
-; longitude and latitude axis as 1D or 2D arrays)
+; @param gridout {in}{type=scalar string or 2 element vector if boxmean method}
+; if boxmean method:
+;  2 elements vector defining the size of the box used to compute the
+;  mean value. It must divide the original grid size.
+; else:
+;  output grid file name (must exist and must contain the
+;  longitude and latitude axis as 1D or 2D arrays)
 ;
 ; @keyword GRIDIN {type=scalar string}{default=set to filein}
@@ -96 +126 @@
 ;
 ; @keyword METHOD {type=scalar string}{default='bilinear'}
-; interpolation method: can be only 'bilinear' (or 'imoms3' if the input grid
+; interpolation method: can be only 'bilinear', 'boxmean' (or 'imoms3' if the input grid
 ; is a "regular" grid). 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.
+; boxmean interpolation is simply the mean value over a box of nx by ny
+; points (see gridout definition).
 ;
 ; @keyword SMOOTH {type=scalar 0 or 1}{default=0}
@@ -248 +280 @@
 ;   IDL> file_interp, filein, fileout, gridout, inxaxisname = 'lo', inyaxisname = 'la', keep = ['lo', 'la', 'cond_sed']
 ;   IDL> file_interp, in, out, gdout, inuseasmask = 'sst', inmissing_value = -1.00000e+30, missing_value = -1000.00
+;   IDL> file_interp,'sst_reg025.nc', 'sst_reg1.nc',[4,4], method = 'boxmean'
 ;
 ; @history
@@ -308 +341 @@
   IF keyword_set(set_yaxisname) THEN inyaxisname = set_yaxisname
 ;
-  ncdf_getaxis, gridout, outdimidx, outdimidy, outlon, outlat, xaxisname = outxaxisname, yaxisname = outyaxisname
-  get_gridparams, outlon, outlat, jpiout, jpjout, 2
+  IF method EQ 'boxmean' THEN BEGIN
+    IF n_elements(gridout) NE 2 THEN stop
+    divx = round(gridout[0])
+    divy = round(gridout[1])
+    IF jpiin MOD divx NE 0 THEN BEGIN 
+      print, 'in boxmean method, the x size ('+strtrim(divx, 1)+') of the box used to average the data must devide the size of the x dimension ('+strtrim(jpiin, 1)+')'
+      return
+    ENDIF
+    IF jpjin MOD divy NE 0 THEN BEGIN 
+      print, 'in boxmean method, the y size ('+strtrim(divy, 1)+') of the box used to average the data must devide the size of the y dimension ('+strtrim(jpjin, 1)+')'
+      return
+    ENDIF
+    jpiout = jpiin / divx
+    jpjout = jpjin / divy
+    outlon = inlon   &   outlon = boxmean_interp(outlon, divx, divy)
+    outlat = inlat   &   outlat = boxmean_interp(outlat, divx, divy)
+    IF jpiout EQ 1 OR jpjout EQ 1 THEN BEGIN   
+      outlon = reform(outlon, jpiout, jpjout, /overwrite)
+      outlat = reform(outlat, jpiout, jpjout, /overwrite)
+    ENDIF 
+  ENDIF ELSE BEGIN 
+    ncdf_getaxis, gridout, outdimidx, outdimidy, outlon, outlat, xaxisname = outxaxisname, yaxisname = outyaxisname
+    get_gridparams, outlon, outlat, jpiout, jpjout, 2
+  ENDELSE 
 ;
 ; masks
@@ -318 +373 @@
   IF size(inmask, /n_dimensions) EQ 2 THEN inmasksz = [inmasksz, 0]
   IF n_elements(inmaskname) EQ 0 THEN inmaskname = 'not defined' ; default definition
-  outmask = ncdf_getmask(maskout, MASKNAME = outmaskname, INVMASK = outinvmask, USEASMASK = outuseasmask $
-                        , MISSING_VALUE = outmissing_value, ADDSCL_BEFORE = outaddscl_before, TESTOP = outtestop)
+  IF method EQ 'boxmean' THEN BEGIN
+    outmask = inmask
+    IF inmask[0] NE -1 THEN outmask = boxmean_interp(outmask, divx, divy)
+  ENDIF ELSE BEGIN 
+    outmask = ncdf_getmask(maskout, MASKNAME = outmaskname, INVMASK = outinvmask, USEASMASK = outuseasmask $
+                           , MISSING_VALUE = outmissing_value, ADDSCL_BEFORE = outaddscl_before, TESTOP = outtestop)
+  ENDELSE 
 ;
 ; irregular grids?
@@ -338 +398 @@
   ENDCASE
 
-  IF inirr AND method NE 'bilinear' THEN stop
+  IF inirr AND method EQ 'imoms3' THEN stop
 ;
 ; Dimensions
@@ -451 +511 @@
                                                   , INIRR = inirr, METHOD = method, SMOOTH = smooth $
                                                   , WEIG = weig, ADDR = addr, MISSING_VALUE = var_missing_value $
+                                                  , DIVX = divx, DIVY = divy $
                                                   , GETHAN = gethan, LETHAN = lethan, _extra = ex)
               IF interp AND n_elements(data) EQ 1 THEN data = replicate(data, jpiout, jpjout)
@@ -457 +518 @@
             3:BEGIN             ; 3D
               FOR k = 0, indimsz[varinq.dim[2]]-1 DO BEGIN
+                IF k MOD 100 EQ 0 THEN print, k
                 incnt = [indimsz[varinq.dim[0: 1]], 1]
                 outcnt = [outdimsz[varinq.dim[0: 1]], 1]
@@ -468 +530 @@
                                                     , INIRR = inirr, METHOD = method, SMOOTH = smooth $
                                                     , WEIG = weig, ADDR = addr, MISSING_VALUE = var_missing_value $
+                                                    , DIVX = divx, DIVY = divy $
                                                     , GETHAN = gethan, LETHAN = lethan, _extra = ex)
                 IF interp AND n_elements(data) EQ 1 THEN data = replicate(data, jpiout, jpjout)
@@ -475 +538 @@
             4:BEGIN             ; 4D
               FOR t = 0, indimsz[varinq.dim[3]]-1 DO BEGIN
+                IF t MOD 100 EQ 0 THEN print, t
                 FOR k = 0, indimsz[varinq.dim[2]]-1 DO BEGIN
                   incnt = [indimsz[varinq.dim[0: 1]], 1, 1]
@@ -486 +550 @@
                                                       , INIRR = inirr, METHOD = method, SMOOTH = smooth $
                                                       , WEIG = weig, ADDR = addr, MISSING_VALUE = var_missing_value $
+                                                      , DIVX = divx, DIVY = divy $
                                                       , GETHAN = gethan, LETHAN = lethan, _extra = ex)
                   IF interp AND n_elements(data) EQ 1 THEN data = replicate(data, jpiout, jpjout)

trunk/SRC/Interpolation/fromirr.pro

@@ -121 +121 @@
   ENDIF
 ;---------------
-; to the interpolation
+; do the interpolation
 ;---------------
-  dataout = total(weig*datain[addr], 1)
+  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

trunk/SRC/Interpolation/fromreg.pro

@@ -126 +126 @@
   ENDIF
 ;
-  dataout = total(weig*datain[addr], 1)
+;---------------
+; 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