diff --git a/src/xtgeo/cube/cube1.py b/src/xtgeo/cube/cube1.py index 1f6b7aa32..074eef23f 100644 --- a/src/xtgeo/cube/cube1.py +++ b/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) diff --git a/src/xtgeo/grid3d/grid.py b/src/xtgeo/grid3d/grid.py index 170df3373..52d083fe5 100644 --- a/src/xtgeo/grid3d/grid.py +++ b/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 diff --git a/src/xtgeo/surface/_regsurf_cube_window_v3.py b/src/xtgeo/surface/_regsurf_cube_window_v3.py index 2754a6d2c..4c1265ffc 100644 --- a/src/xtgeo/surface/_regsurf_cube_window_v3.py +++ b/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 diff --git a/src/xtgeo/surface/regular_surface.py b/src/xtgeo/surface/regular_surface.py index fb6a894f0..d67ac6b54 100644 --- a/src/xtgeo/surface/regular_surface.py +++ b/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,