BUG: fix compute_attributes_in_window() when input cube has trace_id_codes=2 present

The trace_id_codes=2 means "dead traces"

Author: jcrivenaes

src/xtgeo/cube/_cube_window_attributes.py

@@ -63,6 +63,7 @@ class CubeAttrs:
     _upper: RegularSurface | None = None  # upper surf, resampled to cube map resolution
     _lower: RegularSurface | None = None  # lower surf, resampled to cube map resolution
     _min_thickness_mask: RegularSurface | None = None  # mask for min. thickness trunc.
+    _mask_map_by_traceidcode: RegularSurface | None = None  # mask for traceidcode 2
 
     _result_attr_maps: dict = field(default_factory=dict)  # holds the resulting maps
 
@@ -94,6 +95,7 @@ def _process_upper_lower_surface(self) -> None:
             else self.lower_surface
         )
 
+        # the template surface is the topology that defines the resulting attribute maps
         self._template_surface = (
             upper
             if isinstance(self.upper_surface, (float, int))
@@ -184,7 +186,7 @@ def _create_reduced_cube(self) -> None:
         from xtgeo import Cube  # avoid circular import by having this here
 
         cubev = self.cube.values.copy()  # copy, so we don't change the input instance
-        cubev[self.cube.traceidcodes == 2] = np.nan  # set dead traces to nan
+        cubev[self.cube.traceidcodes == 2] = 0.0  # set traceidcode 2 to zero
 
         # Create a boolean mask based on the threshold
         mask = self._depth_array < self._outside_depth
@@ -311,13 +313,20 @@ def _add_to_attribute_map(self, attr_name: str, values: np.ndarray) -> None:
         attr_map = self._upper.copy()
         attr_map.values = np.ma.masked_invalid(values)
 
+        # apply mask for the cube's dead traces (traceidcode 2)
+        attr_map.values = np.ma.masked_where(
+            self.cube.traceidcodes == 2, attr_map.values
+        )
+
         # now resample to the original input map
         attr_map_resampled = self._template_surface.copy()
         attr_map_resampled.resample(attr_map)
 
-        attr_map_resampled.values = np.ma.masked_where(
-            self.upper_surface.values.mask, attr_map_resampled.values
-        )
+        # Use template_surface consistently for masking (it's already set correctly)
+        if hasattr(self._template_surface.values, "mask"):
+            attr_map_resampled.values = np.ma.masked_where(
+                self._template_surface.values.mask, attr_map_resampled.values
+            )
 
         self._result_attr_maps[attr_name] = attr_map_resampled
 

tests/test_surface/test_surface_cube_attrs.py

@@ -1,6 +1,7 @@
 import pathlib
 import warnings
 
+import numpy as np
 import pytest
 
 import xtgeo
@@ -11,6 +12,7 @@
 
 SFILE1 = pathlib.Path("cubes/etc/ib_synth_iainb.segy")
 SFILE2 = pathlib.Path("cubes/reek/syntseis_20030101_seismic_depth_stack.segy")
+SFILE3 = pathlib.Path("cubes/etc/cube_w_deadtraces.segy")
 
 TOP2A = pathlib.Path("surfaces/reek/2/01_topreek_rota.gri")
 TOP2B = pathlib.Path("surfaces/reek/2/04_basereek_rota.gri")
@@ -35,6 +37,13 @@ def fixture_loadsfile2(testdata_path):
     return xtgeo.cube_from_file(testdata_path / SFILE2)
 
 
+@pytest.fixture(name="loadsfile3")
+def fixture_loadsfile3(testdata_path):
+    """Fixture for loading a SFILE3"""
+    logger.info("Load seismic file 3")
+    return xtgeo.cube_from_file(testdata_path / SFILE3)
+
+
 def test_single_slice_yflip_snapxy_both(loadsfile1):
     """Test cube values with exact location of maps, compared to cube samples"""
     cube1 = loadsfile1
@@ -517,3 +526,41 @@ def test_warn_errs_new_module_reek(loadsfile1, loadsfile2, testdata_path):
 
     with pytest.raises(ValueError, match="no valid data in the interval"):
         _ = cube2.compute_attributes_in_window(surf1, surf1, ndiv=4)  # same surface
+
+
+def test_compute_attrs_handling_of_dead_traces(loadsfile3):
+    """Test handling of dead traces in the cube."""
+    cube = loadsfile3
+
+    # Create a surface that is above the dead traces
+    surf = xtgeo.surface_from_cube(cube, 1000.0)
+
+    # Compute attributes in window
+    attrs = cube.compute_attributes_in_window(surf, surf + 10)
+
+    # Check if the result is as expected
+    assert "mean" in attrs
+    assert attrs["mean"].values.mean() == pytest.approx(239.18202, abs=1e-4)
+
+    # count number of dead traces in cube
+    dead_traces = (cube.traceidcodes == 2).sum()
+
+    # count number of masked values in map
+    masked_values = np.ma.count_masked(attrs["mean"].values)
+
+    assert dead_traces == masked_values, (
+        "Number of dead traces does not match masked values"
+    )
+
+
+def test_attribute_around_constant_cube_slices(loadsfile2):
+    """Get attribute around a constant cube slices."""
+
+    mycube = loadsfile2
+
+    level1 = 1710
+    level2 = 1730
+
+    myattrs = mycube.compute_attributes_in_window(level1, level2)
+
+    assert myattrs["mean"].values.mean() == pytest.approx(-0.00293736, abs=1e-5)