diff --git a/src/clib/xtg/surf_cube_attr_intv.c b/src/clib/xtg/surf_cube_attr_intv.c index 25eb47157..72bb2d698 100644 --- a/src/clib/xtg/surf_cube_attr_intv.c +++ b/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; diff --git a/src/clib/xtg/surf_stack_slice_cube.c b/src/clib/xtg/surf_stack_slice_cube.c index ae8425968..25b423dc4 100644 --- a/src/clib/xtg/surf_stack_slice_cube.c +++ b/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); diff --git a/src/xtgeo/surface/regular_surface.py b/src/xtgeo/surface/regular_surface.py index 2bf0b774f..3753088d6 100644 --- a/src/xtgeo/surface/regular_surface.py +++ b/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) + )