Changeset 383

Timestamp:

10/30/08 14:55:33 (16 years ago)

Author:

smasson

Message:

bugfix in file_interp with missing values

File:

• trunk/SRC/Interpolation/file_interp.pro (modified) (4 diffs)

trunk/SRC/Interpolation/file_interp.pro

@@ -15 +15 @@
 ;
 ; take care of NaN values
-  nanmsk = finite(data)
-  IF total(nanmsk) NE n_elements(nanmsk) THEN BEGIN
-    data[where(nanmsk EQ 0b)] = 1.e20 ; put large value to be sure they are removed during the interpolation
-    mask = inmask * temporary(nanmsk)
+  nanmask = finite(data)
+  totnanmask = total(nanmask)
+  IF totnanmask EQ 0 THEN return, !values.f_nan
+  IF totnanmask NE n_elements(nanmask) THEN BEGIN
+    data[where(nanmask EQ 0b)] = 1.e20 ; put large value to be sure they are removed during the interpolation
+    IF inmask[0] NE -1 THEN mask = temporary(nanmask) * inmask ELSE mask = temporary(nanmask)
   ENDIF ELSE mask = inmask
-;
+; take care of missing values
   tpmiss = size(missing_value, /type)
   IF tpmiss NE 0 AND tpmiss NE 7 THEN BEGIN
     CASE 1 OF
-      missing_value GT 1.e6:mask = temporary(mask) * (data LT (missing_value-10.))
-      missing_value LT -1.e6:mask = temporary(mask) * (data GT (missing_value+10.))
-      abs(missing_value) LT 1.e-6:mask = temporary(mask) * (abs(data) GT 1.e-6)
-      ELSE:mask = temporary(mask) * (data NE missing_value)
+      missing_value GT 1.e6:missmask = data LT (missing_value - 10.)
+      missing_value LT -1.e6:missmask = data GT (missing_value + 10.)
+      abs(missing_value) LT 1.e-6:missmask = abs(data) GT 1.e-6
+      ELSE:missmask = data NE missing_value
     ENDCASE
+    IF total(missmask) EQ 0 THEN return, missing_value
+    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)
+  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
@@ -436 +440 @@
                                                   , INIRR = inirr, METHOD = method, SMOOTH = smooth $
                                                   , WEIG = weig, ADDR = addr, MISSING_VALUE = var_missing_value, _extra = ex)
+              IF interp AND n_elements(data) EQ 1 THEN data = replicate(data, jpiout, jpjout)
               ncdf_varput, outid, outvarid[i], temporary(data)
             END
@@ -451 +456 @@
                                                     , INIRR = inirr, METHOD = method, SMOOTH = smooth $
                                                     , WEIG = weig, ADDR = addr, MISSING_VALUE = var_missing_value, _extra = ex)
+                IF interp AND n_elements(data) EQ 1 THEN data = replicate(data, jpiout, jpjout)
                 ncdf_varput, outid, outvarid[i], temporary(data), offset = off, count = outcnt
               ENDFOR
@@ -467 +473 @@
                                                       , INIRR = inirr, METHOD = method, SMOOTH = smooth $
                                                       , WEIG = weig, ADDR = addr, MISSING_VALUE = var_missing_value, _extra = ex)
+                  IF interp AND n_elements(data) EQ 1 THEN data = replicate(data, jpiout, jpjout)
                   ncdf_varput, outid, outvarid[i], temporary(data), offset = off, count = outcnt
                 ENDFOR