diff --git a/src/clib/xtg/surf_cube_attr_intv.c b/src/clib/xtg/surf_cube_attr_intv.c index 72bb2d698..446dc9feb 100644 --- a/src/clib/xtg/surf_cube_attr_intv.c +++ b/src/clib/xtg/surf_cube_attr_intv.c @@ -145,7 +145,7 @@ pillar_attrs(double *pillar, int np, double *pattr) static void pillar_attrs_sums(double *pillar, int np, double *pattr) { - /* SUM attributes are only menaingful on a discrete level*/ + /* SUM attributes are only meaningful on a discrete level*/ int n; for (n = 11; n <= 13; n++) 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..7bc6b48c0 --- /dev/null +++ b/src/xtgeo/surface/_regsurf_cube_window_v3.py @@ -0,0 +1,506 @@ +# -*- coding: utf-8 -*- +"""Regular surface vs Cube, slice a window interval v3, in pure numpy""" + +from __future__ import annotations + +from typing import Type + +import numpy as np +from scipy.interpolate import interp1d + +import xtgeo +from xtgeo.common import XTGeoDialog + +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 + + +# TMP +# pylint: disable=W0613, W0612 + + +def _cut_cube_deadtraces(cube: Type[xtgeo.Cube]) -> Type[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 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") + + return dvalues + + +def _upper_lower_surface( + self, + cube, + zsurf: Type[xtgeo.RegularSurface], + other: Type[xtgeo.RegularSurface], + other_position, + zrange, +) -> list: + """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): + """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): + """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 + + return useattrs + + +def _compute_stats(cref, attr, self, upper): + """Compute the attribute and return the attribute map(s)""" + logger.info("Compute stats...") + + if attr == "mean": + values = np.nanmean(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 == "var": + values = np.var(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) + elif attr == "sumpos": + use = cref.copy() + use[cref < 0] = np.nan + values = np.nansum(use, axis=2) + 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 = values + actual.resample(sampled) + + logger.info("Compute stats... done") + return actual + + +def slice_cube_window( + self, + cube, + zsurf=None, + other=None, + other_position="below", + sampling="nearest", + mask=True, + zrange=10, + ndiv=None, + attribute="max", + maskthreshold=0.1, + snapxy=False, + showprogress=False, + deadtraces=True, +): + """Main entrypoint towards caller""" + logger.info("Inactive snapxy %s", snapxy) + + cvalues = _cut_cube_deadtraces(cube) + upper, lower = _upper_lower_surface( + self, cube, zsurf, other, other_position, zrange + ) + + dvalues = _create_depth_cube(cube) + + cref, dref = _refine_cubes_vertically(cvalues, dvalues, 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 = {} + for attr in use_attrs: + if attr in SUMATTRS: + # use cvalues, which is not refined vertically + res = _compute_stats(cvalues, attr, self, upper) + else: + res = _compute_stats(cref, attr, self, upper) + + 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 + + +# def _slice_cube_window( +# self, +# cube, +# zsurf, +# other, +# other_position, +# sampling, +# mask, +# zrange, +# ndiv, +# attribute, +# maskthreshold, +# showprogress, +# deadtraces, +# ): # pylint: disable=too-many-branches, too-many-statements +# """Slice Cube between surfaces to find attributes + +# New from Sept 2023, simplify using pure numpy operations +# """ +# logger.info("Slice cube window method v3") + +# olddead = None +# if deadtraces: +# olddead = cube.values_dead_traces(xtgeo.UNDEF) + +# optprogress = 0 +# if showprogress: +# optprogress = 1 + +# if not isinstance(attribute, list): +# if attribute == "all": +# attrlist = ALLATTRS +# else: +# attrlist = [attribute] +# else: +# attrlist = attribute + +# if zsurf is not None: +# this = zsurf +# else: +# this = self.copy() + +# if other is not None: +# zdelta = np.absolute(this.values - other.values) +# zdiameter = zdelta.max() +# if other_position.lower() == "below": +# surf1 = this +# surf2 = other +# else: +# surf1 = other # BEWARE self +# surf2 = this +# else: +# surf1 = this.copy() +# surf2 = this.copy() +# surf1.values -= zrange +# surf2.values += zrange +# zdiameter = 2 * zrange + +# if ndiv is None: +# ndiv = int(2 * zdiameter / cube.zinc) +# if ndiv < 1: +# ndiv = 1 +# logger.warning("NDIV < 1; reset to 1") + +# # force/overrule a coarse sampling for sampling option "cube" +# ndivdisc = int((zdiameter) * 1.0001 / cube.zinc) +# if sampling == "cube": +# ndiv = ndivdisc + +# zrinc = zdiameter / ndiv +# logger.debug("zdiameter and cube zinc: %s %s", zdiameter, cube.zinc) +# logger.debug("zrinc and ndiv: %s %s", zrinc, ndiv) + +# optsum = 0 +# if any("sum" in word for word in attrlist): +# optsum = 1 + +# results = _attributes_betw_surfaces( +# self, +# cube, +# surf1, +# surf2, +# sampling, +# mask, +# zrinc, +# ndiv, +# ndivdisc, +# optprogress, +# maskthreshold, +# optsum, +# ) + +# if deadtraces: +# cube.values_dead_traces(olddead) # reset value for dead traces + +# # build the returning result +# alist = {} +# for att in attrlist: +# num = ALLATTRS.index(att) +# alist[att] = self.copy() +# alist[att].values = results[num, :] + +# return alist + + +# def _attributes_betw_surfaces( +# self, +# cube, +# surf1, +# surf2, +# sampling, +# maskopt, +# zrinc, +# ndiv, +# ndivdisc, +# optprogress, +# maskthreshold, +# optsum, +# ): +# """This is the actual lowlevel engine communicating with C code""" + +# logger.info("Attributes between surfaces") + +# results = np.zeros((len(ALLATTRS) * self.ncol * self.nrow), dtype=np.float64) + +# optnearest = 0 +# if sampling in ["nearest", "cube"]: +# optnearest = 1 + +# _cxtgeo.surf_cube_attr_intv( +# cube.ncol, +# cube.nrow, +# cube.nlay, +# cube.zori, +# cube.zinc, +# cube.values, +# surf1.values.data, +# surf2.values.data, +# surf1.values.mask, +# surf2.values.mask, +# zrinc, +# ndiv, +# ndivdisc, +# results, +# optnearest, +# maskopt, +# optprogress, +# maskthreshold, +# optsum, +# ) + +# logger.info("Results updated, with size %s", results.shape) + +# results = results.reshape((len(ALLATTRS), self.ncol * self.nrow), order="C") + +# return results + + +# def _slice_cube_window_resample( +# self, +# cube, +# zsurf, +# other, +# other_position, +# sampling, +# mask, +# zrange, +# ndiv, +# attribute, +# maskthreshold, +# showprogress, +# deadtraces, +# ): +# """Makes a resample from original surfaces first to fit cube topology""" + +# logger.info("Attributes between surfaces, resampling version") + +# scube = xtgeo.surface_from_cube(cube, 0.0) + +# scube.resample(self) + +# szsurf = None +# if zsurf: +# szsurf = scube.copy() +# szsurf.resample(zsurf) + +# sother = None +# if other: +# sother = scube.copy() +# sother.resample(other) + +# attrs = _slice_cube_window( +# scube, +# cube, +# szsurf, +# sother, +# other_position, +# sampling, +# mask, +# zrange, +# ndiv, +# attribute, +# maskthreshold, +# showprogress, +# deadtraces, +# ) + +# # now resample attrs back to a copy of self +# zelf = self.copy() +# for key, _val in attrs.items(): +# zelf.resample(attrs[key]) +# attrs[key] = zelf.copy() + +# return attrs diff --git a/tests/test_surface/test_surface_cube_attrs.py b/tests/test_surface/test_surface_cube_attrs.py index 7f9bc504c..eb36204ec 100644 --- a/tests/test_surface/test_surface_cube_attrs.py +++ b/tests/test_surface/test_surface_cube_attrs.py @@ -2,6 +2,7 @@ from os.path import join import pytest + import xtgeo xtg = xtgeo.common.XTGeoDialog() @@ -155,6 +156,54 @@ def test_various_attrs_algorithm2(loadsfile1): assert surfx.values.mean() == pytest.approx(529.328, abs=0.01) +def test_various_attrs_algorithm3(loadsfile1): + cube1 = loadsfile1 + surf1 = xtgeo.surface_from_cube(cube1, 2540) + surf2 = xtgeo.surface_from_cube(cube1, 2548) + + surfx = surf1.copy() + surfx.slice_cube_window( + cube1, + other=surf2, + other_position="below", + attribute="mean", + sampling="trilinear", + snapxy=True, + ndiv=None, + algorithm=3, + ) + + assert surfx.values.mean() == pytest.approx(177.34, abs=0.1) + + surfx = surf1.copy() + surfx.slice_cube_window( + cube1, + other=surf2, + other_position="below", + attribute="maxneg", + sampling="trilinear", + snapxy=True, + ndiv=None, + algorithm=3, + ) + + assert surfx.values.count() == 0 # no nonmasked elements + + surfx = surf1.copy() + surfx.slice_cube_window( + cube1, + other=surf2, + other_position="below", + attribute="sumabs", + sampling="cube", + snapxy=True, + ndiv=None, + algorithm=3, + ) + + assert surfx.values.mean() == pytest.approx(529.328, abs=0.01) + + def test_avg_surface(loadsfile1): cube1 = loadsfile1 surf1 = xtgeo.surface_from_cube(cube1, 1100.0) @@ -237,8 +286,10 @@ def test_avg_surface2(loadsfile1): @pytest.mark.benchmark(group="cube slicing") -@pytest.mark.parametrize("algorithm", [pytest.param(1, marks=pytest.mark.bigtest), 2]) -def test_avg_surface_large_cube_algorithm1(benchmark, algorithm): +@pytest.mark.parametrize( + "algorithm", [pytest.param(1, marks=pytest.mark.bigtest), 2, 3] +) +def test_avg_surface_large_cube_algorithmx(benchmark, algorithm): cube1 = xtgeo.Cube(ncol=120, nrow=120, nlay=100, zori=200, xinc=12, yinc=12, zinc=4) cube1.values[400:80, 400:80, :] = 12 @@ -270,7 +321,7 @@ def test_attrs_reek(tmpdir, loadsfile2): t2a = xtgeo.surface_from_file(TOP2A) t2b = xtgeo.surface_from_file(TOP2B) - attlist = ["maxpos", "maxneg"] + attlist = ["maxpos", "maxneg", "mean", "rms"] attrs1 = t2a.slice_cube_window( cube2, other=t2b, sampling="trilinear", attribute=attlist, algorithm=1 @@ -278,12 +329,19 @@ def test_attrs_reek(tmpdir, loadsfile2): attrs2 = t2a.slice_cube_window( cube2, other=t2b, sampling="trilinear", attribute=attlist, algorithm=2 ) + attrs3 = t2a.slice_cube_window( + cube2, other=t2b, sampling="trilinear", attribute=attlist, algorithm=3 + ) - for att in attrs1.keys(): + for att, _ in attrs1.items(): srf1 = attrs1[att] srf2 = attrs2[att] + srf3 = attrs3[att] srf1.to_file(join(tmpdir, "attr1_" + att + ".gri")) srf2.to_file(join(tmpdir, "attr2_" + att + ".gri")) + srf3.to_file(join(tmpdir, "attr3_" + att + ".gri")) assert srf1.values.mean() == pytest.approx(srf2.values.mean(), abs=0.005) + assert srf3.values.mean() == pytest.approx(srf2.values.mean(), abs=0.005) + print("\nok")