Merge pull request #440 from jcrivenaes/fix-surf-cube-slice-issue-grane

Fix surf cube slice issue grane

src/clib/xtg/surf_cube_attr_intv.c

@@ -276,8 +276,9 @@ surf_cube_attr_intv(int ncol,
     int i, ic;
     for (i = 0; i < nsurf1; i++) {
         for (ic = 0; ic <= ndiv; ic++) {
-            if (maskv1[i] == 0 && maskv2[i] == 0) {
-                stack[i][ic] = surfsv1[i] + ic * slicezinc;
+            double vadd = surfsv1[i] + ic * slicezinc;
+            if (maskv1[i] == 0 && maskv2[i] == 0 && vadd <= surfsv2[i]) {
+                stack[i][ic] = vadd;
                 rmask[i][ic] = 0;
             } else {
                 stack[i][ic] = UNDEF;

src/clib/xtg/surf_stack_slice_cube.c

@@ -5,7 +5,7 @@
  *    surf_stack_slice_cube.c
  *
  * DESCRIPTION:
- *    Slicing a stack of surfaces and reurn cube values
+ *    Slicing a stack of surfaces and return cube values
  *
  *    In this version, the surface and cube share geometry, which makes calculations
  *    way much simpler and probably also much faster
@@ -14,10 +14,8 @@
  *    ncol, nrow...  i     Cube dimensions and relevant increments
  *    czori, czinc   i     Cube settings
  *    cubevalsv      i     Cube values
- *    stack          i     map array stack with Z values
- *    rmask          i     Length of slice array
- *    surfsv        i/o    map to update
- *    maskv         i/o    mask array for map, may be updated!
+ *    stack         i/o    stacked map array stack with Z values
+ *    rmask          i     stacked map array stack with masks
  *    optnearest     i     If 1 use nerest node, else do interpolation aka trilinear
  *    optmask        i     If 1 then masked cells (undef) are made if cube is UNDEF
  *                         (surface outside cube will not occur in this case)
@@ -68,6 +66,8 @@ surf_stack_slice_cube(int ncol,
                     continue;
 
                 double zval = stack[icmap][nd];
+                if (zval > UNDEF_LIMIT)
+                    continue;
 
                 // find vertical index of node right above
                 int k1 = (int)((zval - czori) / czinc);

src/xtgeo/grid3d/_grid_roxapi.py

@@ -349,6 +349,9 @@ def _export_grid_cornerpoint_roxapi_v1(
     geom = CPG.create(self.dimensions)
 
     logger.info(geom)
+    scopy = self.copy()
+    scopy.make_zconsistent()
+    scopy._xtgformat1()
 
     npill = (self.ncol + 1) * (self.nrow + 1) * 3
     nzcrn = (self.ncol + 1) * (self.nrow + 1) * 4 * (self.nlay + 1)
@@ -358,7 +361,7 @@ def _export_grid_cornerpoint_roxapi_v1(
         self.nrow,
         self.nlay,
         self._coordsv,
-        self._zcornsv,
+        scopy._zcornsv,
         self._actnumsv,
         npill,
         npill,
@@ -421,6 +424,7 @@ def _export_grid_cornerpoint_roxapi_v2(
 
     scopy = self.copy()
     scopy.make_zconsistent()
+    scopy._xtgformat2()
 
     npill = (self.ncol + 1) * (self.nrow + 1) * 3
     nzcrn = (self.ncol + 1) * (self.nrow + 1) * (self.nlay + 1) * 4

src/xtgeo/surface/_regsurf_cube_window_v2.py

@@ -160,7 +160,7 @@ def _slice_cube_window(
             ndiv = 1
             logger.warning("NDIV < 1; reset to 1")
 
-    # force/overrule a coarse sampling for sampling "cube"
+    # force/overrule a coarse sampling for sampling option "cube"
     ndivdisc = int((zdiameter) * 1.0001 / cube.zinc)
     if sampling == "cube":
         ndiv = ndivdisc
@@ -269,7 +269,7 @@ def _slice_cube_window_resample(
     showprogress,
     deadtraces,
 ):
-    """Makes a resample from origonal surfaces first to fit cube topology"""
+    """Makes a resample from original surfaces first to fit cube topology"""
 
     logger.info("Attributes between surfaces, resampling version")
 

tests/test_surface/test_surface_cube_attrs.py

@@ -15,7 +15,10 @@
 TPATH = xtg.testpath
 
 SFILE1 = join(TPATH, "cubes/etc/ib_synth_iainb.segy")
+SFILE2 = join(TPATH, "cubes/reek/syntseis_20030101_seismic_depth_stack.segy")
 
+TOP2A = join(TPATH, "surfaces/reek/2/01_topreek_rota.gri")
+TOP2B = join(TPATH, "surfaces/reek/2/04_basereek_rota.gri")
 
 # ======================================================================================
 # This is a a set of tests towards a synthetic small cube made by I Bush in order to
@@ -30,6 +33,13 @@ def fixture_loadsfile1():
     return xtgeo.Cube(SFILE1)
 
 
+@pytest.fixture(name="loadsfile2")
+def fixture_loadsfile2():
+    """Fixture for loading a SFILE2"""
+    logger.info("Load seismic file 2")
+    return xtgeo.Cube(SFILE2)
+
+
 def test_single_slice_yflip_snapxy_both(loadsfile1):
     """Test cube values with exact location of maps, compared to cube samples"""
     cube1 = loadsfile1
@@ -271,3 +281,30 @@ def test_avg_surface_large_cube():
     print("Algorithm 1: ", xtg.timer(t1))
 
     logger.info("Testing done")
+
+
+def test_attrs_reek(loadsfile2):
+
+    logger.info("Make cube...")
+    cube2 = loadsfile2
+
+    t2a = xtgeo.RegularSurface(TOP2A)
+    t2b = xtgeo.RegularSurface(TOP2B)
+
+    attlist = ["maxpos", "maxneg"]
+
+    attrs1 = t2a.slice_cube_window(
+        cube2, other=t2b, sampling="trilinear", attribute=attlist, algorithm=1
+    )
+    attrs2 = t2a.slice_cube_window(
+        cube2, other=t2b, sampling="trilinear", attribute=attlist, algorithm=2
+    )
+
+    for att in attrs1.keys():
+        srf1 = attrs1[att]
+        srf2 = attrs2[att]
+
+        srf1.to_file(join(TMD, "attr1_" + att + ".gri"))
+        srf2.to_file(join(TMD, "attr2_" + att + ".gri"))
+
+        assert srf1.values.mean() == pytest.approx(srf2.values.mean(), abs=0.005)