Fix bugs #436 and #439

src/clib/xtg/surf_cube_attr_intv.c

@@ -276,8 +276,9 @@ surf_cube_attr_intv(int ncol,
     int i, ic;
     for (i = 0; i < nsurf1; i++) {
         for (ic = 0; ic <= ndiv; ic++) {
-            if (maskv1[i] == 0 && maskv2[i] == 0) {
-                stack[i][ic] = surfsv1[i] + ic * slicezinc;
+            double vadd = surfsv1[i] + ic * slicezinc;
+            if (maskv1[i] == 0 && maskv2[i] == 0 && vadd <= surfsv2[i]) {
+                stack[i][ic] = vadd;
                 rmask[i][ic] = 0;
             } else {
                 stack[i][ic] = UNDEF;

src/clib/xtg/surf_stack_slice_cube.c

@@ -5,7 +5,7 @@
  *    surf_stack_slice_cube.c
  *
  * DESCRIPTION:
- *    Slicing a stack of surfaces and reurn cube values
+ *    Slicing a stack of surfaces and return cube values
  *
  *    In this version, the surface and cube share geometry, which makes calculations
  *    way much simpler and probably also much faster
@@ -14,10 +14,8 @@
  *    ncol, nrow...  i     Cube dimensions and relevant increments
  *    czori, czinc   i     Cube settings
  *    cubevalsv      i     Cube values
- *    stack          i     map array stack with Z values
- *    rmask          i     Length of slice array
- *    surfsv        i/o    map to update
- *    maskv         i/o    mask array for map, may be updated!
+ *    stack         i/o    stacked map array stack with Z values
+ *    rmask          i     stacked map array stack with masks
  *    optnearest     i     If 1 use nerest node, else do interpolation aka trilinear
  *    optmask        i     If 1 then masked cells (undef) are made if cube is UNDEF
  *                         (surface outside cube will not occur in this case)
@@ -68,6 +66,8 @@ surf_stack_slice_cube(int ncol,
                     continue;
 
                 double zval = stack[icmap][nd];
+                if (zval > UNDEF_LIMIT)
+                    continue;
 
                 // find vertical index of node right above
                 int k1 = (int)((zval - czori) / czinc);

src/xtgeo/surface/regular_surface.py

@@ -43,7 +43,7 @@
 import math
 from types import FunctionType
 import warnings
-
+import numbers
 from collections import OrderedDict
 
 import numpy as np
@@ -481,6 +481,11 @@ def ny(self):  # pylint: disable=C0103
         warnings.warn("Deprecated; use nrow instead", DeprecationWarning)
         return self._nrow
 
+    @property
+    def dimensions(self):
+        """2-tuple: The surface dimensions as a tuple of 2 integers (read only)"""
+        return (self._ncol, self._nrow)
+
     @property
     def nactive(self):
         """Number of active map nodes (read only)."""
@@ -2552,13 +2557,13 @@ def translate_coordinates(self, translate=(0, 0, 0)):
     # Private
     # ==================================================================================
 
-    def _ensure_correct_values(self, values):
+    def _ensure_correct_values(
+        self, values
+    ):  # pylint: disable=too-many-branches, too-many-statements
         """Ensures that values is a 2D masked numpy (ncol, nrow), C order.
 
-        This is an improved but private version over ensure_correct_values
-
         Args:
-            values (array or scalar): Values to process.
+            values (array-like or scalar): Values to process.
 
         Return:
             Nothing, self._values will be updated inplace
@@ -2567,7 +2572,7 @@ def _ensure_correct_values(self, values):
         if not self._isloaded:
             return
 
-        if values is None:
+        if values is None or values is False:
             self._values = None
             return
 
@@ -2580,7 +2585,7 @@ def _ensure_correct_values(self, values):
 
             if (
                 currentmask is not None
-                and currentmask.all() == newmask.all()
+                and np.array_equal(currentmask, newmask)
                 and self._values.shape == values.shape
                 and values.flags.c_contiguous is True
             ):
@@ -2591,7 +2596,7 @@ def _ensure_correct_values(self, values):
                 self._values += vals
 
             else:
-                # replace any undef or nan with mask, return a view (copy=False)
+                # replace any undef or nan with mask
                 vals = values.astype(np.float64)
                 vals = ma.masked_greater(vals, self.undef_limit, copy=True)
                 vals = ma.masked_invalid(vals, copy=True)
@@ -2603,13 +2608,20 @@ def _ensure_correct_values(self, values):
                     vals = np.asanyarray(vals, order="C")
                     vals = ma.array(vals, mask=mask, order="C")
                 self._values = vals
-                return
 
-        elif np.isscalar(values):
+        elif isinstance(values, numbers.Number):
             if currentmask is not None:
                 self._values *= 0
-                self._values += values
-                return
+                vals = np.ones(self.dimensions, dtype=np.float64) * values
+                vals = np.ma.array(vals, mask=currentmask)
+
+                # there maybe cases where values scalar input is some kind of UNDEF
+                # which will change the mask
+                vals = ma.masked_greater(vals, self.undef_limit, copy=False)
+                vals = ma.masked_invalid(vals, copy=False)
+
+                self._values += vals
+
             else:
                 vals = ma.zeros((self.ncol, self.nrow), order="C", dtype=np.float64)
                 self._values = vals + float(values)
@@ -2623,7 +2635,35 @@ def _ensure_correct_values(self, values):
                 try:
                     vals = ma.reshape(vals, (self.ncol, self.nrow), order="C")
                 except ValueError as emsg:
-                    xtg.error("Cannot reshape array: {}".format(emsg))
+                    logger.critical("Cannot reshape array: %s", emsg)
                     raise
 
             self._values = vals
+
+        elif isinstance(values, (list, tuple)):  # ie values ~ list-like
+            vals = ma.array(values, order="C", dtype=np.float64)
+            vals = ma.masked_greater(vals, self.undef_limit, copy=True)
+            vals = ma.masked_invalid(vals, copy=True)
+
+            if vals.shape != (self.ncol, self.nrow):
+                try:
+                    vals = ma.reshape(vals, (self.ncol, self.nrow), order="C")
+                except ValueError as emsg:
+                    logger.critical("Cannot reshape array: %s", emsg)
+                    raise
+
+            self._values = vals
+
+        else:
+            raise ValueError("Input values are in invalid format: {}".format(values))
+
+        if self._values is not None:
+
+            if self._ilines is None:
+                self._ilines = np.array(range(1, self._ncol + 1), dtype=np.int32)
+                self._xlines = np.array(range(1, self._nrow + 1), dtype=np.int32)
+
+            if self._values.mask is ma.nomask:
+                self._values = ma.array(
+                    self._values, mask=ma.getmaskarray(self._values)
+                )