WIP

src/xtgeo/cube/cube1.py

@@ -6,6 +6,7 @@
 import pathlib
 import tempfile
 import warnings
+from typing import Any, Optional
 
 import deprecation
 import numpy as np
@@ -130,7 +131,7 @@ def wrapper(cls, *args, **kwargs):
     return wrapper
 
 
-class Cube:  # pylint: disable=too-many-public-methods
+class Cube:  # pylint: disable=too-many-public-methods, W0201
     """Class for a (seismic) cube in the XTGeo framework.
 
     The values are stored as a 3D numpy array (4 bytes; float32 is default),
@@ -948,7 +949,14 @@ def from_roxar(self, project, name, folder=None):  # pragma: no cover
         self._metadata.required = self
 
     def to_roxar(
-        self, project, name, folder=None, domain="time", compression=("wavelet", 5)
+        self,
+        project: Any,
+        name: str,
+        folder: Optional[str] = None,
+        propname: str = "seismic_attribute",
+        domain: str = "time",
+        compression: tuple = ("wavelet", 5),
+        target: str = "seismic",
     ):  # pragma: no cover
         """Export (transfer) a cube from a XTGeo cube object to Roxar data.
 
@@ -958,30 +966,49 @@ def to_roxar(
             will not be saved until the user do an explicit project save action.
 
         Args:
-            project (str or roxar._project): Inside RMS use the magic 'project',
+            project: Inside RMS use the magic 'project',
                 else use path to RMS project, or a project reference
-            name (str): Name of cube (seismic data) within RMS project.
-            folder (str): Cubes may be stored under a folder in the tree, use '/'
+            name: Name of cube (seismic data) within RMS project.
+            folder: Cubes may be stored under a folder in the tree, use '/'
                 to seperate subfolders.
-            domain (str): 'time' (default) or 'depth'
+            propname: Name of grid property; only relevant when target is "grid" and
+                defaults to "seismic_attribute"
+            domain: 'time' (default) or 'depth'
             compression (tuple): description to come...
+            target: Optionally, the seismic cube can be written to the `Grid model`
+                tree in RMS. Internally, it will be convert to a "box" grid with one
+                gridproperty, before it is written to RMS. The ``compression``and
+                ``domain`` are not relevant when writing to grid model.
 
         Raises:
             To be described...
 
         Example::
 
-            zz = xtgeo.cube.Cube('myfile.segy')
-            zz.to_roxar(project, 'reek_cube')
-
-            # alternative
             zz = xtgeo.cube_from_file('myfile.segy')
             zz.to_roxar(project, 'reek_cube')
+            # write cube to "Grid model" tree in RMS instead
+            zz.to_roxar(project, 'cube_as_grid', propname="impedance", target="grid")
 
+        .. versionchanged:: 3.4 Add ``target`` and ``propname`` keys
         """
-        _cube_roxapi.export_cube_roxapi(
-            self, project, name, folder=folder, domain=domain, compression=compression
-        )
+
+        if "grid" in target.lower():
+            _tmpgrd = xtgeo.grid_from_cube(self, propname=name)
+            _tmpprop = _tmpgrd.props[0]
+            _tmpprop.name = propname if propname else "seismic_attribute"
+            _tmpgrd.to_roxar(project, name)
+            _tmpprop.to_roxar(project, name, _tmpprop.name)
+
+        else:
+            _cube_roxapi.export_cube_roxapi(
+                self,
+                project,
+                name,
+                folder=folder,
+                domain=domain,
+                compression=compression,
+            )
 
     @staticmethod
     def scan_segy_traces(sfile, outfile=None):
@@ -1054,15 +1081,15 @@ def _ensure_correct_values(self, values):
 
         if isinstance(values, numbers.Number):
             self._values = np.zeros(self.dimensions, dtype=np.float32) + values
-            return
 
-        if isinstance(values, np.ndarray):
+        elif isinstance(values, np.ndarray):
             values = values.reshape(self.dimensions)
 
             if not values.data.c_contiguous:
                 values = np.ascontiguousarray(values)
+            self._values = values
 
-        if isinstance(values, (list, tuple)):
-            values = np.array(values, dtype=np.float32).reshape(self.dimensions)
+        elif isinstance(values, (list, tuple)):
+            self._values = np.array(values, dtype=np.float32).reshape(self.dimensions)
 
-        self._values = values
+        self._values = self._values.astype(np.float32)

src/xtgeo/grid3d/grid.py

@@ -16,9 +16,19 @@
 import xtgeo
 from xtgeo.common import XTGDescription, _XTGeoFile
 
-from . import (_grid3d_fence, _grid_etc1, _grid_export, _grid_hybrid,
-               _grid_import, _grid_import_ecl, _grid_import_xtgcpgeom,
-               _grid_refine, _grid_roxapi, _grid_wellzone, _gridprop_lowlevel)
+from . import (
+    _grid3d_fence,
+    _grid_etc1,
+    _grid_export,
+    _grid_hybrid,
+    _grid_import,
+    _grid_import_ecl,
+    _grid_import_xtgcpgeom,
+    _grid_refine,
+    _grid_roxapi,
+    _grid_wellzone,
+    _gridprop_lowlevel,
+)
 from ._ecl_grid import Units
 from ._grid3d import _Grid3D
 from .grid_properties import GridProperties
@@ -176,24 +186,24 @@ def create_box_grid(
 
 
 def grid_from_cube(
-    cube: xtgeo.Cube, prop_name: str = "seismics", oricenter: bool = True
+    cube: xtgeo.Cube, propname: str = "seismics", oricenter: bool = True
 ):
     """Create a rectangular 'shoebox' grid from an existing cube.
 
-    The cube values itself will then be stored with name given by property key.
+    The cube values itself will then be stored with name given by ``propname`` key.
 
     Since the cube actually is node centered, while grids are cell oriented,
-    the geometries here are shifted half an increment as default. To avoiud this, use
+    the geometries here are shifted half an increment as default. To avoid this, use
     oricenter=False.
 
     Args:
         cube: The xtgeo Cube instance
-        prop_name: Name of seismic property, if None then only the grid geometry
+        propname: Name of seismic property, if None then only the grid geometry
             will be made
         oricenter: Default is True, to treat seismic nodes as cell center values in
             a grid.
 
-    .. versionadded:: 3.24
+    .. versionadded:: 3.4
     """
 
     grd = create_box_grid(
@@ -204,14 +214,15 @@ def grid_from_cube(
         rotation=cube.rotation,
         flip=cube.yflip,
     )
-    gprop = xtgeo.GridProperty(
-        ncol=cube.ncol,
-        nrow=cube.nrow,
-        nlay=cube.nlay,
-        values=cube.values.copy(),
-        name=prop_name,
-    )
-    grd.props = [gprop]
+    if propname is not None:
+        gprop = xtgeo.GridProperty(
+            ncol=cube.ncol,
+            nrow=cube.nrow,
+            nlay=cube.nlay,
+            values=cube.values.copy(),
+            name=propname,
+        )
+        grd.props = [gprop]
     return grd
 
 

src/xtgeo/surface/_regsurf_cube_window_v3.py

@@ -17,6 +17,7 @@
 
 warnings.filterwarnings(action="ignore", message="All-NaN slice encountered")
 warnings.filterwarnings(action="ignore", message="Mean of empty slice")
+warnings.filterwarnings(action="ignore", message="Degree of freedom")
 
 xtg = XTGeoDialog()
 
@@ -76,19 +77,17 @@ def _get_iso_maskthreshold_surface(
 
 
 def _proxy_surf_outside_cube(
-    _regsurf_: xtgeo.RegularSurface, cube: xtgeo.Cube, deadtraces: bool
+    _regsurf_: xtgeo.RegularSurface,
+    cube: xtgeo.Cube,
 ) -> xtgeo.RegularSurface:
     """Proxy for the part of input surface that is outside the cube area."""
     logger.info("Get a proxy for part of original surface being outside the cube")
     outside = _regsurf_.copy()
-    dummy_cube = cube.copy()
-    dummy_cube.values = xtgeo.UNDEF - 1  # to make cube values very large, but not UNDEF
-
-    zsurf = _regsurf_.copy()
-    zsurf.values = cube.zori  # to ensure that surface hits cube
+    outside.values = 0.0
 
-    outside.slice_cube(zsurf=zsurf, mask=False, deadtraces=deadtraces)
-    outside.values = np.ma.where(outside.values < (xtgeo.UNDEF - 1), 1, 0)
+    tmp_surf = xtgeo.surface_from_cube(cube, 0.0)
+    boundary = tmp_surf.get_boundary_polygons()
+    outside.set_outside(boundary, 1.0)
     return outside  # is 1 outside the cube area, 0 within the cube
 
 
@@ -226,6 +225,7 @@ def _compute_stats(
     upper: xtgeo.RegularSurface,
     masksurf: xtgeo.RegularSurface,
     sampling: str,
+    snapxy: bool,
 ) -> xtgeo.RegularSurface:
     """Compute the attribute and return the attribute map"""
     logger.info("Compute stats...")
@@ -287,10 +287,13 @@ def _compute_stats(
     sampled.values = np.ma.masked_invalid(values)
     sampled.values = np.ma.masked_where(masksurf.values == 0, sampled.values)
 
-    sampling_option = "nearest" if sampling in ("nearest", "cube") else "bilinear"
-    actual.resample(
-        sampled, sampling=sampling_option
-    )  # will be on input map resolution
+    if not snapxy:
+        sampling_option = "nearest" if sampling in ("nearest", "cube") else "bilinear"
+        actual.resample(
+            sampled, sampling=sampling_option
+        )  # will be on input map resolution
+    else:
+        actual = sampled
 
     logger.info("Compute stats... done")
     return actual
@@ -313,8 +316,6 @@ def slice_cube_window(
     deadtraces: bool = True,
 ):
     """Main entry point towards caller"""
-    logger.info("Inactive snapxy, has no role in algorithm 3 %s", snapxy)
-
     if showprogress:
         print("progress: initialising for attributes...")
 
@@ -326,16 +327,17 @@ def slice_cube_window(
 
     outside_proxy = None
     if not mask:
-        outside_proxy = _proxy_surf_outside_cube(_regsurf_, cube, deadtraces)
+        outside_proxy = _proxy_surf_outside_cube(_regsurf_, cube)
 
     masksurf = _get_iso_maskthreshold_surface(upper, lower, maskthreshold)
 
     dvalues = _create_depth_cube(cube)
 
+    apply_ndiv = 1 if sampling == "cube" else ndiv
     if showprogress:
-        print("progress: refine according to ndiv...")
+        print(f"progress: refine according to actual ndiv = {apply_ndiv}...")
 
-    cref, dref = _refine_cubes_vertically(cvalues, dvalues, ndiv)
+    cref, dref = _refine_cubes_vertically(cvalues, dvalues, apply_ndiv)
 
     cref = _filter_cube_values_upper_lower(cref, dref, upper, lower)
 
@@ -344,19 +346,21 @@ def slice_cube_window(
     use_attrs = _expand_attributes(attribute)
 
     attrs = {}
+    if showprogress:
+        print("progress: compute mean, variance, etc attributes...")
     for attr in use_attrs:
         if attr in SUMATTRS:
-            # use cvalues, which is not refined vertically
-            if showprogress:
-                print("progress: compute sum attributes...")
-            res = _compute_stats(cval, attr, _regsurf_, upper, masksurf, sampling)
+            # use cval, which is not refined vertically
+            res = _compute_stats(
+                cval, attr, _regsurf_, upper, masksurf, sampling, snapxy
+            )
         else:
-            if showprogress:
-                print("progress: compute mean, variance, etc attributes...")
-            res = _compute_stats(cref, attr, _regsurf_, upper, masksurf, sampling)
+            res = _compute_stats(
+                cref, attr, _regsurf_, upper, masksurf, sampling, snapxy
+            )
 
-        if outside_proxy:
-            res.values = np.ma.where(outside_proxy == 1, _regsurf_.values, res.values)
+        if outside_proxy and not snapxy:
+            res.values = np.ma.where(outside_proxy > 0, _regsurf_.values, res.values)
 
         attrs[attr] = res
 

src/xtgeo/surface/regular_surface.py

@@ -2638,8 +2638,11 @@ def slice_cube_window(
             maskthreshold (float): Only if two surface; if isochore is less
                 than given value, the result will be masked.
             snapxy (bool): If True (optional), then the map values will get
-                values at nearest Cube XY location. Only relevant to use if
-                surface is derived from seismic coordinates (e.g. Auto4D).
+                values at nearest Cube XY location, and the resulting surfaces layout
+                map will be defined by the seismic layout. Quite relevant to use if
+                surface is derived from seismic coordinates (e.g. Auto4D), but can be
+                useful in other cases also, as long as one notes that map definition
+                may change from input.
             showprogress (bool): If True, then a progress is printed to stdout.
             deadtraces (bool): If True (default) then dead cube traces
                 (given as value 2 in SEGY trace headers), are treated as
@@ -2680,6 +2683,12 @@ def slice_cube_window(
         if other is None and zrange is None:
             zrange = 10
 
+        if algorithm != 3:
+            warnings.warn(
+                "Other algorithms than no. 3 may be deprecated from xtgeo version 4",
+                PendingDeprecationWarning,
+            )
+
         asurfs = _regsurf_cube_window.slice_cube_window(
             self,
             cube,