From 9829ac09c799e79c5a5aad59be2570f0452e8c5c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jan=20C=2E=20Riven=C3=A6s?= Date: Thu, 21 Sep 2023 08:33:40 +0200 Subject: [PATCH] ENH: add algorithm 3 for slice_cube_window() This algorithm is purely in Python which makes debugging significantly easier. For niv > 1, the speed compared with algorithm 2 will be somewhat decreased, but we can later see if it is possible to speedup the refining of an existing cube. --- src/xtgeo/surface/_regsurf_cube_window.py | 27 +- src/xtgeo/surface/_regsurf_cube_window_v3.py | 365 +++++++++++++++++++ src/xtgeo/surface/regular_surface.py | 102 ++++-- 3 files changed, 455 insertions(+), 39 deletions(-) create mode 100644 src/xtgeo/surface/_regsurf_cube_window_v3.py diff --git a/src/xtgeo/surface/_regsurf_cube_window.py b/src/xtgeo/surface/_regsurf_cube_window.py index 37d87022d..572362b1d 100644 --- a/src/xtgeo/surface/_regsurf_cube_window.py +++ b/src/xtgeo/surface/_regsurf_cube_window.py @@ -4,9 +4,11 @@ import numpy as np import numpy.ma as ma -from xtgeo.common import XTGeoDialog -from xtgeo.common import XTGShowProgress -from . import _regsurf_cube_window_v2 as cwv + +from xtgeo.common import XTGeoDialog, XTGShowProgress + +from . import _regsurf_cube_window_v2 as cwv2 +from . import _regsurf_cube_window_v3 as cwv3 xtg = XTGeoDialog() @@ -69,8 +71,25 @@ def slice_cube_window( deadtraces=deadtraces, ) + elif algorithm == 2: + attrs = cwv2.slice_cube_window( + self, + cube, + zsurf=zsurf, + other=other, + other_position=other_position, + sampling=sampling, + mask=mask, + zrange=zrange, + ndiv=ndiv, + attribute=attribute, + maskthreshold=maskthreshold, + snapxy=snapxy, + showprogress=showprogress, + deadtraces=deadtraces, + ) else: - attrs = cwv.slice_cube_window( + attrs = cwv3.slice_cube_window( self, cube, zsurf=zsurf, diff --git a/src/xtgeo/surface/_regsurf_cube_window_v3.py b/src/xtgeo/surface/_regsurf_cube_window_v3.py new file mode 100644 index 000000000..20e23b737 --- /dev/null +++ b/src/xtgeo/surface/_regsurf_cube_window_v3.py @@ -0,0 +1,365 @@ +# -*- coding: utf-8 -*- +"""Regular surface vs Cube, slice a window interval v3, in pure numpy.""" + +from __future__ import annotations + +import warnings +from typing import Optional, Tuple, Union + +import numpy as np +from scipy.interpolate import interp1d + +import xtgeo +from xtgeo.common import XTGeoDialog + +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() + +logger = xtg.functionlogger(__name__) + + +STATATTRS = [ + "min", + "max", + "mean", + "var", + "rms", + "maxpos", + "maxneg", + "maxabs", + "meanabs", + "meanpos", + "meanneg", +] +SUMATTRS = [ + "sumpos", + "sumneg", + "sumabs", +] +ALLATTRS = STATATTRS + SUMATTRS + +# self ~ RegularSurface() instance + + +def _cut_cube_deadtraces(cube: xtgeo.Cube, deadtraces: bool) -> np.ndarray: + """Take the cube numpy values and filter away dead traces as np.nan.""" + logger.info("Assign dead traces") + dvalues = cube.values.copy() + + if deadtraces and 2 in cube.traceidcodes: + dvalues[cube.traceidcodes == 2] = np.nan + logger.info("Dead traces encountered in this cube, set values to np.nan") + else: + logger.info("No dead traces encountered in this cube, or deadtraces is False") + + return dvalues + + +def _get_iso_maskthreshold_surface( + upper: xtgeo.RegularSurface, + lower: xtgeo.RegularSurface, + maskthreshold: float, +) -> xtgeo.RegularSurface: + """Return a surface with value 0 where isochore <= threshold""" + logger.info("Maskthreshold based on isochore") + result = upper.copy() + result.fill() + diff = lower - upper + result.values = np.ma.where(diff.values <= maskthreshold, 0, 1) + return result + + +def _proxy_surf_outside_cube( + self, + 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 = self.copy() + outside.values = 0.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 + + +def _upper_lower_surface( + self, + cube: xtgeo.Cube, + zsurf: xtgeo.RegularSurface, + other: xtgeo.RegularSurface, + other_position: str, + zrange: float, +) -> Tuple[xtgeo.RegularSurface, xtgeo.RegularSurface]: + """Return upper and lower surface, sampled to cube resolution.""" + + logger.info("Define surfaces to apply...") + this = zsurf if zsurf is not None else self.copy() + + if other is not None: + if other_position.lower() == "below": + surf1 = this + surf2 = other + else: + surf1 = other # avoid changing self instance + surf2 = this + else: + surf1 = this.copy() + surf2 = this.copy() + surf1.values -= zrange + surf2.values += zrange + + # get the surfaces on cube resolution + upper = xtgeo.surface_from_cube(cube, 0) + lower = upper.copy() + upper.resample(surf1) + lower.resample(surf2) + + logger.info( + "Return resmapled surfaces, avg for upper and lower is %s, %s", + upper.values.mean(), + lower.values.mean(), + ) + return upper, lower + + +def _create_depth_cube(cube: xtgeo.Cube) -> np.ndarray: + """Create a cube (np array) where values are cube depths; to be used as filter.""" + logger.info("Create a depth cube...") + dcube = cube.values.copy() + darr = [cube.zori + n * cube.zinc for n in range(dcube.shape[2])] + dcube[:, :, :] = darr + + logger.info("Created a depth cube starting from %s", np.mean(dcube)) + + return dcube + + +def _refine_cubes_vertically(cvalues, dvalues, ndiv): + """Refine the cubes vertically for better resolution""" + if not ndiv: + ndiv = 2 # default + logger.info("Resampling vertically, according to ndiv = %s", ndiv) + if ndiv <= 1: + logger.info("ndiv is less or equal to 1; no refinement done") + return cvalues, dvalues + + logger.info("Original shape is %s", cvalues.shape) + cref = np.random.rand(cvalues.shape[0], cvalues.shape[1], ndiv * cvalues.shape[2]) + dref = cref.copy() + + num_points = cref.shape[-1] + # Create interpolation function for cube values + interp_func1 = interp1d( + np.arange(cvalues.shape[-1]), + cvalues, + axis=-1, + kind="linear", + fill_value="extrapolate", + ) + interp_func2 = interp1d( + np.arange(dvalues.shape[-1]), + dvalues, + axis=-1, + kind="linear", + fill_value="extrapolate", + ) + # Resample array2 to match the number of points in array1 + cref = interp_func1(np.linspace(0, cvalues.shape[-1] - 1, num_points)) + dref = interp_func2(np.linspace(0, dvalues.shape[-1] - 1, num_points)) + + logger.info("Resampling done, new shape is %s", cref.shape) + return cref, dref + + +def _filter_cube_values_upper_lower(cvalues, dvalues, upper, lower): + """Filter the cube (cvalues) based on depth interval.""" + + nnans = np.count_nonzero(np.isnan(cvalues)) + logger.info("Filter cube in depth... number of nans is %s", nnans) + + upv = np.expand_dims(upper.values, axis=2) + lov = np.expand_dims(lower.values, axis=2) + + cvalues[dvalues < upv] = np.nan + cvalues[dvalues > lov] = np.nan + + nnans = np.count_nonzero(np.isnan(cvalues)) + ndefi = np.count_nonzero(~np.isnan(cvalues)) + logger.info("Filter cube in depth done, updated number of nans is %s", nnans) + logger.info("Filter cube in depth done, remaining is %s", ndefi) + return cvalues + + +def _expand_attributes(attribute: Union[str, list]) -> list: + """The 'attribute' may be a name, 'all', or a list of attributes""" + useattrs = None + if isinstance(attribute, str): + if attribute == "all": + useattrs = ALLATTRS + else: + useattrs = [attribute] + else: + useattrs = attribute + + if not all(item in ALLATTRS for item in useattrs): + raise ValueError( + f"One or more values are not a valid, input list is {useattrs}, " + f"allowed list is {ALLATTRS}" + ) + return useattrs + + +def _compute_stats( + cref: np.ndarray, + attr: str, + self: xtgeo.RegularSurface, + upper: xtgeo.RegularSurface, + masksurf: xtgeo.RegularSurface, + sampling: str, + snapxy: bool, +) -> xtgeo.RegularSurface: + """Compute the attribute and return the attribute map""" + logger.info("Compute stats...") + + if attr == "mean": + values = np.nanmean(cref, axis=2) + elif attr == "var": + values = np.nanvar(cref, axis=2) + elif attr == "max": + values = np.nanmax(cref, axis=2) + elif attr == "min": + values = np.nanmin(cref, axis=2) + elif attr == "rms": + values = np.sqrt(np.nanmean(np.square(cref), axis=2)) + elif attr == "maxneg": + use = cref.copy() + use[cref >= 0] = np.nan + values = np.nanmin(use, axis=2) + elif attr == "maxpos": + use = cref.copy() + use[cref < 0] = np.nan + values = np.nanmax(use, axis=2) + elif attr == "maxabs": + use = cref.copy() + use = np.abs(use) + values = np.nanmax(use, axis=2) + elif attr == "meanabs": + use = cref.copy() + use = np.abs(use) + values = np.nanmean(use, axis=2) + elif attr == "meanpos": + use = cref.copy() + use[cref < 0] = np.nan + values = np.nanmean(use, axis=2) + elif attr == "meanneg": + use = cref.copy() + use[cref >= 0] = np.nan + values = np.nanmean(use, axis=2) + elif attr == "sumneg": + use = cref.copy() + use[cref >= 0] = np.nan + values = np.nansum(use, axis=2) + values = np.ma.masked_greater_equal(values, 0.0) # to make undefined map areas + elif attr == "sumpos": + use = cref.copy() + use[cref < 0] = np.nan + values = np.nansum(use, axis=2) + values = np.ma.masked_less_equal(values, 0.0) # to make undefined map areas + elif attr == "sumabs": + use = cref.copy() + use = np.abs(use) + values = np.nansum(use, axis=2) + else: + raise ValueError(f"The attribute name {attr} is not supported") + + actual = self.copy() + sampled = upper.copy() + + sampled.values = np.ma.masked_invalid(values) + sampled.values = np.ma.masked_where(masksurf.values == 0, sampled.values) + + 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 + + +def slice_cube_window( + self, + cube: xtgeo.Cube, + zsurf: Optional[xtgeo.RegularSurface] = None, + other: Optional[xtgeo.RegularSurface] = None, + other_position: str = "below", + sampling: str = "nearest", + mask: bool = True, + zrange: float = 10.0, + ndiv: Optional[int] = None, + attribute: str = "max", + maskthreshold: float = 0.1, + snapxy: bool = False, + showprogress: bool = False, + deadtraces: bool = True, +): + """Main entry point towards caller""" + if showprogress: + print("progress: initialising for attributes...") + + cvalues = _cut_cube_deadtraces(cube, deadtraces) + + upper, lower = _upper_lower_surface( + self, cube, zsurf, other, other_position, zrange + ) + + outside_proxy = None + if not mask: + outside_proxy = _proxy_surf_outside_cube(self, 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(f"progress: refine according to actual ndiv = {apply_ndiv}...") + + cref, dref = _refine_cubes_vertically(cvalues, dvalues, apply_ndiv) + + cref = _filter_cube_values_upper_lower(cref, dref, upper, lower) + + cval = _filter_cube_values_upper_lower(cvalues, dvalues, upper, lower) + + 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 cval, which is not refined vertically + res = _compute_stats(cval, attr, self, upper, masksurf, sampling, snapxy) + else: + res = _compute_stats(cref, attr, self, upper, masksurf, sampling, snapxy) + + if outside_proxy and not snapxy: + res.values = np.ma.where(outside_proxy > 0, self.values, res.values) + + attrs[attr] = res + + # if attribute is str, self shall be updated and None returned, + # otherwise a dict of attributes objects shall be returned + if isinstance(attrs, dict) and len(attrs) == 1 and isinstance(attribute, str): + self.values = attrs[attribute].values + return None + + return attrs diff --git a/src/xtgeo/surface/regular_surface.py b/src/xtgeo/surface/regular_surface.py index 2348772fc..ea81dd6b9 100644 --- a/src/xtgeo/surface/regular_surface.py +++ b/src/xtgeo/surface/regular_surface.py @@ -32,6 +32,7 @@ # -------------------------------------------------------------------------------------- # pylint: disable=too-many-public-methods +from __future__ import annotations import functools import io @@ -42,7 +43,7 @@ from collections import OrderedDict from copy import deepcopy from types import FunctionType -from typing import List, Optional, Tuple, Type, Union +from typing import Dict, List, Literal, Optional, Tuple, Type, Union import deprecation import numpy as np @@ -69,6 +70,22 @@ xtg = xtgeo.common.XTGeoDialog() logger = xtg.functionlogger(__name__) +# valid argumentts for seismic attributes +ValidAttrs = Literal[ + "all", + "max", + "min", + "rms", + "mean", + "var", + "maxpos", + "maxneg", + "sumpos", + "sumneg", + "meanabs", + "meanpos", + "meanneg", +] # ====================================================================================== # METHODS as wrappers to class init + import @@ -2541,21 +2558,21 @@ def slice_cube( def slice_cube_window( self, - cube, - zsurf=None, - other=None, - other_position="below", - sampling="nearest", - mask=True, - zrange=None, - ndiv=None, - attribute="max", - maskthreshold=0.1, - snapxy=False, - showprogress=False, - deadtraces=True, - algorithm=2, - ): + cube: xtgeo.Cube, + zsurf: Optional[xtgeo.RegularSurface] = None, + other: Optional[xtgeo.RegularSurface] = None, + other_position: str = "below", + sampling: Literal["nearest", "cube", "trilinear"] = "nearest", + mask: bool = True, + zrange: Optional[float] = None, + ndiv: Optional[int] = None, + attribute: Union[List[ValidAttrs], ValidAttrs] = "max", + maskthreshold: float = 0.1, + snapxy: bool = False, + showprogress: bool = False, + deadtraces: bool = True, + algorithm: Literal[1, 2, 3] = 2, + ) -> Optional[Dict[xtgeo.RegularSurface]]: """Slice the cube within a vertical window and get the statistical attrubutes. The statistical attributes can be min, max etc. Attributes are: @@ -2590,44 +2607,49 @@ def slice_cube_window( available. Args: - cube (Cube): Instance of a Cube() - zsurf (RegularSurface): Instance of a depth (or time) map, which + cube: Instance of a Cube() here + zsurf: Instance of a depth (or time) map, which is the depth or time map (or...) that is used a slicer. If None, then the surface instance itself is used a slice criteria. Note that zsurf must have same map defs as the surface instance. - other (RegularSurface): Instance of other surface if window is + other: Instance of other surface if window is between surfaces instead of a static window. The zrange input is then not applied. - sampling (str): 'nearest'/'trilinear'/'cube' for nearest node (default), + sampling: 'nearest'/'trilinear'/'cube' for nearest node (default), or 'trilinear' for trilinear interpolation. The 'cube' option is only available with algorithm = 2 and will overrule ndiv and sample at the cube's Z increment resolution. - mask (bool): If True (default), then the map values outside + mask: If True (default), then the map values outside the cube will be undef, otherwise map will be kept as-is - zrange (float): The one-sided "radius" range of the window, e.g. 10 + zrange: The one-sided "radius" range of the window, e.g. 10 (10 is default) units (e.g. meters if in depth mode). The full window is +- zrange (i.e. diameter). - If other surface is present, zrange is computed based on that. - ndiv (int): Number of intervals for sampling within zrange. None + If other surface is present, zrange is computed based on those + two surfaces instead. + ndiv: Number of intervals for sampling within zrange. None means 'auto' sampling, using 0.5 of cube Z increment as basis. If - algorithm = 2 and sampling is 'cube', the cube Z increment + algorithm = 2/3 and sampling is 'cube', the cube Z increment will be used. - attribute (str or list): The requested attribute(s), e.g. + attribute: The requested attribute(s), e.g. 'max' value. May also be a list of attributes, e.g. ['min', 'rms', 'max']. By such, a dict of surface objects is - returned. Note 'all' will make a list of possible attributes + returned. Note 'all' will make a list of all possible attributes. 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). - showprogress (bool): If True, then a progress is printed to stdout. - deadtraces (bool): If True (default) then dead cube traces + snapxy: If True (optional), then the map values will get + 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: If True, then a progress is printed to stdout. + deadtraces: If True (default) then dead cube traces (given as value 2 in SEGY trace headers), are treated as - undefined, nad map will be undefined at dead trace location. - algorithm (int): 1 for legacy method, 2 (default) for new faster - and more precise method available from xtgeo version 2.9. + undefined, and map will be undefined at dead trace location. + algorithm: 1 for legacy method, 2 (default) for new faster + and more precise method available from xtgeo version 2.9, and + algorithm 3 as new implementation from Sept. 2023 (v3.4) Example:: @@ -2657,10 +2679,20 @@ def slice_cube_window( .. versionchanged:: 2.9 Added ``algorithm`` keyword, default is now 2, while 1 is the legacy version + + .. versionchanged:: 3.4 Added ``algorithm`` 3 which is more robust and + hence recommended! + """ 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,