Merge pull request #2 from jcrivenaes/surf-slice-cube-window

Surf slice cube window

Author: jcrivenaes

src/xtgeo/cube/_cube_roxapi.py

@@ -61,6 +61,9 @@ def _roxapi_cube_to_xtgeo(self, rcube):
     self._xlines = np.array(range(xlstart, self._nrow + xlstart, xlincr),
                             dtype=np.int32)
 
+    # roxar API does not store traceid codes, assume 1
+    self._traceidcodes = np.ones((self._ncol, self._nrow), dtype=np.int32)
+
     if rcube.is_empty:
         xtg.warn('Cube has no data; assume 0')
     else:

src/xtgeo/grid3d/_grid_etc1.py

@@ -440,6 +440,58 @@ def collapse_inactive_cells(self):
     return self
 
 
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# Do cropping
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+def do_cropping(self, spec):
+    """Do cropping of geometry (and params)."""
+
+    (ic1, ic2), (jc1, jc2), (kc1, kc2) = spec
+
+    if ic1 < 1 or ic2 > self.ncol or jc1 < 1 or jc2 > self.nrow or \
+       kc1 < 1 or kc2 > self.nlay:
+
+        raise ValueError('Boundary for tuples not matching grid'
+                         'NCOL, NROW, NLAY')
+
+    # compute size of new cropped grid
+    nncol = ic2 - ic1 + 1
+    nnrow = jc2 - jc1 + 1
+    nnlay = kc2 - kc1 + 1
+
+    ntot = nncol * nnrow * nnlay
+    ncoord = (nncol + 1) * (nnrow + 1) * 2 * 3
+    nzcorn = nncol * nnrow * (nnlay + 1) * 4
+
+    new_num_act = _cxtgeo.new_intpointer()
+    new_p_coord_v = _cxtgeo.new_doublearray(ncoord)
+    new_p_zcorn_v = _cxtgeo.new_doublearray(nzcorn)
+    new_p_actnum_v = _cxtgeo.new_intarray(ntot)
+
+    _cxtgeo.grd3d_crop_geometry(self.ncol, self.nrow, self.nlay,
+                                self._p_coord_v,
+                                self._p_zcorn_v,
+                                self._p_actnum_v,
+                                new_p_coord_v,
+                                new_p_zcorn_v,
+                                new_p_actnum_v,
+                                ic1, ic2, jc1, jc2, kc1, kc2,
+                                new_num_act,
+                                0,
+                                xtg_verbose_level)
+
+    self._p_coord_v = new_p_coord_v
+    self._p_zcorn_v = new_p_zcorn_v
+    self._p_actnum_v = new_p_actnum_v
+
+    self._nactive = _cxtgeo.intpointer_value(new_num_act)
+    self._ncol = nncol
+    self._nrow = nnrow
+    self._nlay = nnlay
+
+    # TODO: subgrid
+
+
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 # Reduce grid to one layer
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

src/xtgeo/grid3d/grid.py

@@ -715,6 +715,25 @@ def collapse_inactive_cells(self):
 
         self = _grid_etc1.collapse_inactive_cells(self)
 
+    def do_cropping(self, croprange):
+        """Reduce the grid size by cropping, the grid will have new dimensions.
+
+        Args:
+            croprange (tuple): Crop range on the form
+                ((i1, i2), (j1, j2), (k1, k2)), from i1, i2 means from i1 to
+                i2, etc (inclusive).
+
+        Example::
+
+            >>> from xtgeo.grid3d import Grid
+            >>> gf = Grid('gullfaks2.roff')
+            >>> gf.do_cropping((3, 6), (4, 20), (1, 10))
+            >>> gf.to_file('gf_reduced.roff')
+
+        """
+
+        _grid_etc1.do_cropping(self, croprange)
+
     def reduce_to_one_layer(self):
         """Reduce the grid to one single single layer.
 

src/xtgeo/surface/_regsurf_cube.py

@@ -95,7 +95,7 @@ def slice_cube_window(self, cube, zsurf=None, other=None,
                       zrange=10, ndiv=None, attribute='max',
                       maskthreshold=0.1, snapxy=False,
                       showprogress=False, deadtraces=True,
-                      deletecube=False):
+                      deletecube=False, algorithm=1):
 
     """Slice Cube with a window and extract attribute(s)
 
@@ -140,12 +140,18 @@ def slice_cube_window(self, cube, zsurf=None, other=None,
     # This will run slice in a loop within a window. Then, numpy methods
     # are applied to get the attributes
 
-    if other is None:
+    if other is None and algorithm == 1:
         attvalues = _slice_constant_window(this, cube, sampling, zrange,
                                            ndiv, mask, attrlist, snapxy,
                                            showprogress=showprogress,
                                            deadtraces=deadtraces,
                                            deletecube=deletecube)
+    elif other is None and algorithm == 2:
+        attvalues = _slice_constant_window2(this, cube, sampling, zrange,
+                                            ndiv, mask, attrlist, snapxy,
+                                            showprogress=showprogress,
+                                            deadtraces=deadtraces,
+                                            deletecube=deletecube)
     else:
         attvalues = _slice_between_surfaces(this, cube, sampling, other,
                                             other_position, zrange,
@@ -157,19 +163,20 @@ def slice_cube_window(self, cube, zsurf=None, other=None,
 
     results = dict()
 
-    for attr in attrlist:
-        scopy = self.copy()
-        scopy.values = attvalues[attr]
-        results[attr] = scopy
-
-    # for backward compatibility
-    if qattr_is_string:
-        self.values = attvalues[attrlist[0]]
-        return None
+    if algorithm != 2:
+        for attr in attrlist:
+            scopy = self.copy()
+            scopy.values = attvalues[attr]
+            results[attr] = scopy
 
-    return results
+        # for backward compatibility
+        if qattr_is_string:
+            self.values = attvalues[attrlist[0]]
+            return None
 
-    logger.info('Mean of cube attribute is {}'.format(self.values.mean()))
+        return results
+    else:
+        return None
 
 
 def _slice_constant_window(this, cube, sampling, zrange,
@@ -224,6 +231,77 @@ def _slice_constant_window(this, cube, sampling, zrange,
     return attvalues  # this is dict with numpies, one per attribute
 
 
+def _slice_constant_window2(this, cube, sampling, zrange,
+                            ndiv, mask, attrlist, snapxy, showprogress=False,
+                            deadtraces=True, deletecube=False):
+    """Slice a window, (constant in vertical extent); faster and better
+    algorithm (algorithm2).
+    """
+
+    zincr = zrange / float(ndiv)
+    ztmp = this.copy()
+    ztmp.values = ztmp.values - zincr * (ndiv + 1)
+
+    if mask:
+        opt2 = 0
+    else:
+        opt2 = 1
+
+    if deadtraces:
+        # set dead traces to cxtgeo UNDEF -> special treatment in the C code
+        olddead = cube.values_dead_traces(_cxtgeo.UNDEF)
+
+    cubeval1d = np.ravel(cube.values, order='C')
+
+    usesampling = 0
+    if sampling == 'trilinear':
+        usesampling = 1
+        if snapxy:
+            usesampling = 2
+
+    # allocate the attribute maps
+    nattr = 5  # predefined attributes
+    nsurf = ztmp.ncol * ztmp.nrow * nattr
+
+    istat, attrmaps = _cxtgeo.surf_slice_cube_window(
+        cube.ncol,
+        cube.nrow,
+        cube.nlay,
+        cube.xori,
+        cube.xinc,
+        cube.yori,
+        cube.yinc,
+        cube.zori,
+        cube.zinc,
+        cube.rotation,
+        cube.yflip,
+        cubeval1d,
+        ztmp.ncol,
+        ztmp.nrow,
+        ztmp.xori,
+        ztmp.xinc,
+        ztmp.yori,
+        ztmp.yinc,
+        ztmp.yflip,
+        ztmp.rotation,
+        ztmp.get_values1d(),
+        zincr,
+        ndiv * 2,
+        nsurf,
+        nattr,
+        usesampling,
+        opt2,
+        xtg_verbose_level)
+
+    if istat != 0:
+        logger.warning('Problem, ISTAT = {}'.format(istat))
+
+    if deadtraces:
+        cube.values_dead_traces(olddead)  # reset value for dead traces
+
+    return istat
+
+
 def _slice_between_surfaces(this, cube, sampling, other, other_position,
                             zrange, ndiv, mask, attrlist, mthreshold,
                             snapxy, showprogress=False, deadtraces=True,

src/xtgeo/surface/regular_surface.py

@@ -1325,7 +1325,7 @@ 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):
+                          deadtraces=True, algorithm=1):
         """Slice the cube within a vertical window and get the statistical
         attrubute.
 
@@ -1374,6 +1374,7 @@ def slice_cube_window(self, cube, zsurf=None, other=None,
             deadtraces (bool): 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 old method, 2 for new alternative.
 
         Example::
 
@@ -1411,7 +1412,8 @@ def slice_cube_window(self, cube, zsurf=None, other=None,
                                                  maskthreshold=maskthreshold,
                                                  snapxy=snapxy,
                                                  showprogress=showprogress,
-                                                 deadtraces=deadtraces)
+                                                 deadtraces=deadtraces,
+                                                 algorithm=algorithm)
         return asurfs
 
     # =========================================================================

tests/test_grid3d/test_grid_operations.py

@@ -131,6 +131,7 @@ def test_refine_vertically():
     grd.to_file('TMP/test_refined_by_3.roff')
 
 
+
 def test_refine_vertically_per_zone():
     """Do a grid refinement vertically, via a dict per zone."""
 
@@ -151,3 +152,18 @@ def test_refine_vertically_per_zone():
     grd.refine_vertically(refinement, zoneprop=zone)
 
     grd.to_file('TMP/test_refined_by_dict.roff')
+
+
+def test_crop_grid():
+    """Crop a grid."""
+
+    logger.info('Read grid...')
+
+    grd = Grid(emegfile)
+    logger.info('Read grid... done, NLAY is {}'.format(grd.nlay))
+    logger.info('Read grid...NCOL, NROW, NLAY is {} {} {}'
+                .format(grd.ncol, grd.nrow, grd.nlay))
+
+    grd.do_cropping(((30, 60), (20, 40), (1, 46)))
+
+    grd.to_file('TMP/grid_cropped.roff')

tests/test_surface/test_regular_surface_vs_cube.py

@@ -186,6 +186,48 @@ def test_slice_attr_window_max():
                   infotext='Method: trilinear, window')
 
 
+@tsetup.skipsegyio
+@tsetup.skipifroxar
+def test_slice_attr_window_max_algorithm2():
+    """Slice a cube within a window, get max, using trilinear interpol,
+    alg 2.
+    """
+
+    logger.info('Loading surface')
+    xs1 = RegularSurface(rtop1)
+    xs2 = xs1.copy()
+    xs3 = xs1.copy()
+
+    logger.info('Loading cube')
+    cc = Cube(rsgy1)
+
+    t1 = xtg.timer()
+    xs1.slice_cube_window(cc, attribute='min', sampling='trilinear',
+                          algorithm=2)
+    t2 = xtg.timer(t1)
+    # logger.info('Window slicing... {} secs'. format(t2))
+
+    # xs1.quickplot(filename=td + '/surf_slice_cube_window_min.png',
+    #               colortable='seismic', minmax=(-1, 1),
+    #               title='Reek Minimum',
+    #               infotext='Method: trilinear, window')
+
+    # xs2.slice_cube_window(cc, attribute='max', sampling='trilinear',
+    #                       showprogress=True)
+
+    # xs2.quickplot(filename=td + '/surf_slice_cube_window_max.png',
+    #               colortable='seismic', minmax=(-1, 1),
+    #               title='Reek Maximum',
+    #               infotext='Method: trilinear, window')
+
+    # xs3.slice_cube_window(cc, attribute='rms', sampling='trilinear')
+
+    # xs3.quickplot(filename=td + '/surf_slice_cube_window_rms.png',
+    #               colortable='jet', minmax=(0, 1),
+    #               title='Reek rms (root mean square)',
+    #               infotext='Method: trilinear, window')
+
+
 @tsetup.skipifroxar
 def test_cube_attr_mean_two_surfaces():
     """Get cube attribute (mean) between two surfaces."""