Merge pull request equinor#458 from jcrivenaes/wells-shoulderbed-filter

ENH: add mask_shoulderbeds() for Well(), issue equinor#457

.pylintrc

@@ -8,7 +8,7 @@ good-names=logger, version, xtg, i, j, k, x, y, z, _
 additional-builtins=_x, _y, _z, _tmp1, _tmp2
 variable-rgx=^[a-z_][_a-z0-9]+((_[a-z0-9]+)*)?$
 argument-rgx=^[a-z_][_a-z0-9]+((_[a-z0-9]+)*)?$
-dummy-variables-rgx=^_+[a-z0-9]*?$)|dummy
+dummy-variables-rgx=^_+[a-z0-9]*?$|dummy
 
 [TYPECHECK]
 generated-members=np.*, numpy.*, pd.*, pandas.*, cxtgeo.*, matplotlib.*

src/xtgeo/well/_well_oper.py

@@ -1,6 +1,4 @@
-# -*- coding: utf-8 -*-
-"""Operations along a well, private module"""
-
+"""Operations along a well, private module."""
 
 import copy
 from distutils.version import StrictVersion
@@ -20,8 +18,8 @@
 
 def delete_log(self, lname):
     """Delete/remove an existing log, or list of logs."""
-
     self._ensure_consistency()
+
     if not isinstance(lname, list):
         lname = [lname]
 
@@ -49,11 +47,10 @@ def delete_log(self, lname):
 
 
 def rescale(self, delta=0.15, tvdrange=None):
-    """Rescale by using a new MD increment
+    """Rescale by using a new MD increment.
 
     The rescaling is technically done by interpolation in the Pandas dataframe
     """
-
     pdrows = pd.options.display.max_rows
     pd.options.display.max_rows = 999
 
@@ -122,8 +119,7 @@ def rescale(self, delta=0.15, tvdrange=None):
 
 
 def make_zone_qual_log(self, zqname):
-    """Make a flag log based on stratigraphic relations"""
-
+    """Make a flag log based on stratigraphic relations."""
     if zqname in self.dataframe:
         logger.warning("Quality log %s exists, will be overwritten", zqname)
 
@@ -186,7 +182,7 @@ def make_zone_qual_log(self, zqname):
 
 
 def make_ijk_from_grid(self, grid, grid_id="", algorithm=1, activeonly=True):
-
+    """Make an IJK log from grid indices."""
     logger.info("Using algorithm %s in %s", algorithm, __name__)
 
     if algorithm == 1:
@@ -198,11 +194,9 @@ def make_ijk_from_grid(self, grid, grid_id="", algorithm=1, activeonly=True):
 
 
 def _make_ijk_from_grid_v1(self, grid, grid_id=""):
-    """
-    Getting IJK from a grid and make as well logs.
+    """Getting IJK from a grid and make as well logs.
 
     This is the first version, using _cxtgeo.grd3d_well_ijk from C
-
     """
     logger.info("Using algorithm 1 in %s", __name__)
 
@@ -274,14 +268,12 @@ def _make_ijk_from_grid_v1(self, grid, grid_id=""):
 
 
 def _make_ijk_from_grid_v2(self, grid, grid_id="", activeonly=True):
-    """
-    Getting IJK from a grid and make as well logs.
+    """Getting IJK from a grid and make as well logs.
 
     This is a newer version, using grid.get_ijk_from_points which in turn
     use the from C method x_chk_point_in_hexahedron, while v1 use the
     x_chk_point_in_cell. This one is believed to be more precise!
     """
-
     # establish a Points instance and make points dataframe from well trajectory X Y Z
     wpoints = xtgeo.Points()
     wpdf = self.dataframe.loc[:, ["X_UTME", "Y_UTMN", "Z_TVDSS"]].copy()
@@ -314,8 +306,7 @@ def _make_ijk_from_grid_v2(self, grid, grid_id="", activeonly=True):
 
 
 def get_gridproperties(self, gridprops, grid=("ICELL", "JCELL", "KCELL"), prop_id=""):
-    """Getting gridproperties as logs"""
-
+    """Getting gridproperties as logs."""
     if not isinstance(gridprops, (xtgeo.GridProperty, xtgeo.GridProperties)):
         raise ValueError('"gridprops" not a GridProperties or GridProperty instance')
 
@@ -438,3 +429,67 @@ def report_zonation_holes(self, threshold=5):
         clm = ["INDEX", "X_UTME", "Y_UTMN", "Z_TVDSS", "Zone", "Well"]
 
     return pd.DataFrame(wellreport, columns=clm)
+
+
+def mask_shoulderbeds(self, inputlogs, targetlogs, nsamples, strict):
+    """Mask targetlogs around discrete boundaries."""
+    logger.info("Mask shoulderbeds for some logs...")
+
+    # check that inputlogs exists and that they are discrete, and targetlogs
+    useinputs = []
+    for inlog in inputlogs:
+        if inlog not in self._wlogtypes.keys() and strict is True:
+            raise ValueError(f"Input log {inlog} is missing and strict=True")
+        if inlog in self._wlogtypes.keys() and self._wlogtypes[inlog] != "DISC":
+            raise ValueError(f"Input log {inlog} is not of type DISC")
+        if inlog in self._wlogtypes.keys():
+            useinputs.append(inlog)
+
+    usetargets = []
+    for target in targetlogs:
+        if target not in self._wlogtypes.keys() and strict is True:
+            raise ValueError(f"Target log {target} is missing and strict=True")
+        if target in self._wlogtypes.keys():
+            usetargets.append(target)
+
+    maxlen = len(self._df) // 2
+    if not isinstance(nsamples, int) or nsamples < 1 or nsamples > maxlen:
+        raise ValueError(f"Keyword nsamples must be an int > 1 and < {maxlen}")
+
+    if not useinputs or not usetargets:
+        logger.info("Mask shoulderbeds for some logs... nothing done")
+        return False
+
+    for inlog in useinputs:
+        inseries = self._df[inlog]
+        bseries = _get_bseries(inseries, nsamples)
+        for target in usetargets:
+            self._df.loc[bseries, target] = np.nan
+
+    logger.info("Mask shoulderbeds for some logs... done")
+    return True
+
+
+def _get_bseries(inseries, nsamples):
+    """Private function for creating a bool filter, returning a bool series."""
+    if not isinstance(inseries, pd.Series):
+        raise RuntimeError("Bug, input must be a pandas Series() instance.")
+
+    if len(inseries) == 0:
+        return pd.Series([], dtype=bool)
+
+    # nsmaples < 1 or input series with <= 1 element will not be prosessed
+    if nsamples < 1 or len(inseries) <= 1:
+        return pd.Series(inseries, dtype=bool).replace(True, False)
+
+    def _growfilter(bseries, nleft):
+        if not nleft:
+            return bseries
+
+        return _growfilter(bseries | bseries.shift(-1) | bseries.shift(1), nleft - 1)
+
+    # make a tmp mask log (series) based input logs and use that for mask filterings
+    transitions = inseries.diff().abs() > 0
+    bseries = transitions | transitions.shift(-1)
+
+    return _growfilter(bseries, nsamples - 1)

src/xtgeo/well/well1.py

@@ -4,7 +4,7 @@
 import sys
 from copy import deepcopy
 from distutils.version import StrictVersion
-from typing import Union, Optional
+from typing import Union, Optional, List
 from pathlib import Path
 import io
 from collections import OrderedDict
@@ -1427,6 +1427,46 @@ def get_fraction_per_zone(
 
         return dfr
 
+    def mask_shoulderbeds(
+        self,
+        inputlogs: List[str],
+        targetlogs: List[str],
+        nsamples: Optional[int] = 2,
+        strict: Optional[bool] = False,
+    ) -> bool:
+        """Mask data around zone boundaries or other discrete log boundaries.
+
+        This operates on number of samples, hence the actual distance which is masked
+        depends on the sampling interval. In future versions, operating on distances
+        (e.g. in TVD (true vertical depth) or MD (measured depth) will be available.
+
+        Args:
+            inputlogs: List of input logs, must be of discrete type.
+            targetlogs: List of logs where mask is applied.
+            nsample: Number of samples around boundaries to filter, per side, i.e.
+                value 2 means 2 above and 2 below, in total 4 samples.
+            strict: If True, will raise Exception of any of the input or target log
+                names are missing.
+
+        Returns:
+            True if any operation has been done. False in case nothing has been done,
+                 e.g. no targetlogs for this particular well and ``strict`` is False.
+
+        Raises:
+            ValueError: Input log {inlog} is missing and strict=True.
+            ValueError: Input log {inlog} is not of type DISC.
+            ValueError: Target log {target} is missing and strict=True.
+            ValueError: Keyword nsamples must be an int > 1 and < {maxlen} (where
+                maxlen is half of number of log samples).
+
+        Example:
+            >>> mywell.mask_shoulderbeds(["ZONELOG", "FACIES"], ["PHIT", "KLOGH"])
+
+        """
+        return _well_oper.mask_shoulderbeds(
+            self, inputlogs, targetlogs, nsamples, strict
+        )
+
     def get_surface_picks(self, surf):
         """Return :class:`.Points` obj where well crosses the surface (horizon picks).
 

tests/test_well/test_well.py

@@ -43,13 +43,13 @@
 
 @pytest.fixture(name="loadwell1")
 def fixture_loadwell1():
-    """Fixture for loading a well (pytest setup)"""
+    """Fixture for loading a well (pytest setup)."""
     return Well(WFILE)
 
 
 @pytest.fixture(name="loadwell3")
 def fixture_loadwell3():
-    """Fixture for loading a well (pytest setup)"""
+    """Fixture for loading a well (pytest setup)."""
     return Well(WELL3)
 
 
@@ -310,8 +310,7 @@ def test_make_hlen(loadwell1):
 
 
 def test_make_zqual_log(loadwell3):
-    """Make a zonelog FLAG quality log"""
-
+    """Make a zonelog FLAG quality log."""
     mywell = loadwell3
     mywell.zonelogname = "ZONELOG"
 
@@ -323,6 +322,68 @@ def test_make_zqual_log(loadwell3):
         print(mywell.dataframe)
 
 
+@pytest.mark.parametrize(
+    "logseries, nsamples, expected",
+    [
+        ([], 0, []),
+        ([1], 0, [False]),
+        ([1], 1, [False]),  # intentional
+        ([1, 2], 0, [False, False]),
+        ([0, 0], 1, [False, False]),
+        ([1, 2], 1, [True, True]),
+        ([1, 1, 1, 2, 2, 2], 1, [False, False, True, True, False, False]),
+        ([1, 1, 1, 2, 2, 2], 2, [False, True, True, True, True, False]),
+        ([10, 10, 10, 2, 2, 2], 1, [False, False, True, True, False, False]),
+        ([1, 1, 1, 2, 2, 2], 10, [True, True, True, True, True, True]),
+        ([np.nan, 1, 1, np.nan, 2, 2], 1, [False, False, False, False, False, False]),
+        ([np.nan, 1, 1, 2, np.nan, 2], 1, [False, False, True, True, False, False]),
+    ],
+)
+def test_mask_shoulderbeds_get_bseries(logseries, nsamples, expected):
+    """Test for corner cases in _mask_shoulderbeds_logs."""
+    from xtgeo.well._well_oper import _get_bseries
+
+    logseries = pd.Series(logseries, dtype="float64")
+    expected = pd.Series(expected, dtype="bool")
+
+    results = _get_bseries(logseries, nsamples)
+    if logseries.empty:
+        assert expected.empty is True
+    else:
+        assert results.equals(expected)
+
+
+def test_mask_shoulderbeds(loadwell3, loadwell1):
+    """Test masking shoulderbeds effects."""
+    mywell = loadwell3
+
+    usewell = mywell.copy()
+    usewell.mask_shoulderbeds(["ZONELOG"], ["GR"], nsamples=3)
+
+    assert not np.isnan(mywell.dataframe.at[1595, "GR"])
+    assert np.isnan(usewell.dataframe.at[1595, "GR"])
+
+    # another well set with more discrete logs and several logs to modify
+    mywell = loadwell1
+    usewell = mywell.copy()
+    usewell.mask_shoulderbeds(["Zonelog", "Facies"], ["Perm", "Poro"], nsamples=2)
+    assert np.isnan(usewell.dataframe.at[4763, "Perm"])
+    assert np.isnan(usewell.dataframe.at[4763, "Poro"])
+
+    # corner cases
+    with pytest.raises(ValueError):
+        usewell.mask_shoulderbeds(
+            ["Zonelog", "Facies"], ["NOPerm", "Poro"], nsamples=2, strict=True
+        )
+
+    with pytest.raises(ValueError):
+        usewell.mask_shoulderbeds(["Perm", "Facies"], ["NOPerm", "Poro"], nsamples=2)
+
+    assert usewell.mask_shoulderbeds(["Zonelog"], ["Perm", "Poro"]) is True
+    assert usewell.mask_shoulderbeds(["Zonelog"], ["Dummy"]) is False
+    assert usewell.mask_shoulderbeds(["Dummy"], ["Perm", "Poro"]) is False
+
+
 def test_rescale_well(loadwell1):
     """Rescale (resample) a well to a finer increment"""